Practice 09 - BPR for SLIM and MF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems 2020/21\n",
    "\n",
    "### Practice - BPR for SLIM and MF\n",
    "\n",
    "### State of the art machine learning algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A few info about gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": [
    "### Gradient Descent"
   ]
  },
  {
   "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>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: x**3-2*x**2+2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVyVdd7/8dfnsIiACAi44S7u5oZ7ltli61hN+2LLNFZ3zdRMzdZMzXL3q2bmbpk2zZapqaYas8wpzWxxX0FFRQQRUJAdlFUQON/fHxzvmyGQAx64zvJ5Ph7n4YFzcXxzgW+v8z3f63uJMQallFLez2Z1AKWUUl1DC18ppXyEFr5SSvkILXyllPIRWvhKKeUjtPCVUspHtFn4IhIkIjtEJElEkkXkjy1sIyLyooiki8heEZncOXGVUkp1lL8T29QC84wxlSISAGwSkdXGmG1NtrkMiHPcpgOLHX8qpZRyE20e4ZtGlY4PAxy35mdrLQD+4dh2GxAuIn1dG1UppdTZcOYIHxHxAxKB4cArxpjtzTbpD2Q3+TjH8bm8Zs+zCFgEEBISMmXUqFEdjK2UUtYyBpLzyogMCaRfz+5d9vcmJiYWG2OiO/K1ThW+MaYBmCgi4cCnIjLOGLO/ySbS0pe18DxLgaUA8fHxJiEhoQORlVLKeruPHueaV7ew+NbJXDa+6wY0RORIR7+2XbN0jDEngHXApc0eygEGNPk4FsjtaCillHJ3CVnHAZgyOMLiJM5zZpZOtOPIHhHpDlwEHGy22UpgoWO2zgygzBiTh1JKeamdWaUM7hVMTI8gq6M4zZkhnb7AO45xfBvwL2PM5yJyH4AxZgmwCrgcSAeqgbs6Ka9SSlnOGEPCkePMGxVjdZR2abPwjTF7gUktfH5Jk/sGeMC10ZRSyj1lFFdRWnWKqR40nAN6pq1SSrVbQlYpAPGDIy1O0j5a+Eop1U47s44TGRLI0KgQq6O0ixa+Ukq1U0JWKfGDIhBpaUa6+9LCV0qpdiisqCGrpJqpHjacA1r4SinVLomO+ffxHvaGLWjhK6VUu+zMOk5QgI2x/XpaHaXdtPCVUqodEo6UMnFAOIH+nlefnpdYKaUsUlVbT3JuuUeO34MWvlJKOW330RM02A1TBnne+D1o4SullNO2Z5bgZxOPO+HqNC18pZRy0raMEsb370loN6dWlnc7WvhKKeWEk6ca2JN9ghlDe1kdpcO08JVSygm7jh6nrsEwY6hnDueAFr5SSjllW4Znj9+DFr5SSjnF08fvQQtfKaXa5A3j96CFr5RSbfKG8XvQwldKqTZ5w/g9aOErpVSbvGH8HrTwlVLqjLxl/B608JVS6oy8ZfwetPCVUuqMvGX8HrTwlVLqjLxl/B608JVSqlVVtfVeM34PWvhKKdWqHVml1DUYzh0eZXUUl9DCV0qpVmw+VEygv80jL1jeEi18pZRqxab0YqYOjiAowM/qKC7RZuGLyAAR+U5EUkQkWUQeamGbuSJSJiJ7HLcnOieuUkp1jaKKWg7mVzDbS4ZzAJx527keeMQYs0tEegCJIrLWGHOg2XYbjTFXuj6iUkp1vS2HiwGYMzza4iSu0+YRvjEmzxizy3G/AkgB+nd2MKWUstKmQ8WEBwcwpl+Y1VFcpl1j+CIyGJgEbG/h4ZkikiQiq0VkrAuyKaWUJYwxbE4vZtawXvjZxOo4LuN04YtIKLAceNgYU97s4V3AIGPMBOAlYEUrz7FIRBJEJKGoqKijmZVSqlNlFFeRW1bjVeP34GThi0gAjWX/vjHmk+aPG2PKjTGVjvurgAAR+d6eMsYsNcbEG2Pio6O9Z1xMKeVdNqc3jt97y/z705yZpSPAm0CKMea5Vrbp49gOEZnmeN4SVwZVSqmusulQMbER3RkYGWx1FJdyZpbObOB2YJ+I7HF87jFgIIAxZglwHXC/iNQDJ4GbjDGmE/IqpVSnqm+wszWjhCvG98VxHOs12ix8Y8wm4IzftTHmZeBlV4VSSimr7D1WRkVNvdeN34OeaauUUv9hQ1oRImjhK6WUt1uXWsSE2HAiQwKtjuJyWvhKKeVQWnWKpJwTnD/CO2cRauErpZTDxkNFGANzR2rhK6WUV1ufWkREcADnxIZbHaVTaOErpRRgtxs2HCpiTly0Vy2n0JQWvlJKAcm55RRXnvLa4RzQwldKKQDWpxUCMCdOC18ppbzautQixvfvSXSPblZH6TRa+Eopn1dWXceuo8e9djrmaVr4Simftym9GLsXT8c8TQtfKeXz1qUWEhbkz8QB3jkd8zQtfKWUT7PbDd+lFnL+yBj8/by7Er37u1NKqTbsyTlBceUpLhodY3WUTqeFr5TyaV8fKMDPJswdoYWvlFJe7ZuUQqYOjqBncIDVUTqdFr5Symdll1aTWlDBRaN7Wx2lS2jhK6V81tcpBQBa+Eop5e2+TilgeEwog6NCrI7SJbTwlVI+qbymju0ZpVzoA7NzTtPCV0r5pPWpRdTbjc8M54AWvlLKR32TUkBEcACTB0ZYHaXLaOErpXxOXYOd71KLuGBUjNde7KQlWvhKKZ+zPaOUspN1XDKmj9VRupQWvlLK56zen0f3AD+vXw65OS18pZRPabAb1iQXcMGoaLoH+lkdp0tp4SulfErikeMUV9Zy6bi+Vkfpclr4Simf8uX+fAL9bMwb5Tvz709rs/BFZICIfCciKSKSLCIPtbCNiMiLIpIuIntFZHLnxFVKqY4zxrAmOZ85cVGEdvO3Ok6Xc+YIvx54xBgzGpgBPCAiY5ptcxkQ57gtAha7NKVSSrnA3pwyjp04yaXjfGt2zmlt/hdnjMkD8hz3K0QkBegPHGiy2QLgH8YYA2wTkXAR6ev4WmWBytp69h8r42BeOakFFeSX1VBUWcvxqjqMMRigm7+NiJBAIoMDGRAZzLCYUEbEhHJObLjPvZmlfMPq/fn424SLx/jO2bVNtes1jYgMBiYB25s91B/IbvJxjuNz/1H4IrKIxlcADBw4sH1JVZsyi6tYtS+P9WlF7DpynHq7ASAiOIDYiGCiQ7sxIqYHNpsgQE29nRPVp8gtq2FrRgnVpxoA8LcJ42N7Mn1ILy4eE8OkARHYfOjkFOWdjDF8uT+PmcN6ER4caHUcSzhd+CISCiwHHjbGlDd/uIUvMd/7hDFLgaUA8fHx33tctV9NXQMrdh9jWWIOiUeOAzCufxg/Pm8o04ZEMrZvGNE9uiFy5sI2xpBXVsPB/HJ2Zh1nZ2Ypb27KYMn6w0T36MalY/tw3ZRYzont2eZzKeWOUgsqyCqp5sfnDbU6imWcKnwRCaCx7N83xnzSwiY5wIAmH8cCuWcfT7WmoqaO97cf5Y2NmRRX1hIXE8pvLhvF1ZP60zssqN3PJyL0C+9Ov/DuzBvV+HK3vKaO7w4WsiY5n2WJ2by77Qhj+oZx8/SB/HByf4IDfe9NL+W5Pk/Kwyb43Nm1TbX5L1YaD+feBFKMMc+1stlK4EER+RCYDpTp+H3naLAb/pWQzf+sSaWk6hRz4qK4f+5EZg7t5fIj77CgABZM7M+Cif0pr6njs93H+OeObB5fsZ/n16Zx16zBLJw52CcuDac8mzGGlUm5zBoWRXSPblbHsYwzh2izgduBfSKyx/G5x4CBAMaYJcAq4HIgHagG7nJ9VJWUfYLHPt1Hcm450wZH8uado5k4ILxL/u6woABunzmY22YMIuHIcRavO8yza9NYsv4wd587hEXnDaVHkBa/ck9JOWUcLa3mwQuGWx3FUs7M0tlEy2P0TbcxwAOuCqX+U12DnZe+TeeV79KJDu3GSzdP4spz+loyli4iTB0cydQ7I0nJK+fl79J56dt0/rn9KA9dFMfN0wYS4Kfn8yn3snJPLoF+Nub76HTM0/RfppvLLq3m2le38OI3h1gwsR9rfnYeV03o5xZvnI7uG8Yrt0zmswdmMzwmlCc+S2b+CxvYcrjY6mhK/a8Gu+HzvbnMHRlNz+6+/SpUC9+NbTpUzFUvb+JISRWLb53MczdMdMtf2AkDwvlw0QzevCOe+gbDLa9v5+cf7aG4stbqaEqxPbOEwopafjCxn9VRLKfTLNzUGxszeGpVCsNjQll6e7zbX2RZRLhwdG9mD4/ile/SWbL+MN8cLOTxK8fww8n93eIVifJN/07KJSTQjwtH+ebJVk3pEb6bsdsNT61K4ckvUrhkTB8+/a/Zbl/2TQUF+PHIJSNZ/dAcRvbuwaPLkrj33UQ92leWOFVvZ9W+fC4e01vPHkcL363UNdj5xcd7WbohgztmDuLVWycT4qELPA2P6cGHi2bwuytGsy6tiPnPb+Cr5HyrYykfs/FQEWUn63Q4x0EL303UN9h56MPdLN+Vw88vHsEffjDW45czsNmEe+YM5d8PnkvvsCAWvZvI71bso6auwepoykd8tieX8OAAzh3uW1e2ao0WvhtosBseXZbEqn35/O6K0fz0wjivGvMe2acHKx6YzaLzhvLetqNct2QLR0qqrI6lvFx5TR1fHcjnivF9CfTXqgMtfMvZ7YbHPtnHij25/GL+SO6Z453rfAT623js8tG8sTCe7NKTXPniJr7crydjq86zam8eNXV2ro8f0PbGPkIL32LPrk3lo4RsfjJvOA/4wFmAF43pzec/OZehMaHc994unl6dQoNd19FTrvdxYg7DY0KZENvT6ihuQwvfQv/amc0r3x3m5mkD+PnFI6yO02UGRAaz7N6Z3Dp9IK+tz+Ced3ZSXlNndSzlRTKLq0g4cpwfTo71quHRs6WFb5FNh4p57NN9zImL4k8LxvncL2Wgv43/d814nrx6HBsPFXP1K5vJKKq0OpbyEp/sysEmcM2k/lZHcSta+BY4UlLF/e8nMjwmlFdvnezTa8/cNmMQ798znRPVdSx4ZTPr04qsjqQ8nN1uWJ6Yw5y4aPr0bP9S4d7Md5vGIidPNXDvu4n42YTXF8brCpPA9KG9WPngbGIjgrn77Z18uOOo1ZGUB9uaUUJuWQ3XTYm1Oorb0cLvQsYYfvvpPlILKnjhxokMiAy2OpLbiI0IZtl9M5k9PIpff7KPZ79KpXERVqXa5+PEHHoE+fvsdWvPRAu/C32wI5tPdh/j4QtHMHdkjNVx3E5oN3/evCOeG+MH8NK36TzyryRO1dutjqU8SHlNHV/uz+eqCf0ICtClFJrzzPP2PVB6YSV/+jyZOXFR/GSe90+/7KgAPxvP/HA8sRHdeXZtGvnlNSy+bYpbrhKq3M9nu49xsq6BG3XufYv0CL8LnKq38/BHuwkO9OfZ6yd4/JIJnU1E+MmFcTx3wwR2ZJZy42tbKayosTqWcnPGGN7ffpSx/cI4R+fet0gLvws8tzaN/cfKeeba8cR04ALjvuraybH8/a6pHC2t5volW8kurbY6knJju7NPcDC/glumD/S5ac7O0sLvZAlZpby2ofHkqkvG+vbl1TpiTlw07zmmbV63ZAtpBRVWR1Ju6p/bjxIS6MeCiTr3vjVa+J2opq6BXy3fS7+e3fndFWOsjuOxJg+M4KN7Z2A3cMNrW9mTfcLqSMrNlJ2s4/O9uSyY1J9QD11SvCto4XeiV75L53BRFU9fO95j17V3F6P6hLH8vln0CPLn1te3sSVdr5ur/s+nu3KoqbNzy7SBVkdxa1r4nSQlr5zF6w5z7eT+nDdC1+J2hYG9gvn4vlnERgRz5993svZAgdWRlBswxvDPHUeZENuTcf31zdoz0cLvBA12w6+X76Vn9wAe16Ecl+odFsRH985gdL8w7n8vkdX7dIllX7cto5S0gkpunT7I6ihuTwu/E/xz+xGScsp44qoxRIQEWh3H64QHB/Luj6YxYUA4D36wm38n5VodSVno75sziQgO0MsYOkEL38WOV53if75KY9awXvxggv4CdpawoADeuXsaUwZF8NCHu/lkV47VkZQFskurWZtSwC3TB+qZtU7Qwnex59amUVlbz++vGqtzgTtZaDd/3r5rKjOG9uKRZUn8a2e21ZFUF3tnSxZ+Itw+Y7DVUTyCFr4LHcgt5/3tR7h9xiBG9ulhdRyfEBzoz1t3TuXc4VH8cvle3t9+xOpIqotU1dbzUUI2l43vq8sgO6nNwheRt0SkUET2t/L4XBEpE5E9jtsTro/p/owx/OHfyfTsHsDPLvKdq1e5g6AAP15fGM+8UTH89tP9vL050+pIqgss35VDRU09d80ebHUUj+HMEf7bwKVtbLPRGDPRcfvT2cfyPGuS89mRWcqj80fSM1gX+upqQQF+LLltCpeM6c0f/n2ANzZmWB1JdSK73fD2liwmDAhn8sAIq+N4jDYL3xizASjtgiweq77Bzl/WpBIXE8pNU/XED6sE+tt45dbJXD6+D09+kcLSDYetjqQ6yTcHC8koquJuPbpvF1eN4c8UkSQRWS0iY1vbSEQWiUiCiCQUFXnPpeyWJeaQUVTFL+aPxE9XwrRUgJ+Nv900iSvO6ctTqw6yeJ2WvrcxxvDqunRiI7pzxfi+VsfxKK44338XMMgYUykilwMrgLiWNjTGLAWWAsTHx3vF5YxOnmrgha/TmDIoQq+w4yYC/Gz87caJ+Inw5y8PYjeGBy7QaxB4ix2Zpew+eoL/XjAWfx++HnRHnHXhG2PKm9xfJSKvikiUMcYnFjt5e0sWBeW1vHTzZJ2G6Ub8/Ww8d8MEbAJ/XZOK3W74yYUtHocoD7N4/WF6hQRyvV7kpN3OuvBFpA9QYIwxIjKNxmGikrNO5gHKqutYvC6deaNimDYk0uo4qhl/PxvP3jARmwjPrk2jwRge1hlUHu1AbjnrUot49JIReqJVB7RZ+CLyATAXiBKRHOD3QACAMWYJcB1wv4jUAyeBm4yPXH36zc2ZlNfU84v5I62OolrhZxP+6rjK2AtfH8Ju4GcXxemrMQ+1ZP1hQgL99ESrDmqz8I0xN7fx+MvAyy5L5CHKTtbx982ZXDauD6P7hlkdR52Bn034yw/PwSbw4jeHMMbw84tHaOl7mKziKj7fm8uPzh2iU587SBdp76B3tmRRUVPPg3pBco9gswnPXHsONhFe+jadBrvhF/NHaul7kBe/PUSgv40fnzfU6igeSwu/Aypq6nhzUyYXj+nN2H66/ransNmEp64Zj80mvLruMHYDv7pUS98THC6qZMXuY9wzZygxPXQZhY7Swu+Af2w9QtnJOn46T2d9eBqbTXhywThs0jgebDeG31w2Skvfzb34zSGCAvy4V4/uz4oWfjtV1dbzxsYMLhgZzfhYPbr3RDab8N8LxmETYemGDBrsht9dMVpL300dKqhgZVIu950/jF6h3ayO49G08Nvpgx1HOV5dp3O6PZyI8McfjMUmwpubMrEbwxNXjtHSd0MvfHOI4AA/Fs3Ro/uzpYXfDnUNdv6+OYtpQyJ1wSYvICL8/qox2ER4a3MmxsDvr9LSdyf7j5Xxxd48HrxguF49zgW08Nth1b48jp04yR9/0OpyQcrDiAiPXzkaPxu8vjGTBrtpPPLXNZEsZ4zhqVUpRAQHsOh8Pbp3BS18JxljWLohg2HRIcwbFWN1HOVCIsJjl4/GZhNeW5+B3ZjGMX4tfUutSy1iy+ES/nDVGMKCdN69K2jhO2nr4RKSc8t55trxWgReSET49aWjsImweF3j7J3/d7X+rK1S32DnqVUpDIkK4Zbpg6yO4zW08J302oYMokK7cfWk/lZHUZ1ERPjl/JH4ifDyd+nY7fC0/gdviWWJORwqrGTJbZMJ9NcVMV1FC98JqfkVrE/TBZt8gYjwyCUjsNmEF785RIMx/PmH5+h1DrpQZW09z61NI35QBPPH9rE6jlfRwnfCO1uz6OZv41Z9aekTRISfXzwCm+BYcM3w1+smaOl3kb99nUZxZS2vL4zXGVMupoXfhrKTdXy66xgLJvbTaWE+5uGLRmAT4bm1aRgD/3O9ln5nS82v4K3NWdw0dQATB4RbHcfraOG3YXliDifrGlg4c7DVUZQFfnphXOMSy2tSsRvDs9dP0KssdRJjDI+v2E9YkD+/nD/K6jheSQv/DOx2w7vbjjB5YDjj+usyCr7qgQuGIwJ/+TIVu4Hnb9DS7wyf7DrGjqxSnrl2vL6a7iRa+GewKb2YzOIqHrpxotVRlMX+a+5w/ER4evVB7HbDCzdNJEBL32XKqut4enUKkwaGc4NeurDTaOGfwT+2ZhEVGshl43WmgIJ7zx+Gn0148osU7Mbw4s2TtPRd5I+fJ3O8uo537tYT3jqT/ra2Iru0mm8OFnLT1IF089epmKrRPXOG8viVY1i9P58H/7mLU/V2qyN5vG9SCvhk1zEemDtMry/RybTwW/H+9qPYRLhl+kCroyg386Nzh/D7q8awJrmARe8mUH2q3upIHqusuo7HPt3HqD49eFCvL9HptPBbUNdg5+PEbOaNiqFfeHer4yg3dNfsITx97Xg2pBVxy+vbKa06ZXUkj/Snzw9QXHmKv143Qc+o7QK6h1vwTUohxZWnuGmqvnmkWnfztIEsvm0KKXnlXLdkC9ml1VZH8igrk3JZviuHB+YO04sJdREt/Bb8KyGbmB7dOH9EtNVRlJubP7YP790zneKKWn64eAspeeVWR/IIR0uq+e0n+5g8MJyf6sWEuowWfjP5ZTWsSy3k+vhYnWutnDJ1cCQf3z8LP5tww5KtbD1cYnUkt1bXYOenH+4Ggb/dNEn/nXUh3dPNLN+Vg92gc4FVu4zo3YPl98+iT88g7nhrB8sTc6yO5Lb+8uVB9mSf4Jlrz2FAZLDVcXyKFn4Tdrvho53ZzBgayaBeIVbHUR6mX3h3Pr5vFlOHRPDIsiT+8mXjSVrq/3y25xivb8xk4cxBXHFOX6vj+Bwt/Ca2ZZZwtLSaG/XNWtVBPYMDePuuadw8bSCvrjvMf72/S6dtOuw/Vsavlu9l2pBIHr9yjNVxfJIWfhP/2plNjyB/LhunRx6q4wL8bDx1zTh+d8Vo1hzI58bXtpFfVmN1LEsVVdRy77uJRAQH8uqtk/UMZYu0uddF5C0RKRSR/a08LiLyooiki8heEZns+pidr7ymjtX781kwsZ9e5ESdNRHhnjlDeWNhPBlFlVz50ia2Z/jmm7lVtfXc/fZOSqtO8drtU4gK7WZ1JJ/lzH+zbwOXnuHxy4A4x20RsPjsY3W9L/flU1tv57opOpyjXOfC0b359IHZhAX5c8sb23ljYwbG+M64fl2Dnfvf38WBvHJeuXUS58TqGvdWarPwjTEbgNIzbLIA+IdptA0IFxGPGxP5dPcxhkSFMEFPAFEuNqJ3Dz57cDYXjY7hyS9SePCfu6ms9f5xfbvd8KuP97IhrYinrhnHvFG9rY7k81wxkNYfyG7ycY7jc98jIotEJEFEEoqKilzwV7tG7omTbMss4eqJ/fWSaqpT9AgKYMltU/jNZaNYvT+PBS9vIjm3zOpYncZuN/xq+V4+2X2MRy4ewY1TdU0qd+CKwm+pIVt8zWqMWWqMiTfGxEdHu89ZrCuTcjEGFkzsZ3UU5cVEhHvPH8Z790ynoqaea17ZwhsbM7xu6ubpsl+WmMNDF8bxEz2T1m24ovBzgKYD37FArguet8us2H2MSQPDGRylc+9V55s1LIovHz6P80dG8+QXKSx8awcF5d4xi6e2voGHP9rzv2X/s4tHWB1JNeGKwl8JLHTM1pkBlBlj8lzwvF0iJa+cg/kVXDOpxVEopTpFZEggS2+fwlPXjCfhSCmXPL+BjxNzPPoN3bKTddzx1g5WJuXyy0tHatm7oTaveCUiHwBzgSgRyQF+DwQAGGOWAKuAy4F0oBq4q7PCdoYVe47hbxOuGO9x7zMrDyeO6y1MHxrJLz/ey6PLkvhszzGeuma8xy05cKSkih//I4HM4ipeuHEiV+sBlFsSq44o4uPjTUJCgiV/92l2u2HWM98ytl8Yb9451dIsyrfZ7Yb3th/hz6sPYjfw8EVx3Dl7sEdcbW1Ncj6PLktCgCW3TWHW8CirI3k1EUk0xsR35Gt9+nS3bZkl5JfX6NGIspzNJiycOZi1Pz+f2cN78fTqg1zy/AbWJOe77TBPTV0DT35+gHvfTWRIVAhf/HSOlr2b8+nC/3xvHsGBflw0WucHK/fQL7w7b9wxlXfunkagn417303k5te3sevocauj/YeErFIuf3Ejb2zK5LYZA1l230yPG4byRW2O4Xur+gY7X+7P58LRveke6P4vm5VvOX9ENLMfmsMHO7N5fm0a1766hTlxUfxkXhzThkRalquwoobn1x7iw51H6dezO+/+aBpz4txnirU6M58t/O2ZpZRWneKK8X2sjqJUi/z9bNw+YxDXTurP+9uPsHRDBje8tpX4QREsnDWYS8f26bLrwJZV1/HW5kxe35jBqXo7d84azKOXjCSkm89WiEfy2Z/W6eGcuSNjrI6i1BmFdPNn0XnDuH3GYD7YcZS3t2Tx0w92ExXajeumxHLVhL6M6RvWKWeJZxRV8s6WLJYl5lB9qoErxvflF/NH6jkrHsonC7++wc6a5MbhHF0ZU3mK7oF+3H3uEO6cNZj1h4p4b+sRXt+YwZL1hxkaHcL8sX2YPSyK+MERHf69NsZwpKSar1MK+HdSLkk5ZQT4CT+Y0J8fnTuEMf3CXPxdqa7kk4W/LeP0cI7OvVeex2YTLhgZwwUjYyitOsXq/Xl8npTH6xsyWLzuMIH+Ns7p35PRfcMY1bcHAyKC6dMziKjQbgT62/C3CQ12Q9nJOspO1nGkpJqM4koO5lWwI7OUfMdZv+P6h/Gby0ZxzaT+xIQFWfxdK1fwycL/Yl8uIYF+zB2pbzYpzxYZEsit0wdx6/RBVNbWsyOzhM3pJezNOcGK3ceo2Ob8qpx9woKIHxzBjKG9OHd4lA7beCGfK/yms3N0OEd5k9Bu/swb1ft/lyE2xnDsxElyT9RQUF5DcWUtdQ126hoMfjahZ/cAenYPIDaiO0OiQugRFGDxd6A6m88V/taMEo5X1+kFlJXXExFiI4KJjdD58aqRz514tWpfHiGBfpw/QodzlFK+xacKv8FuWHuggAtGxehwjlLK5/hU4e8+epziylPMH6snWymlfI9PFf5XBwoI8BOdnaOU8kk+U/jGGNYk5zNrWLGvuYwAAAkQSURBVJTORlBK+SSfKfxDhZUcKanmkrG6MqZSyjf5TOF/lZwPwMW6FLJSykf5TOGvPVDApIHheoq4Uspn+UTh55WdJCmnjEvG6OwcpZTv8onC//pAAQAXj9HhHKWU7/KJwv/qQAFDo0MYHhNqdRSllLKM1xd+eU0dWw+X6NG9UsrneX3hb0wrpt5udHaOUsrneX3hf3uwkJ7dA5g4INzqKEopZSmvLny73bA+rZDzR0Tj7+fV36pSSrXJq1tw37EyiitPMW+UXqhcKaW8uvC/PViICLr2vVJK4WThi8ilIpIqIuki8usWHp8rImUissdxe8L1Udvvu9RCJg0IJyIk0OooSilluTYvcSgifsArwMVADrBTRFYaYw4023SjMebKTsjYIUUVtezNKePRS0ZYHUUppdyCM0f404B0Y0yGMeYU8CGwoHNjnb11qYUAXKDj90opBThX+P2B7CYf5zg+19xMEUkSkdUiMtYl6c7Cd6mF9A7rxpi+YVZHUUopt+BM4UsLnzPNPt4FDDLGTABeAla0+EQii0QkQUQSioqK2pe0Heoa7GxMK+aCkTGItBRfKaV8jzOFnwMMaPJxLJDbdANjTLkxptJxfxUQICJRzZ/IGLPUGBNvjImPju68mTMJWcepqK3X4RyllGrCmcLfCcSJyBARCQRuAlY23UBE+ojjUFpEpjmet8TVYZ21Pq2IAD9h9vDv/Z+jlFI+q81ZOsaYehF5EFgD+AFvGWOSReQ+x+NLgOuA+0WkHjgJ3GSMaT7s02U2Hipi8sAIQru1+e0ppZTPcKoRHcM0q5p9bkmT+y8DL7s2WscUV9aSnFvOL+aPtDqKUkq5Fa8703ZzejEAc+J0OEcppZryusLfkFZMRHAAY/v1tDqKUkq5Fa8qfGMMGw8VMXt4FH42nY6plFJNeVXhpxVUUlhRy3lxuliaUko151WFv/FQ48lc5+r4vVJKfY9XFf6GQ8UMjwmlX3h3q6MopZTb8ZrCr6lrYHtGic7OUUqpVnhN4SdkHae23q7j90op1QqvKfyN6Y3LKUwfGml1FKWUckveU/hpxcQPiiQ4UJdTUEqplnhF4RdX1nIgr1xn5yil1Bl4ReFvPdy4MKeujqmUUq3zjsLPKKFHN3/G9dOrWymlVGu8o/APlzBtSCT+fl7x7SilVKfw+IbML6shs7iKmcN6WR1FKaXcmscX/taMxuWQtfCVUurMPL7wt6SXEB4cwOg+On6vlFJn4vGFvzWjhOlDIrHpcshKKXVGHl342aXV5Bw/yaxhOh1TKaXa4tGFf3r+vY7fK6VU2zy78DNKiAoNJC4m1OooSinl9jy28I0xbD1cwoyhvRDR8XullGqLxxZ+ZnEV+eU1OpyjlFJO8tjC35rROH6vb9gqpZRzPLfwD5fQJyyIwb2CrY6ilFIewSML3xjD9sxSZgyN1PF7pZRykkcWflZJNUUVtUwbouP3SinlLI8s/J2ZpQBMGxJhcRKllPIcThW+iFwqIqkiki4iv27hcRGRFx2P7xWRya6P+n92ZJUSGRLIsGidf6+UUs5qs/BFxA94BbgMGAPcLCJjmm12GRDnuC0CFrs453/YkVlK/KAIHb9XSql2cOYIfxqQbozJMMacAj4EFjTbZgHwD9NoGxAuIn1dnBWAgvIajpZWM21IZGc8vVJKeS1/J7bpD2Q3+TgHmO7ENv2BvKYbicgiGl8BANSKyP52pW3ix3+GH3f0i10jCii2NsJZ0fzW8uT8npwdPD//yI5+oTOF39K4ienANhhjlgJLAUQkwRgT78Tf75Y0v7U0v3U8OTt4R/6Ofq0zQzo5wIAmH8cCuR3YRimllIWcKfydQJyIDBGRQOAmYGWzbVYCCx2zdWYAZcaYvOZPpJRSyjptDukYY+pF5EFgDeAHvGWMSRaR+xyPLwFWAZcD6UA1cJcTf/fSDqd2D5rfWprfOp6cHXw4vxjzvaF2pZRSXsgjz7RVSinVflr4SinlI7qs8EXkehFJFhG7iLQ6JaqtZRysIiKRIrJWRA45/mxxIR8RyRKRfSKy52ymT7mCuy2J0V5O5J8rImWOfb1HRJ6wImdrROQtESls7XwTd97/TmR3930/QES+E5EUR+881MI27rz/ncnf/p+BMaZLbsBoGk8YWAfEt7KNH3AYGAoEAknAmK7K2Eb+vwC/dtz/NfDnVrbLAqLcIG+b+5LGN9pX03gexQxgu9W525l/LvC51VnP8D2cB0wG9rfyuDvv/7ayu/u+7wtMdtzvAaR52O+/M/nb/TPosiN8Y0yKMSa1jc2cWcbBKguAdxz33wGutjCLM9xqSYwOcOffBacYYzYApWfYxG33vxPZ3ZoxJs8Ys8txvwJIofHs/6bcef87k7/d3G0Mv7UlGtxBb+M4t8DxZ0wr2xngKxFJdCwlYRVn9qU7729ns80UkSQRWS0iY7smmsu48/53hkfsexEZDEwCtjd7yCP2/xnyQzt/Bs4srdCeYF8DfVp46LfGmM+ceYoWPtdl80bPlL8dTzPbGJMrIjHAWhE56Dha6mouWxLDIs5k2wUMMsZUisjlwAoaV2z1FO68/9viEfteREKB5cDDxpjy5g+38CVutf/byN/un4FLC98Yc9FZPoWlSzScKb+IFIhIX2NMnuNlX2Erz5Hr+LNQRD6lcWjCisL39CUx2szW9B+AMWaViLwqIlHGGE9ZGMud9/8ZecK+F5EAGsvyfWPMJy1s4tb7v638HfkZuNuQjjPLOFhlJXCH4/4dwPdesYhIiIj0OH0fuATo8IqgZ8nTl8RoM7+I9BFpvCiCiEyj8fe5pMuTdpw77/8zcvd978j2JpBijHmulc3cdv87k79DP4MufNf5Ghr/R60FCoA1js/3A1Y12e5yGt+RPkzjUJDl75g7cvUCvgEOOf6MbJ6fxhklSY5bstX5W9qXwH3AfY77QuPFbQ4D+2hl9pQb53/QsZ+TgG3ALKszN8v/AY1LhNc5fvd/5Cn734ns7r7vz6VxeGYvsMdxu9yD9r8z+dv9M9ClFZRSyke425COUkqpTqKFr5RSPkILXymlfIQWvlJK+QgtfKWU8hFa+Eop5SO08JVSykf8f3qv7hKsciQHAAAAAElFTkSuQmCC\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": 4,
   "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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAADSCAYAAABuMkW8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZfbA8e8hgFQB6SJNRRRQUJBiAXSxYUERFUVRV+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+uCgADot2D4ZuCy4XwFoG9yvF3xhLwZKAVWBFsG2F4CPguNVBD4GHg+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+I7B1FpBxWvZ2XJcDbqnpNIWPPK6a8LMcKtDnB7/VyHftXVW2U54urzgMuFpESQFdgqIhUDZ7XOo+nrMYKj6aq+ntkb2H3Q+7Dc5yLBS8rEqisCJLS7CbPVQAiUhYYCbygqp/m2H0Z1gxLsJ9g5+F3YCdQV0RK5Eio6gG/UHj5nkPyP//RcATwbJRfs9h5M18CUNV1wIPAKyLSTUQqiEgJEWmBtbnnpyKwNigcWwOX5Ng2FDhLRE4QkdLYFdvePhPvAGeLyGkikiYiZYL5RQ6K4C2swpoE8ptr5H3gbhGpErzmP3Js+xbYICJ3ikjZ4PjNRORYABG5VESqB4VI9hXNTqwKuZOIXCgiJUWkqoi0CPZ7HXg+e04TEakjIqdF8F7ArugKnDfFuVjzsiIhy4rR7GpOBZse4CdVfSqP932miPxNREphnd63YZ3mp2D9qe4QkVIi0hE4G3gvwjhzyvcckv/5LxIRqYP1xfqmoH3jjSdTCSL4Yt2KdTpciX1JX8Pm5Pg6n6deBzwkIhuxDoXv53jNOcD12BXpcmx0Rp6TrwVt412wDo6rsKuX24ngM6SqW4BHga+Cvhtt89jtQay6+FdsdM/bOZ6/EysYWgTbVwNvYCOXwDrdzhGRTdiQ8O5BH5DFWD+N27BOjzOw+VzAztt84BsR2QB8TuT9HAYATYL3MjLC5zgXE15WJFxZ0R/oIbvaEbsD54lNdpl9O1FVf8b6o70UvK+zgbODfk3bsZF+ZwTbXgF6qupPEcb5lwjO4V7PfxRcAgwO+uUlFFH1FgvnnHMuLCLyLvC+qqbsxZnY3FI/AO1VdWXY8RSWJ1POOeecc0VQYLVr0N79rYj8ICJzROTBPPYREekjIvNFZKaIHFM84TrnXOS8/HLOxUIko/m2ASer6qag09uXIvJpjvk5wNppGwW3NtjMum2iHq1zzhWOl1/OuWIXSYdAVdVNwa+lglvutsEuwFvBvt8AlUWkdnRDdc65wvHyyzkXCxGN5guGRs7ARoaMVdUpuXapw+6TeC1l98nCnHMuFF5+OeeKW0STdgZDJVuISGVghIg0U9XZOXaRvJ6W+wER6YWtNUT58uVbHn744fsQsnOuOP3xB/z+OzRvDiWjPK3vtGnTVqtq9ei+av6iVX6Bl2HOpbL8yq9CFZWquk5EJmBzdeQsjJZiM6JmOwibrTX38/tjc2rQqlUrnTp1amEO75yLgTPPhP33hxkzov/aIlLYpSeipqjlV/AaXoY5l6LyK78iGc1XPbiiy57mvhO2inVOHwE9g1ExbYH1qrq8CDE750Kwcyd89RWceGLYkUSHl1/OuViIpGaqNjBYRNKw5Ot9Vf1ERHoDqGo/bDr8ztgssVuAK4spXudcMZozB9avhxNOCDuSqPHyyzlX7ApMplR1JnB0Ho/3y3FfsaUGnHMJbFKwLnyy1Ex5+eWciwVfm88595cvv4Q6daB+/bAjcc65xOHJlHMOAFWrmTrxRJC8xrc555zLkydTzjkAFi2yKRGSqL+Uc87FhCdTzjkg+fpLOedcrHgy5ZwDrL9UpUrQtGnYkTjnXGLxZMo5B8DEiXD88ZCWFnYkzjmXWDyZcs6xfDn89BN07Bh2JM45l3g8mXLOMWGC/TzppFDDcM65hOTJlHOO8eOtv9TRe0xv6ZxzriCeTDnnmDAB2rf3/lLOObcvPJlyLsX9/jvMm+f9pZxzbl95MuVcihs/3n56fynnnNs3nkw5l+LGj4cqVaB587Ajcc65xOTJlHMpbvx46NABSnhp4Jxz+8SLT+dS2KJF8Ouv3sTnnHNF4cmUcynM55dyzrmi82TKuRQ2fjxUrerr8TnnXFF4MuVcilKFceNsSgTvL+Wcc/uuwCJUROqKyHgRmSsic0Tkpjz26Sgi60VkRnC7v3jCdc5Fy88/w5IlcOqpYUdSfLz8cs7FQskI9skEblPV6SJSEZgmImNV9cdc+01S1bOiH6JzrjiMHWs/Tzkl3DiKmZdfzrliV2DNlKouV9Xpwf2NwFygTnEH5pwrXp99BoceCg0bhh1J8fHyyzkXC4XqKSEiDYCjgSl5bG4nIj+IyKci4t1ZnYtj27db5/NkbuLLzcsv51xxiaSZDwARqQAMA25W1Q25Nk8H6qvqJhHpDIwEGuXxGr2AXgD16tXb56Cdc0UzeTJs3pw6yVQ0yq/gdbwMc87tIaKaKREphRVE6ao6PPd2Vd2gqpuC+6OBUiJSLY/9+qtqK1VtVb169SKG7pzbV2PHQlpaaixuHK3yK9juZZhzbg+RjOYTYAAwV1Wf28s+tYL9EJHWweuuiWagzrno+ewzaNsWKlUKO5Li5eWXcy4aPvoo/+2RNPMdD1wGzBKRGcFj9wD1AFS1H9ANuFZEMoEMoLuq6j7G7JwrRmvWwNSp8O9/hx1JTHj55ZwrkoED4eqr89+nwGRKVb8EpIB9+gJ9CxOccy4c48bZhJ2p0F/Kyy/nXFE8/TTccYdNIZM9nUxefN5j51LMZ59Z816rVmFH4pxz8UnVkqg77oCLLoJPPsl//4hH8znnEp8qjBkDJ58MJf3b75xze8jMhF69rHnvuuugTx8bsJMfr5lyLoXMnm1LyHTuHHYkzjkXfzIyoFs3S6QeeAD69i04kQKvmXIupYwaZT89mXLOud2tXg1dutg8fC+9BDfcEPlzPZlyLoWMGgXHHAMHHhh2JM45Fz8WLoTTT4fFi+H99612qjC8mc+5FLFmDXz9NZx5ZtiROOdc/Pj2W5t3b80a+PzzwidS4MmUcyljzBjIyvJkyjnnsn38sa0EUaGCXWyecMK+vY4nU86liFGjoHp1OPbYsCNxzrnwvfoqnHsuNG1q/aQaN9731/JkyrkUkJkJ//2vdTwv4d9651wKy8qCu+6yaQ86d4YJE6BmzaK9pndAdy4FfPMNrF3rTXzOudS2ZQtcfjkMHQq9e9uovWjMuefJlHMpYNQoKzBSYQkZ55zLy++/29QH06fDM8/ArbeC5LvYVOQ8mXIuBYwaZR0rK1UKOxLnnIu9adPgnHNgwwb46CM466zovr73nnAuyS1cCLNmwdlnhx2Jc87F3tChcOKJVjv/1VfRT6TAkynnkt6IEfbzvPPCjcM552JJFR59FC64AFq0sPmkjjqqeI7lzXzOJbkRI6wgadgw7Eiccy42MjLgmmsgPR169IA33oAyZYrveF4z5VwS++MPm4jOa6Wcc6liyRJr1ktPh0cegbffLt5ECrxmyrmk9uGHVtXtyZRzLhV88YU1623bZh3NY9VX1GumnEtiI0bAIYdAs2ZhR+Kcc8VHFfr0gb/9DapWtf5RsRx048mUc0lq/XoYN85qpaI1l4pzzsWbjAy44gq46SYbqTdlStGWhtkXBSZTIlJXRMaLyFwRmSMiN+Wxj4hIHxGZLyIzReSY4gnXOReR9HRGHXojO3bAeW93tc4DKcjLL+eSUHo6NGgAJUqw5KB2nNh0DW+9BQ8+CMOHw/77xz6kSPpMZQK3qep0EakITBORsar6Y459zgAaBbc2wKvBT+dcrKWnQ69eDN8ymFosp+2KkdBrjG3r0SPc2GLPyy/nkklQvrFlC2PpRI/f09lGST6+bQJn3d8xtLAKTKZUdTmwPLi/UUTmAnWAnIVRF+AtVVXgGxGpLCK1g+e6AixbBjNmwJw5Nvpq1SpbmDYtDSpUgAMPhLp1bXh7kyZQunTYEbu4du+9bNxSglGcyVUMoARqC1Lde2/KJVNefjmXZO69l51btvII9/MgD9CEHxnG+TQeug2e+S20sAo1mk9EGgBHA1NybaoDLMnx+9Lgsd0KIxHpBfQCqFevXuEiTSKqNuJgxAgYPRrmz9+1rWxZqF7dEqasLJv6fvXqXdtLl4bWrW2l686dbQIy7w/jdrN4MR/Tna2UpTvv7fZ4Kitq+RW8hpdhzoVlzRpWLdpMDz5lLKfSk8G8wnWUZwssDvcfYcQd0EWkAjAMuFlVN+TenMdTdI8HVPuraitVbVW9evXCRZoENm+G55+Hww+Hk06C/v3hsMPssS++gLVrrQJh0SKYNw8WLLBaqq1b4eef4b334MYbbZ977rGaqqOOgmefhRUrwn53Lm4cdBDv0Z06LOU4vt71eAr/849G+QVehjkXilWr4K67+Kpud47meybSnte5mkFcYYkUhF6+RZRMiUgprCBKV9XheeyyFKib4/eDgGVFDy85bNsGL7wABx9sq1RXrw6DB8OaNbYA7c03Q/v2UKVK3s/fbz9Lui66CJ5+2hZsXL4c+vWDihXhn/+0ZsCrroKfforte3Px58+jOvBfTuci/mNNfADlytm6CinIyy/nEtSKFXD77Wj9Bjz7ZCYdtv6XMlXL802Zk7iaAbuuguKgfItkNJ8AA4C5qvrcXnb7COgZjIppC6z3/gbmyy+tBumWW+DII2026i+/hJ497e+/r2rVgv/7P3u9uXOtP96771qfqgsusJotl4JWr2bk5xXYQWkuqjXR2oDr17dq0BTrLwVefjmXkJYvt5qHhg1Z++xAulb9gn/yDF3OS2Pagiq0eOMGK9fiqXxT1XxvwAlYlfdMYEZw6wz0BnoH+wjwMrAAmAW0Kuh1W7Zsqcls61bVG29UBdUGDVRHjy7+Y65cqfqvf6mWL69asqTqDTfYYy6F3HabnsZ/teFB2zQrK+xg9gRM1QLKhmjeiqv80hQow5yLuaVL7R9nmTKqaWk64dRH9aBa27VUKdXnn9fQy7T8yq+YFWq5b8lcEC1cqNqqlZ3dG29U3bQptsdfvly1d2/VtDTVypVV33gj/A+hi4Hff9eV+x2kaZKpd98ddjB5i3UyVZy3ZC7DnIupxYtVr79edb/9VNPSdMflV+l916/REiVUGzVSnTo17ABNfuWXz4AeZV99Ba1aWTPbiBHw4otQvnxsY6hVC159FWbNsg7qV19tHd5//jm2cbgYe+QRhu04h52axkUXhR2Mc84VYNEiuPZaW/PqtdfgsstYNH4hHea9wcMvH0DPnjB9OrRsGXagBfNkKoqGDbN1gapVs07i554bbjxHHAHjx8Prr8MPP0Dz5pbcaZ7jlFxCW7gQXn+dt6rfRpMmlkQ751xc+vVXuOYaOPRQGDDARk/Nn88Hp75O87PrMXu29QEeONDmWkwEnkxFyeDB1vH7mGOsduqQQ8KOyJQoYTVTc+fCKafYyMHOnX0qhaTz4IPMK9GYySsO5vLLfe4x51wcWrAA/v53aNQI3nrLRlEtWMDGp17l6ofrc+GFVgkwYwZcfHHYwRaOJ1NR8M47cOWV0KkTfP651UzFm1q14KOP4OWXYcIEG1n46adhR+Wi4scf4Z13eOuYFyhRAi69NOyAnHMuh19+gcsvt9WHhwyB66+32vS+fZn0W12aN7daqHvugYkToWHDsAMuPE+mimjoUPuMnHQSjBxZtOkOipsIXHcdTJ0KtWvDmWfCww/bTOsugd1/P1nlKvDW0pM45RRbfsg550I3d65d3R1xBHzwgc06vXAhvPgiW6vW4fbboUMHa0GZONGmiipVKuyg940nU0Xw1Vf2OWnXzmp94jmRyqlpU5g82WK//37r27VuXdhRuX0yfToMG8YXXV9k8dI0Lr887ICccylvzhxrp2va1EZi3Xqr9ZN67jmoXZvvv7eBWs88Yy19M2bA8ceHHXTReDK1j+bNgy5dbAb7Dz+M/Yi9oipXzvp5vfSSNfcde6y1FrkE869/wQEHMDjzEvbfP/xBD865FDZzJlx4ofUj+eQTuPNO+O03W7qjZk0yM+Gxx6BNG1s+bfRoG3meKJ3M8+PJ1D7YsAHOPtuazUaPhqpVw45o34jADTfYiL+NG+G446zPl0sQX34Jn37KplvuY+iHpbnwQlso2znnYmrGDOja1YaM//e/1vnpt9/g8cdt/TSssur44+Hee23X2bPhjDPCDTuaPJkqJFUbHTd/vvWXOvTQsCMquhNOgClTbH2/M86wkaouzqlaqVSrFh9Uu5bNm/EmPudcbE2bZk00Rx8N48ZZv5HffoNHHvmrlmH7dnjoIdtl4ULrf/7ee3DAAeGGHm0lww4g0bz0kvWje/JJ6ziXLOrXt4qOCy+0ZHHBAvs+lPB0Oz6NHWs9Nvv25bVB+3HEEYnf58A5lyC+/dYypFGjoHJlePBB61xeufJuu02dalNIzZwJ3btDnz5/VVQlHf9XWQjffgu33WaJ+O23hx1N9FWqZM3cvXpZ7WyPHnZV4eJMdq1UgwbMaN2LKVPsb+ZzSznnitXkydZ80aaN3X/kEZvF/P77d0ukMjKsu1SbNrBqlfUrHjIkeRMp8JqpiG3ZApddZsPOBw1K3n9cpUpBv3426eidd1onwWHDkqODYNIYOdIu+QYO5LU3S1GmDPTsGXZQzrmk9eWXVhM1dqxNpPjEEzbPTsWKe+w6caJNbv7LL1Yr9cwze1RYJSWvmYrQ3Xfbh2PgwOT/YIjAHXfAm29ah/ROnWDNmrCjcgDs3An33QeNG7Pp3EtJT7em2WTrf+CciwNffAEnnwwnnmhrkj39tE1xcOedeyRSq1fb5OYdOliLxtix8MYbyf//MpsnUxEYN87aev/xD/tcpYorr7RaqRkzoH17+P33sCNyDBliw2IeeoghH5Rk40bo3TvsoJxzSUPV/ul16AAdO9rEm889Z0nUP/+5RzNFVpZdeDduDG+/bXnW7Nl2EZ5KPJkqwObNu5YSeuKJsKOJvXPPtXmoliyxDs7z5oUdUQrbsQMeeABatIBu3ejXz6Zzads27MCccwlP1aqT2reHv/3Nhqy/+KINwbvlljxnpZ4923Kuq66CJk3g++/t/2SizbsYDZ5MFSC7f90bbyTODOfRdtJJNhfV5s02jcLs2WFHlKIGDrSC7ZFH+ObbEkyfbrMHJ2v/PedcDKja3FDHHQennmpTG/Tta0O6b7wxz8nrNm+Gu+6y6Q5+/NGm0/niC2jWLPbhxwtPpvLx44/Wee6KKyxZT2UtW8KkSVCypNX8zpgRdkQpZutW6wDarh107szzz9voS59byjm3T1Rt+HabNjZCb9kym458/nxbiLhMmTyf8t57cPjhNj1Qz57w88/WepPq0+ik+NvfO9VdgxWeeirsaOLD4Yfb1Ue5ctZ3bOrUsCNKIa++ap3WHnuMxUuEYcNsOgQfZemcKxRVm6ugVStbymPVKnj9devD0bs37Ldfnk/7/ntr0rv4YpviYNIkq5GqVi3G8cepApMpEXlTRFaKSJ6NOyLSUUTWi8iM4HZ/9MOMvSFDLHF48snknhujsA491Ia+VqpkzeqTJ4cdUQrYuNEm/urUCTp25KWX7OF//CPcsBJFqpZhzu0mKwuGD7e2uXPPhfXrref4L7/YTM2lS+f5tFWrrDtBy5bWF71/f/juO+vy4XaJpGZqEHB6AftMUtUWwe2hoocVrowMmwqhZUvrWOd216CBJVQ1algT+8SJYUeU5F580Uq0Rx9l0ya7iOzWzZb/cREZRIqVYc79JSvLlu1o0QLOP98mTRw8GH76yYZslyqV59N27LCi57DDLOe66SarvLrmGkhLi/F7SAAFJlOqOhFYG4NY4kafPrB4sfWXSvV24L2pW9dq7g46yJrb//e/sCNKUn/+aR/ELl2gdWsGDrQLyltuCTuwxJGKZZhz7NxpHZyOPNImo9uxA955xzoD9+xpHWDzoGrzAjdrBjffDK1b23Iwzz+fOnNG7YtopQrtROQHEflURJrubScR6SUiU0Vk6qpVq6J06OhatQoeewzOOcc6Wru9O/BAmDABDj4YzjoLxowJO6Ik9PTTsGEDPPwwmZnwwgvWB71Nm7ADSzpJU4a5FJeZCenplg1dfLE9NmSIDcPu0WOvSRTAN9/YYKvzzrPap48+soF+RxwRo9gTWDSSqelAfVVtDrwEjNzbjqraX1VbqWqr6nHaEemhh2zY55NPhh1JYqhZ06ZNOPxwqzwZPTrsiJLIH39YPXv37nDkkfznPzYzwh13hB1Y0kmqMsylqMxMa75r0gQuvdSa7z74AGbNsjIkn7a5efPgggvsQm3ePHjtNauNOvtsn3olUkVOplR1g6puCu6PBkqJSEL27//tN1uX7pprLDlwkalWzZr5mja1K5qPPw47oiTx+OOwbRs8+CBZWVZjeuSRVmvqoieZyjCXgnbs2DUF+RVX2IyZw4fb/DXduuXbV2XlSptKqkkTm5z53/+2mRF69cq3AsvlocjJlIjUErHcVURaB6+ZkCu5PfqoJe//+lfYkSSeAw6wdfyaN7c+jiNGhB1Rglu82DL7K6+ERo0YOdK6Otxzj/fji7ZkKsNcCtm+3UajHHaYjZSqUsWmPJg+3a5q8yko1q61suTgg+GVV+zp8+fbAgs+3cq+KTD3FJEhQEegmogsBR4ASgGoaj+gG3CtiGQCGUB3VdVii7iY/PorDBoE114LdeqEHU1iqlLFViM4/XTr7zhkiF0YuX3w8MP28777ULWZ+Bs1sqp4VzipUoa5FLFtm62G8PjjdtF17LE2Y3nnzgW2yW3YYP0un33WZly56CKrjWrcODahJ7MCkylVvbiA7X2BvlGLKCTZtVJ33hl2JImtUiXriN65szXTp6fbF9YVwrx5Vlhefz3Uq8foUTZh3sCBPiR5X6RKGeaS3NatNkvmE0/A0qW2KOdrr8FppxWYRG3eDC+/bH2B1661iqsHH7RuAy46vMEA69Q7aJC1E3utVNHtv7+NADn+eLjkEkuoXCE88IDNQnzPPWRlwb33WnV8jx5hB+aci7mMDJuv55BD4IYbbKK/sWPh66+tGSCfRGrTJquFOuQQqyho08Ym3Bw+3BOpaPMuZlitVMmStnCji44KFWxk39lnw2WX2UATX0cuAjNn2twwd90FNWvy/nvwww82Pcxe5tZzziWjLVus5umpp2xkb4cOVhB07FhgTdS6dVYT9fzzsGaNrVYxbJhd4LrikfLJ1LJl8PbbNoLvwAPDjia5lC9v62h26WL9qDMzfUb5At13n1Xt3X47O3bYYIijjto1XYxzLslt2mRrcT7zjA23O/lku8Dq0KHAp65ebejGpzQAABquSURBVH2iXnrJ+keddZbVbLdtG4O4U1zKJ1N9+thEsbfeGnYkyalcOZv4rWtXW/4pM9PWeXJ5mDLFTtYjj0CVKgzoBwsWWELqI/icS3IbN1p10rPPWlZ0yilw//0RLYK3dKnVQvXrZ62C559vo/WOPjoGcTsgxZOpjRvtw9e1q7Upu+JRtqxNldCtmy1KnplpfatdLv/6l62qfdNNbNxoHUSPP9468zvnktT69TYa77nnrHf4GWdYDXW7dgU+9fvvLff6z39sGZhLLrF1ZX3G8thL6WRqwAD7HP/zn2FHkvzKlLE2+4susj6UmZm2cKYLjB9vE3U99xxUqMBjd1s3iREjfAZi55LSunXWNPL883b/rLMsiWrdOt+nqdoAn2eegXHjrH/qP/5h5Wn9+jGK3e0hZZOpHTvsM3ziib7OWazstx+8/771/7n5Zkuobrst7KjigKp1bDjoILj2WubPt5yqZ0/v6+Bc0lm71paJevFFu5rv0sWSqJYt831aRga8+66VDT/+aCPPn3rK+vv6AsThS9lkauhQm++sr88uE1OlS1tfyh49rEYwM9Pn9mL0aJg82UbulCnDrbfaeXriibADc85FzZo1dgXfp4/1Mena1ZKoFi3yfdqvv1p/9AEDLA9r3twGTV14oZUTLj6kbDLVp4/Nwn/mmWFHknpKlbIrrOzpKLJHraWk7ImkDjkErrySTz+1tQ2fegpq1w47OOdcka1aZR2bXn7ZZs+84AIr8PKZ6Ckry6aS6tsXRo2yASjnnmt9TSOYGcGFICWTqenT4ZtvbAipj5IKR8mSdnVVsqRdnGVm2lyVKVdIDB3610RSm7aV4tprbWkH70/mXIJbscI6Nr3yirXRde9uF05Nm+71KatXW7n4yiu2Vl6NGvaU//s/6wXg4ldKJlOvvmpD9n0SyXClpe1aIuXBBy2hevjhFEqoMjNt6HPTptC9O/feak3PkyZ59b1zCWv5cnj6aRsqvm2bDbG79144/PA8d8/KsrEnAwbAyJG2fnG7dlYmnn++9TV18S/lkqk//7TlTS691DvtxYO0NCtESpa0mejXr7caw5RYg+7tt+Hnn2H4cCZ/m8ZLL1k1vs9S7FwC+v13a5/v39/6Llx6qU32dNhhee6+aJFdTA4caBdRBxxgU8dcdZVN1OsSS8olU4MHW43rddeFHYnLVqKE9b2uVMm6FqxYYXlGUl+Rbdtml56tWpFx2rlc1Qrq1oXHHgs7MOdcoSxZYisIv/GGzQDds6clUXlMXpiRYfPyvvmm9YkC6NTJcrAuXWwKGZeYUiqZysqytuh27QocQOFirEQJ615Qu7aN8lu1yqq8K1UKO7Ji8sYbdmnavz+33yHMnQtjxkDFimEH5pyLyKJFNuR2wAD7/YorbMbMhg132y0z0+aDSk+3eeM2brQLp/vus2W2GjSIeeSuGKRUMjVuHMybZx2dXXy67TaoWdMKmQ4d4NNPk3BU25YttmRM+/Z8uOUUXn7ZljM69dSwA3POFejXX60KedAguwq8+mqb3yXHjJmq8N13lkD95z9W216pkg3k69HDyraU6MqQQlIqmRowAKpUsU59Ln5deqmtqnL++XDccZZQ7aXvZmLq2xf++IOlL3/I368Sjj7am/eci3sLFljHzrfesk6evXvDHXdYNROWQM2caSs9DBlio/FKl7aJzXv0sGWhvBkveaVMMvXnn1bFes01/oFOBKedBhMm2Dxgbdva1d1pp4UdVRSsXw9PPsnWU8+h6xOt2b7dCt6k7h/mXCL75RdLotLTbZK8G26wJOrAA8nKgimTYfhwuy1caJVVHTtai1/Xrj7QKVWkTDI1ZIj1+f3738OOxEWqVSv49ls45xy7qnv+eVuDKqGnTnjuOXTtWnqXfpPvvrMEv3HjsINyzu1h7lxLorKvdm66CW6/ncxqtfjiCxj+qH1/ly+3HKtTJ0ugzjnH5odyqaXAKStF5E0RWSkis/eyXUSkj4jMF5GZInJM9MMsujfftGn4jz467EhcYdSvD199BWefbWVZ79426jghrV4Nzz3HC0cNZPAnVXngAZvV2BWvZCnDXIzMmWMLiDZtatnSbbex4ttFDD7qWS66qRbVq1viNGiQdUNIT7cBM6NHW/cpT6RSUyQ1U4OAvsBbe9l+BtAouLUBXg1+xo2ZM2HaNFtX0iWeChWsCv3ee23wzI8/2vp+deqEHVkhPfkk720+m1tnXkHXrjZfp4uJQSR4GeZiYOZMGxgydCg7y1Vk6mV9GV39ckaPL8/Up22XWrWs6e6ss6zbQbly4Ybs4keByZSqThSRBvns0gV4S1UV+EZEKotIbVVdHqUYi2zgQKuGveSSsCNx+6pECXj8cZvM7pprbGqL9PQEGQGXng533snY34+gJ6Nof/gK0tNr+lJGMZIMZZiLovR0uzJbvBjq1YNevdDvpvLryBmMK3Mm45rMYOwfzVj9VglKlLA+m488Yl0NWrRI8G4GrthEo89UHWBJjt+XBo/tURCJSC+gF0C9evWicOiCbd8O77xjE6JVqxaTQ7pidPHF1lTbrRucfrqtF/rAA3E8zDg9HXr14n9b2nIuIzmCuXy46AzKDHvahvi4eBDXZZiLouD7yJYtLKM24xadwLh7azFOnmcR9WEr1FoDp59hydOpp0LVqmEH7RJBNJKpvPJ0zWtHVe0P9Ado1apVnvtE25gx1lXliiticTQXC4cfDlOm2NIrDz8M//uf9V9o1CjsyPJwzz2M3tKBrgznMH7hM06lcsZKuzL2ZCpexHUZ5opo50503nzmf7aQr+/6kq8ynmci7fkZm2/lANZwUpkp3PFMfU4+2QaEeO2TK6xoJFNLgbo5fj8IWBaF142Kd9+1K4uEaA5yEStf3hKoTp1shF/z5raiw/XXEz/NZ6tWMXhxR67hdY5kFp9xKlVZa9sWLw41NLebuC7DXCGsXQuzZpHx3WymTdjIVz9U4OtlDfg6qw2raQycQWX+5Di+5hpe52TG0ZwfKLEVuC4r7OhdAotGMvURcIOIvId12lwfL30NNm2CDz+0WqlSpcKOxhWHSy+Fk06yflQ33mgT5r38sg3ECdPOTz7lnosW8BSD+RufM5RuVGb9rh28iSiexG0Z5vYiM9Pmf5o5k23T5zDr641Mm1OGaesOZhotmcU17KA0AI0qr+SsZus4rmMGx73xd474Yxwlclc81qufx0Gci1yByZSIDAE6AtVEZCnwAFAKQFX7AaOBzsB8YAtwZXEFW1gffmgLS158cdiRuOJUpw6MGmXTX9x+u9VS/eMf1pcq5hPmbdnC79c9yuWDT+J/3EDv1tPpM+t8SmVs2LVPuXI2f42LiUQuwxzWT2PmTJg5k3Xf/sKcaVuZtbAc0zKbM42WzKbrX4lTlbIZHNN4C7cev53jTilFu+OE6tVrAMF8BYdfCb0m25JO2fz76KJBVUO5tWzZUotb586qdeuq7txZ7IdycWLVKtVevVRFVGvUUH3hBdUtW4rpYO+8o1q/vh2sfn3Neuhhfa/OLXoAq7Vcya36+ivbNStrz/30nXeKKaD4B0zVkMqcaN9iUYYlhUg//9u3q86apZqerhtv/pdOaXujvln5Zr2VZ/Q0PtU6LFFbtMVuB5Tfqqe03aB33Z6pH3ygunCh2vctWvE4l0t+5ZfY9thr1aqVTp06tdhef/VqWyD31lutL41LLdOm2aLJX3xhn4O777bZ78uXj9IBcowKAphNU26kD+M5mVaHrSf940ocdliUjpVERGSaqrYKO45oKO4yLCnk+p4AVhP09NNsqduY+ROWMn/aeub9ksW8P/Znvh7MPBqxjF2TyJUpuYMmDTNo2rwUzY4tS7Nm0KyZLYnnHcVdLOVXfiXtcjJDh1qzus8tlZpatrS1/SZMsOa+G2+E++6Dyy+3WdSPOKKAF7juOnjtNcjK1Sk1LQ127vzr5zSO4XHuZjhdqcKfvFL5Hq6Z8xglk/ab5eJe9jxKixbt+rxmj+9fu9b662U3a+Wcb+nRR/ccYZp7Tqa89slNFV2/gTU/rWLxjYNYsqUTi6m367alHouur89yDtztaTXKbaTRQRmcekQpDj1mJ02PSqNZM2jYsBRpad7p1cW3pK2Zat8e1qyB2bP96sXZkjSvvAIffGDL0Rx1lM1V1bmz9bHaLflJS9szicrhN+ozknN5h0uZRisqsY7reZlbeY6q8me+z011yVQzJdJKa9WaSo0a7HGrXh2qVLFb5cq7bpUqFdOcaFWqwLp1ke9furS1luVcm6lcOejff1eylKNWSYHNlGfdfrVY1/NGVtU6khWLt7FyeSYrVggr1pZixYayrMioyMrtVVhBDbZSdrdD7sfWv1Kquizh0L+359C21WnUcn8OPRT237/op8G54pRf+ZWUydSyZdYp+aGHrDbCuWwrVth0GcOGwddf2/+TihWhTRurrTq0323U3rGYCmyiNNvtHwiVWcjB/ExjJtOORTQAoCVTuZR3uJKBVCLoYF6/Pvz2W2jvL94lUzJVu3YrPfvsqaxcyW63jRvzf17FirsnV+XLWx6T361MGUv4S5a0ZGy3+xd2peTmdQhKFiXIogQ7Scvz/nZKk0HZvd/SKrKxWkPWbS3Dug3Cet2fdVRmPZXYuZeGjDQyqVFyLTXLbqBGhQxqVtlOzepZHFQX6n38MvXWz6Iei6nG6l0Tevn3xCWglGvmGzHCfp5/frhxuPhTsybccovdli+3PlWTJsE339i8VRt3PJvn84Qs6rOIlkzjnzzDKYylMb/svpOPCkopdepYRU5uGRm28O26dfDnn/Yz9y378fXr7f7vv1u3ouxbRkZhFvQeXqT3UZIdu9KpnRlUXJ9J5dJbOFAX04Q5VGYdlVhPZdb9db/6gCep0agSNRtX5oBqJSlRIseIuZzST4FeH/joOZf0kjKZGjbMZslu0iTsSFw8q10bune3G1gt1coSNVlFdTZRge2UpgKbqMhG6rKEMmzb80WymwQj7U/ikl7ZsvZxKOpUYjt2WFKVnVzt3Gn9QDMzd7+f2fb4v2qf0tgZ1EFl5Xm/VM7EKbiVZOeug+asMWpwjvW7yq1+ffh7w8jeRPb3obD9rpxLMEmXTK1ebbUNd98ddiQu0YhATVZSk5WRPSF3HxPnoqhUKbsV3Jfo68K/+F99pnIkUrlrjB59NO+ReIWtVerRw78jLunFy8IbUfPhh1ZR4E18bp+ULZv/9uzew/XreyLl4kMkM9NWrWo3EfvsvvkmDBxo97Mfy/157tHDHstvH+cckIQ1U8OGQcOG0KJF2JG4hLRli119Z2Ts/vi119pwQOfizZ9/5j2ar3x5m94jv+SnoMTIa5Wci0hSJVPr1sHnn8NNN/l0CK4IcjZrOJcI/vwz7AicS2lJ1cz3ySfWadOb+JxzzjkXK0mVTA0bZsOVW7cOOxLnnHPOpYqkSaa2bIExY+C886BE0rwr55xzzsW7pEk7/vc/6zPcpUvYkTjnnHMulSRNMvXJJ7ZMQ/v2YUfinHPOuVSSFMmUqiVTp51mc9E555xzzsVKUiRT339vixufdVbYkTjnnHMu1SRFMvXxxzavVOfOYUfinHPOuVQTUTIlIqeLyM8iMl9E7spje0cRWS8iM4Lb/dEPde8++QTatoXq1WN5VOdcIoj38ss5l/gKnAFdRNKAl4FTgKXAdyLykar+mGvXSaoa84a2Zctg6lR47LFYH9k5F+/ivfxyziWHSGqmWgPzVXWhqm4H3gPiZgKCUaPsp/eXcs7lIa7LL+dccogkmaoDLMnx+9LgsdzaicgPIvKpiDSNSnQR+OQTW8y8WbNYHdE5l0DiuvxyziWHSJKpvJYM1ly/Twfqq2pz4CVgZJ4vJNJLRKaKyNRVq1YVLtI8bN1qCxufdZYvbOycy1PUyi+IfhnmnEsOkSRTS4G6OX4/CFiWcwdV3aCqm4L7o4FSIlIt9wupan9VbaWqrapHobf4pEm2jMwZZxT5pZxzySlq5VewPaplmHMuOUSSTH0HNBKRhiJSGugOfJRzBxGpJWJ1QyLSOnjdNdEONrcxY2ySzo4di/tIzrkEFbfll3MueRQ4mk9VM0XkBmAMkAa8qapzRKR3sL0f0A24VkQygQygu6rmrkqPus8+gxNOgPLli/tIzrlEFM/ll3MueRSYTMFfVd+jcz3WL8f9vkDf6IaWv2XLYNYsePLJWB7VOZdo4rH8cs4ll4SdAf2zz+znaaeFG4dzzjnnUlvCJlNjxkDNmnDUUWFH4pxzzrlUlpDJ1M6dMHYsnHqqT4ngnHPOuXAlZDI1fTqsWeNNfM4555wLX0ImU9n9pU45Jdw4nHPOOecSMpkaMwaOOQZq1Ag7Euecc86luoRLpjZsgMmTvYnPOeecc/Eh4ZKpSZMgMxM6dQo7Euecc865BEymxo+H/faDdu3CjsQ555xzLgGTqXHj4LjjoGzZsCNxzjnnnEuwZGrNGpgxA046KexInHPOOedMQiVTX3wBqnDyyWFH4pxzzjlnEiqZGj8eypeHY48NOxLnnHPOOZNQydS4cXDiiVC6dNiROOecc86ZhEmm/vgDfvzRm/icc845F18SJpkaP95+ejLlnHPOuXiSUMlU5crQokXYkTjnnHPO7ZIwydS4cdChA6SlhR2Jc84559wuESVTInK6iPwsIvNF5K48touI9Am2zxSRY6IZ5JIlsGCBzy/lnCu8sMsv51zyKzCZEpE04GXgDKAJcLGINMm12xlAo+DWC3g1mkFOnGg/O3SI5qs655JdPJRfzrnkF0nNVGtgvqouVNXtwHtAl1z7dAHeUvMNUFlEakcryEmTYP/94cgjo/WKzrkUEXr55ZxLfpEkU3WAJTl+Xxo8Vth99tmkSXD88d5fyjlXaKGXX8655Fcygn0kj8d0H/ZBRHph1egA20RkdgTHB2yOKcnrKNFTDVhdrEcovHiLyePJX7zFA/EXU+MYHy9q5RcUrQxLQPH22Ym2ZH9/kPzvMdbvr/7eNkSSTC0F6ub4/SBg2T7sg6r2B/oDiMhUVW0VwfFjIt7igfiLyePJX7zFA/EXk4hMjfEho1Z+QXyXYdHm7y/xJft7jKf3F0kz33dAIxFpKCKlge7AR7n2+QjoGYyKaQusV9XlUY7VOecKy8sv51yxK7BmSlUzReQGYAyQBrypqnNEpHewvR8wGugMzAe2AFcWX8jOORcZL7+cc7EQSTMfqjoaK3ByPtYvx30Fri/ksfsXcv/iFm/xQPzF5PHkL97igfiLKebxFFP5BfF3bqPN31/iS/b3GDfvT6wccc4555xz+yJhlpNxzjnnnItHMUumROQCEZkjIlkistfe9wUt/RDFeA4QkbEiMi/4WWUv+/0mIrNEZEZxjESKx6UuIoipo4isD87JDBG5vxhjeVNEVu5tCHqsz08E8cTs3ATHqysi40VkbvD9uimPfWJ9jiKJKabnaV9F8PfuEZzTmSLytYg0j3WMRVHQ+8ux37EislNEusUqtmiI5P0Fn8UZwWf1i1jGFw0RfEYricjHIvJD8B4Tqk9gPJZxeVLVmNyAI7A5ZiYArfayTxqwADgYKA38ADQppnieAu4K7t8FPLmX/X4DqhVTDAW+X6xj7KfYXDhtgSnF/HeKJKaOwCcx+ty0B44BZu9le6zPT0HxxOzcBMerDRwT3K8I/BIHn6FIYorpeSrGv/dxQJXg/hnFfW5j/f6CfdKAcVi/s25hxxzlv19l4EegXvB7jbBjLob3eE/2/zegOrAWKB123IV4f3FXxuV1i1nNlKrOVdWfC9gtkqUfoqULMDi4Pxg4t5iOk594XOoiln+DAqnqROzLvzcxPT8RxBNTqrpcVacH9zcCc9lz9u5Yn6NIYkoIBf29VfVrVf0z+PUbbI6qhBHh5/kfwDBgZfFHFF0RvL9LgOGqujjYPxnfowIVRUSACsG+mbGILRrisYzLS7z1mYrlsg41NZhLJvhZYy/7KfCZiEwTm/04muJxqYtIj9cuqDb+VESaFmM8BYnHpUBCOTci0gA4GpiSa1No5yifmCB+PkPRchV2dZw0RKQOcB7Qr6B9E9RhQBURmRCU8T3DDqgY9MVahpYBs4CbVDUr3JD2TTyWcdkimhohUiLyOVArj033quqHkbxEHo/t83DD/OIpxMscr6rLRKQGMFZEfgquBKIhqktdREkkx5sO1FfVTSLSGRgJNCrGmPIT6/NTkFDOjYhUwGoPblbVDbk35/GUYj9HBcQUT5+hIhORk7Bk6oSwY4myF4A7VXWnFPN6XiEpCbQE/gaUBSaLyDeq+ku4YUXVacAM4GTgEOz/2KQ8vpNxLR7LuJyimkypaqcivkTEyzoUNR4RWSEitVV1eVAdmGf1rqouC36uFJERWDNYtJKpqC51EauYcn6QVXW0iLwiItVUNYw1oGJ9fvIVxrkRkVJYIZOuqsPz2CXm56igmOLsM1QkInIU8AZwhqquCTueKGsFvBckUtWAziKSqaojww0rapYCq1V1M7BZRCYCzbF+OcniSuAJtc5F80XkV+Bw4Ntww4pcPJZxucVbM18kSz9Ey0fA5cH9y4E9as5EpLyIVMy+D5wKRHNh03hc6qLAmESkVtD+joi0xj5HYf0TiaulQGJ9boJjDQDmqupze9ktpucokpji7DO0z0SkHjAcuCzJajMAUNWGqtpAVRsAQ4HrkiiRAiv3TxSRkiJSDmiD9clJJouxmjdEpCY2EGxhqBEVQjyWcXmJas1UfkTkPOAlbDTBKBGZoaqniciBwBuq2ln3svRDMYX0BPC+iFyFfdguCOL8Kx6gJjAiKPNLAu+q6n+jFcDe3q+EuNRFhDF1A64VkUwgA+geXPVEnYgMwUZ+VRORpcADQKkcscT0/EQQT8zOTeB44DJglojMCB67B6iXI6ZYL5cSSUyxPk/7JIK/9/1AVeCVoJzI1DhZeDUSEby/hFbQ+1PVuSLyX2AmkIWV/dG8YC52EfwNHwYGicgsrDnszgSrAY7HMm4PPgO6c84551wRxFszn3POOedcQvFkyjnnnHOuCDyZcs4555wrAk+mnHPOOeeKwJMp55xzzrki8GTKOeecc64IPJlyzjnnnCsCT6acc84554rg/wEs8ErAJnRTXwAAAABJRU5ErkJggg==\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": [
    "## Recap on BPR\n",
    "S.Rendle et al. BPR: Bayesian Personalized Ranking from Implicit Feedback. UAI2009\n",
    "\n",
    "The usual approach for item recommenders is to predict a personalized score $\\hat{x}_{ui}$ for an item that reflects the preference of the user for the item. Then the items are ranked by sorting them according to that score.\n",
    "\n",
    "Machine learning approaches are tipically fit by using observed items as a positive sample and missing ones for the negative class. A perfect model would thus be useless, as it would classify as negative (non-interesting) all the items that were non-observed at training time. The only reason why such methods work is regularization.\n",
    "\n",
    "BPR use a different approach. The training dataset is composed by triplets $(u,i,j)$ representing that user u is assumed to prefer i over j. For an implicit dataset this means that u observed i but not j:\n",
    "$$D_S := \\{(u,i,j) \\mid i \\in I_u^+ \\wedge j \\in I \\setminus I_u^+\\}$$\n",
    "\n",
    "### BPR-OPT\n",
    "A machine learning model can be represented by a parameter vector $\\Theta$ which is found at fitting time. BPR wants to find the parameter vector that is most probable given the desired, but latent, preference structure $>_u$:\n",
    "$$p(\\Theta \\mid >_u) \\propto p(>_u \\mid \\Theta)p(\\Theta) $$\n",
    "$$\\prod_{u\\in U} p(>_u \\mid \\Theta) = \\dots = \\prod_{(u,i,j) \\in D_S} p(i >_u j \\mid \\Theta) $$\n",
    "\n",
    "The probability that a user really prefers item $i$ to item $j$ is defined as:\n",
    "$$ p(i >_u j \\mid \\Theta) := \\sigma(\\hat{x}_{uij}(\\Theta)) $$\n",
    "Where $\\sigma$ represent the logistic sigmoid and $\\hat{x}_{uij}(\\Theta)$ is an arbitrary real-valued function of $\\Theta$ (the output of your arbitrary model).\n",
    "\n",
    "\n",
    "To complete the Bayesian setting, we define a prior density for the parameters:\n",
    "$$p(\\Theta) \\sim N(0, \\Sigma_\\Theta)$$\n",
    "And we can now formulate the maximum posterior estimator:\n",
    "$$BPR-OPT := \\log p(\\Theta \\mid >_u) $$\n",
    "$$ = \\log p(>_u \\mid \\Theta) p(\\Theta) $$\n",
    "$$ = \\log \\prod_{(u,i,j) \\in D_S} \\sigma(\\hat{x}_{uij})p(\\Theta) $$\n",
    "$$ = \\sum_{(u,i,j) \\in D_S} \\log \\sigma(\\hat{x}_{uij}) + \\log p(\\Theta) $$\n",
    "$$ = \\sum_{(u,i,j) \\in D_S} \\log \\sigma(\\hat{x}_{uij}) - \\lambda_\\Theta ||\\Theta||^2 $$\n",
    "\n",
    "Where $\\lambda_\\Theta$ are model specific regularization parameters.\n",
    "\n",
    "### BPR learning algorithm\n",
    "Once obtained the log-likelihood, we need to maximize it in order to find our obtimal $\\Theta$. As the crierion is differentiable, gradient descent algorithms are an obvious choiche for maximization.\n",
    "\n",
    "Gradient descent comes in many fashions, you can find an overview on Cesare Bernardis thesis https://www.politesi.polimi.it/bitstream/10589/133864/3/tesi.pdf on pages 18-19-20. A nice post about momentum is available here https://distill.pub/2017/momentum/\n",
    "\n",
    "The basic version of gradient descent consists in evaluating the gradient using all the available samples and then perform a single update. The problem with this is, in our case, that our training dataset is very skewed. Suppose an item i is very popular. Then we habe many terms of the form $\\hat{x}_{uij}$ in the loss because for many users u the item i is compared against all negative items j.\n",
    "\n",
    "The other popular approach is stochastic gradient descent, where for each training sample an update is performed. This is a better approach, but the order in which the samples are traversed is crucial. To solve this issue BPR uses a stochastic gradient descent algorithm that choses the triples randomly.\n",
    "\n",
    "The gradient of BPR-OPT with respect to the model parameters is: \n",
    "$$\\frac{\\partial BPR-OPT}{\\partial \\Theta} = \\sum_{(u,i,j) \\in D_S} \\frac{\\partial}{\\partial \\Theta} \\log \\sigma (\\hat{x}_{uij}) - \\lambda_\\Theta \\frac{\\partial}{\\partial\\Theta} || \\Theta ||^2$$\n",
    "$$ =  \\sum_{(u,i,j) \\in D_S} \\frac{-e^{-\\hat{x}_{uij}}}{1+e^{-\\hat{x}_{uij}}} \\frac{\\partial}{\\partial \\Theta}\\hat{x}_{uij} - \\lambda_\\Theta \\Theta $$\n",
    "\n",
    "### BPR-MF\n",
    "\n",
    "In order to practically apply this learning schema to an existing algorithm, we first split the real valued preference term: $\\hat{x}_{uij} := \\hat{x}_{ui} − \\hat{x}_{uj}$. And now we can apply any standard collaborative filtering model that predicts $\\hat{x}_{ui}$.\n",
    "\n",
    "The problem of predicting $\\hat{x}_{ui}$ can be seen as the task of estimating a matrix $X:U×I$. With matrix factorization teh target matrix $X$ is approximated by the matrix product of two low-rank matrices $W:|U|\\times k$ and $H:|I|\\times k$:\n",
    "$$X := WH^t$$\n",
    "The prediction formula can also be written as:\n",
    "$$\\hat{x}_{ui} = \\langle w_u,h_i \\rangle = \\sum_{f=1}^k w_{uf} \\cdot h_{if}$$\n",
    "Besides the dot product ⟨⋅,⋅⟩, in general any kernel can be used.\n",
    "\n",
    "We can now specify the derivatives:\n",
    "$$ \\frac{\\partial}{\\partial \\theta} \\hat{x}_{uij} = \\begin{cases}\n",
    "(h_{if} - h_{jf}) & \\text{ if } \\theta=w_{uf}, \\\\\n",
    "w_{uf} & \\text{ if } \\theta = h_{if}, \\\\\n",
    "-w_{uf} & \\text{ if } \\theta = h_{jf}, \\\\\n",
    "0 & \\text{ else }\n",
    "\\end{cases} $$\n",
    "\n",
    "Which basically means: user $u$ prefer $i$ over $j$, let's do the following:\n",
    "- Increase the relevance (according to $u$) of features belonging to $i$ but not to $j$ and vice-versa\n",
    "- Increase the relevance of features assigned to $i$\n",
    "- Decrease the relevance of features assigned to $j$\n",
    "\n",
    "We're now ready to look at some code!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's implement SLIM BPR \n",
    "\n",
    "Once we have understood how BPR works, we can take the general formulation and adapt it to other models.\n",
    "In order to implement SLIM BPR we simply need to replace $\\hat{x}_{uij}$ with the corresponding prediction computed by the model that we want to optimize.\n",
    "\n",
    "$ \\hat{x}_{uij} = \\hat{r}_{ui} - \\hat{r}_{uj}= \\sum_{k=0}^{I} r_{uk}S_{ki} - \\sum_{k=0, k}^{I} r_{uk}S_{kj} = \\sum_{k=0}^{I} r_{uk}(S_{ki} - S_{kj})$ \n",
    "\n",
    "#### Notice that we have two parameter sets for each sample, the weights of the positive item and those of the negative one, we have to update both and so we compute two gradients:\n",
    "\n",
    "$\\frac{\\partial BPR-OPT}{\\partial \\Theta} = \\sum_{(u,i,j) \\in D_S} \\frac{-e^{-\\hat{x}_{uij}}}{1+e^{-\\hat{x}_{uij}}} \\frac{\\partial}{\\partial \\Theta}\\hat{x}_{uij} - \\lambda_\\Theta \\Theta$\n",
    "\n",
    "$\\frac{\\partial}{\\partial S_{:i}}\\hat{x}_{uij} = r_{u:}$\n",
    "\n",
    "$\\frac{\\partial}{\\partial S_{:j}}\\hat{x}_{uij} = -r_{u:}$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What do we need for a SLIM BPR?\n",
    "\n",
    "* Item-Item similarity matrix\n",
    "* Computing prediction\n",
    "* Update rule\n",
    "* Training loop and some patience\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Movielens10M: Verifying data consistency...\n",
      "Movielens10M: Verifying data consistency... Passed!\n",
      "DataReader: current dataset is: <class 'Data_manager.Dataset.Dataset'>\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_tags, Value range: 1.00 / 69.00, Num features: 10217, feature occurrences: 108563, density 9.95E-04\n",
      "\tICM name: ICM_all, Value range: 1.00 / 69.00, Num features: 10237, feature occurrences: 130127, density 1.19E-03\n",
      "\n",
      "\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": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: We create a dense similarity matrix, initialized as zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_users, n_items = URM_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_item_S = np.zeros((n_items, n_items), dtype = np.float)\n",
    "item_item_S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: We sample a triplet\n",
    "\n",
    "#### Create a mask of positive interactions. How to build it depends on the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<69878x10681 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 4707743 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URM_mask = URM_train.copy()\n",
    "URM_mask.data[URM_mask.data <= 3] = 0\n",
    "\n",
    "URM_mask.eliminate_zeros()\n",
    "URM_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22782"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_id = np.random.choice(n_users)\n",
    "user_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get user seen items and choose one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  25,   27,   96,  101,  105,  108,  302,  311, 1003, 1300, 1487,\n",
       "       1490, 1551, 1601, 3161], dtype=int32)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_seen_items = URM_mask.indices[URM_mask.indptr[user_id]:URM_mask.indptr[user_id+1]]\n",
    "user_seen_items"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "101"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pos_item_id = np.random.choice(user_seen_items)\n",
    "pos_item_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To select a negative item it's faster to just try again then to build a mapping of the non-seen items"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4760"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "neg_item_selected = False\n",
    "\n",
    "# It's faster to just try again then to build a mapping of the non-seen items\n",
    "while (not neg_item_selected):\n",
    "    neg_item_id = np.random.randint(0, n_items)\n",
    "\n",
    "    if (neg_item_id not in user_seen_items):\n",
    "        neg_item_selected = True\n",
    "\n",
    "neg_item_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Computing prediction\n",
    "\n",
    "#### The prediction depends on the model: SLIM, Matrix Factorization... \n",
    "#### Note that here the data is implicit so we do not multiply for the user rating, because it is always 1, we just sum the similarities of the seen items."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_ui is 0.0000, x_uj is 0.0000\n"
     ]
    }
   ],
   "source": [
    "x_ui = item_item_S[pos_item_id, user_seen_items].sum()\n",
    "x_uj = item_item_S[neg_item_id, user_seen_items].sum()\n",
    "\n",
    "print(\"x_ui is {:.4f}, x_uj is {:.4f}\".format(x_ui, x_uj))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3 - Computing gradient\n",
    "\n",
    "#### The gradient depends on the objective function: RMSE, BPR... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_uij = x_ui - x_uj\n",
    "x_uij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The original BPR paper uses the logarithm of the sigmoid of x_ij, whose derivative is the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigmoid_item = 1 / (1 + np.exp(x_uij))\n",
    "sigmoid_item"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 - Update model\n",
    "\n",
    "#### How to update depends on the model itself, here we have just one paramether, the similarity matrix, so we perform just one update. In matrix factorization we have two.\n",
    "\n",
    "#### We need a learning rate, which influences how fast the model will change. Small ones lead to slower convergence but often higher results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-3\n",
    "\n",
    "item_item_S[pos_item_id, user_seen_items] += learning_rate * sigmoid_item\n",
    "item_item_S[pos_item_id, pos_item_id] = 0\n",
    "\n",
    "item_item_S[neg_item_id, user_seen_items] -= learning_rate * sigmoid_item\n",
    "item_item_S[neg_item_id, neg_item_id] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Usually there is no relevant change in the scores over a single iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_i is 0.0070, x_j is -0.0075\n"
     ]
    }
   ],
   "source": [
    "x_i = item_item_S[pos_item_id, user_seen_items].sum()\n",
    "x_j = item_item_S[neg_item_id, user_seen_items].sum()\n",
    "\n",
    "print(\"x_i is {:.4f}, x_j is {:.4f}\".format(x_i, x_j))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now we put everything in a training loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_triplet():\n",
    "    \n",
    "    non_empty_user = False\n",
    "    \n",
    "    while not non_empty_user:\n",
    "        user_id = np.random.choice(n_users)\n",
    "        user_seen_items = URM_mask.indices[URM_mask.indptr[user_id]:URM_mask.indptr[user_id+1]]\n",
    "        \n",
    "        if len(user_seen_items)>0:\n",
    "            non_empty_user = True\n",
    "\n",
    "    pos_item_id = np.random.choice(user_seen_items)\n",
    "\n",
    "    neg_item_selected = False\n",
    "\n",
    "    # It's faster to just try again then to build a mapping of the non-seen items\n",
    "    while (not neg_item_selected):\n",
    "        neg_item_id = np.random.randint(0, n_items)\n",
    "\n",
    "        if (neg_item_id not in user_seen_items):\n",
    "            neg_item_selected = True\n",
    "\n",
    "    return user_id, pos_item_id, neg_item_id    \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_one_epoch(item_item_S, learning_rate):\n",
    "\n",
    "    start_time = time.time()\n",
    "    for sample_num in range(n_users):\n",
    "\n",
    "        # Sample triplet\n",
    "        user_id, pos_item_id, neg_item_id = sample_triplet()\n",
    "        \n",
    "        user_seen_items = URM_mask.indices[URM_mask.indptr[user_id]:URM_mask.indptr[user_id+1]]\n",
    "\n",
    "        # Prediction\n",
    "        x_ui = item_item_S[pos_item_id, user_seen_items].sum()\n",
    "        x_uj = item_item_S[neg_item_id, user_seen_items].sum()\n",
    "        \n",
    "        # Gradient\n",
    "        x_uij = x_ui - x_uj\n",
    "\n",
    "        sigmoid_item = 1 / (1 + np.exp(x_uij))\n",
    "        \n",
    "        # Update\n",
    "        item_item_S[pos_item_id, user_seen_items] += learning_rate * sigmoid_item\n",
    "        item_item_S[pos_item_id, pos_item_id] = 0\n",
    "\n",
    "        item_item_S[neg_item_id, user_seen_items] -= learning_rate * sigmoid_item\n",
    "        item_item_S[neg_item_id, neg_item_id] = 0\n",
    "\n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 50000 == 0 or (sample_num +1) == n_users:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = (sample_num +1)/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, samples_per_second))\n",
    "         \n",
    "            \n",
    "    return item_item_S, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 50000 in 2.33 seconds. Samples per second 21429.37\n",
      "Iteration 69878 in 3.18 seconds. Samples per second 21974.59\n",
      "Iteration 50000 in 2.07 seconds. Samples per second 24114.87\n",
      "Iteration 69878 in 2.87 seconds. Samples per second 24363.31\n",
      "Iteration 50000 in 2.20 seconds. Samples per second 22715.52\n",
      "Iteration 69878 in 2.97 seconds. Samples per second 23512.38\n",
      "Iteration 50000 in 2.06 seconds. Samples per second 24321.61\n",
      "Iteration 69878 in 2.91 seconds. Samples per second 24004.01\n",
      "Iteration 50000 in 2.22 seconds. Samples per second 22537.08\n",
      "Iteration 69878 in 3.01 seconds. Samples per second 23247.39\n"
     ]
    }
   ],
   "source": [
    "learning_rate = 1e-6\n",
    "    \n",
    "item_item_S = np.zeros((n_items, n_items), dtype = np.float)\n",
    "\n",
    "for n_epoch in range(5):\n",
    "    item_item_S, samples_per_second = train_one_epoch(item_item_S, learning_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated time with the previous training speed is 3441.25 seconds, or 57.35 minutes\n"
     ]
    }
   ],
   "source": [
    "estimated_seconds = 8e6 * 10 / samples_per_second\n",
    "print(\"Estimated time with the previous training speed is {:.2f} seconds, or {:.2f} minutes\".format(estimated_seconds, estimated_seconds/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Common mistakes in using ML (based on last year's presentations)\n",
    "\n",
    "* Use default parameters and then give up when results are not good\n",
    "* Train for just 1 or 2 epochs\n",
    "* Use huge learning rate or regularization parameters: 1, 50, 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BPR for MF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What do we need for BPRMF?\n",
    "\n",
    "* User factor and Item factor matrices\n",
    "* Computing prediction\n",
    "* Update rule\n",
    "* Training loop and some patience\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: We create the dense latent factor matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_factors = 10\n",
    "\n",
    "user_factors = np.random.random((n_users, num_factors))\n",
    "item_factors = np.random.random((n_items, num_factors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.7659779 , 0.3143565 , 0.26425909, ..., 0.78056121, 0.86843444,\n",
       "        0.5185093 ],\n",
       "       [0.67924292, 0.49266765, 0.02975383, ..., 0.53037521, 0.82714536,\n",
       "        0.08431127],\n",
       "       [0.1162911 , 0.14537291, 0.47275693, ..., 0.03857484, 0.68371271,\n",
       "        0.38612273],\n",
       "       ...,\n",
       "       [0.65260609, 0.55514278, 0.69390257, ..., 0.41161263, 0.46967491,\n",
       "        0.02719706],\n",
       "       [0.50513679, 0.23330894, 0.66084978, ..., 0.50573155, 0.21205575,\n",
       "        0.81472249],\n",
       "       [0.92378248, 0.6938683 , 0.8175108 , ..., 0.20739599, 0.91706012,\n",
       "        0.64776784]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.52519126, 0.20423293, 0.43566583, ..., 0.0459627 , 0.80736068,\n",
       "        0.34986463],\n",
       "       [0.71346571, 0.80265357, 0.2766895 , ..., 0.80716581, 0.95323986,\n",
       "        0.93221669],\n",
       "       [0.54155147, 0.71543082, 0.09198144, ..., 0.61615629, 0.10739024,\n",
       "        0.7833634 ],\n",
       "       ...,\n",
       "       [0.78140293, 0.88435851, 0.29889247, ..., 0.67243557, 0.48467276,\n",
       "        0.1449844 ],\n",
       "       [0.54257934, 0.47111501, 0.44009355, ..., 0.36953606, 0.71166651,\n",
       "        0.93331242],\n",
       "       [0.85966759, 0.16708281, 0.38813289, ..., 0.37159929, 0.68552976,\n",
       "        0.11156988]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "item_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Computing prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(69855, 36, 3306)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "user_id, pos_item_id, neg_item_id = sample_triplet()\n",
    "(user_id, pos_item_id, neg_item_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_ui is 2.0311, x_uj is 2.0837\n"
     ]
    }
   ],
   "source": [
    "x_ui = np.dot(user_factors[user_id,:], item_factors[pos_item_id,:])\n",
    "x_uj = np.dot(user_factors[user_id,:], item_factors[neg_item_id,:])\n",
    "\n",
    "print(\"x_ui is {:.4f}, x_uj is {:.4f}\".format(x_ui, x_uj))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3 - Computing gradient\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.05262646490107015"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_uij = x_ui - x_uj\n",
    "x_uij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5131535805794873"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigmoid_item = 1 / (1 + np.exp(x_uij))\n",
    "sigmoid_item"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 - Update model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "regularization = 1e-4\n",
    "learning_rate = 1e-2\n",
    "\n",
    "H_i = item_factors[pos_item_id,:]\n",
    "H_j = item_factors[neg_item_id,:]\n",
    "W_u = user_factors[user_id,:]\n",
    "\n",
    "\n",
    "user_factors[user_id,:] += learning_rate * (sigmoid_item * ( H_i - H_j ) - regularization * W_u)\n",
    "item_factors[pos_item_id,:] += learning_rate * (sigmoid_item * ( W_u ) - regularization * H_i)\n",
    "item_factors[neg_item_id,:] += learning_rate * (sigmoid_item * (-W_u ) - regularization * H_j)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_i is 2.0481, x_j is 2.0707\n"
     ]
    }
   ],
   "source": [
    "x_ui = np.dot(user_factors[user_id,:], item_factors[pos_item_id,:])\n",
    "x_uj = np.dot(user_factors[user_id,:], item_factors[neg_item_id,:])\n",
    "\n",
    "print(\"x_i is {:.4f}, x_j is {:.4f}\".format(x_ui, x_uj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.022670430897814953"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_uij = x_ui - x_uj\n",
    "x_uij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_one_epoch(user_factors, item_factors, learning_rate):\n",
    "\n",
    "    start_time = time.time()\n",
    "    for sample_num in range(n_users):\n",
    "\n",
    "        # Sample triplet\n",
    "        user_id, pos_item_id, neg_item_id = sample_triplet()\n",
    "        \n",
    "        # Prediction\n",
    "        x_ui = np.dot(user_factors[user_id,:], item_factors[pos_item_id,:])\n",
    "        x_uj = np.dot(user_factors[user_id,:], item_factors[neg_item_id,:])\n",
    "        \n",
    "        # Gradient\n",
    "        x_uij = x_ui - x_uj\n",
    "\n",
    "        sigmoid_item = 1 / (1 + np.exp(x_uij))\n",
    "                \n",
    "        H_i = item_factors[pos_item_id,:]\n",
    "        H_j = item_factors[neg_item_id,:]\n",
    "        W_u = user_factors[user_id,:]\n",
    "\n",
    "\n",
    "        user_factors[user_id,:] += learning_rate * (sigmoid_item * ( H_i - H_j ) - regularization * W_u)\n",
    "        item_factors[pos_item_id,:] += learning_rate * (sigmoid_item * ( W_u ) - regularization * H_i)\n",
    "        item_factors[neg_item_id,:] += learning_rate * (sigmoid_item * (-W_u ) - regularization * H_j)\n",
    "\n",
    "        # Print some stats\n",
    "        if (sample_num +1)% 50000 == 0 or (sample_num +1) == n_users:\n",
    "            elapsed_time = time.time() - start_time\n",
    "            samples_per_second = (sample_num +1)/elapsed_time\n",
    "            print(\"Iteration {} in {:.2f} seconds. Samples per second {:.2f}\".format(sample_num+1, elapsed_time, samples_per_second))\n",
    "         \n",
    "        \n",
    "    return user_factors, item_factors, samples_per_second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 50000 in 1.77 seconds. Samples per second 28252.49\n",
      "Iteration 69878 in 2.66 seconds. Samples per second 26300.70\n",
      "Iteration 50000 in 1.86 seconds. Samples per second 26826.71\n",
      "Iteration 69878 in 2.60 seconds. Samples per second 26907.05\n",
      "Iteration 50000 in 1.85 seconds. Samples per second 26957.57\n",
      "Iteration 69878 in 2.56 seconds. Samples per second 27334.92\n",
      "Iteration 50000 in 1.69 seconds. Samples per second 29620.48\n",
      "Iteration 69878 in 2.35 seconds. Samples per second 29699.76\n",
      "Iteration 50000 in 1.95 seconds. Samples per second 25631.02\n",
      "Iteration 69878 in 2.68 seconds. Samples per second 26077.49\n"
     ]
    }
   ],
   "source": [
    "learning_rate = 1e-6\n",
    "num_factors = 10\n",
    "\n",
    "user_factors = np.random.random((n_users, num_factors))\n",
    "item_factors = np.random.random((n_items, num_factors))\n",
    "\n",
    "for n_epoch in range(5):\n",
    "    user_factors, item_factors, samples_per_second = train_one_epoch(user_factors, item_factors, learning_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final notes on BPR\n",
    "- BPR is very sensitive to the popularity bias of the dataset, so in unbalanced scenarios can lead to trivial recommendations (popular takes all)\n",
    "- More advanced strategies exist that fine-tune the sampling process (WARP)\n",
    "- The value range of your model predictions loses all meaning, you may have predictions >100 as well as negative numbers. This is not an issue if you only have to rank items but it becomes one if you want to buil hybrid methods and combine the prediction of multiple models. If you add a BPR model to the hybrid you likely will need to normalize the predictions in order to ensure they are within the same range as the other models, or you will need to adopt other strategies that build a hybrid based on the ranking and not on the absolute prediction value"
   ]
  },
  {
   "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
}