Practice extra 09-old - SLIM BPR and SLIM EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems 2020/21\n",
    "\n",
    "### Practice 5 - SLIM BPR\n",
    "\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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VGW+x/HPb1IISYAkJKEFCCWUANJCF0Us2LH3rou6uuquW91dd9fr1S1X3bUiltVV14JYUEFEpfeEHkJCCpCQQgqkkpBknvtHhnuzMWUSJjlTfu/XKy8mMyeT70zCN2eeec5zxBiDUkop72ezOoBSSqmuoYWvlFI+QgtfKaV8hBa+Ukr5CC18pZTyEVr4SinlI9osfBEJEpGtIrJLRJJF5E/NbNNNRD4UkXQR2SIisZ0RVimlVMc5s4dfA8w1xowHJgAXisj0JtvcDRwzxgwHngP+4tqYSimlTlebhW8aVDg+DXB8ND1aaz7wtuPyx8C5IiIuS6mUUuq0+TuzkYj4AUnAcOAlY8yWJpsMALIBjDF1IlIK9AaKmtzPAmABQEhIyORRo0adXnqllLKIMZCcV0pESCD9e3Xvsu+blJRUZIyJ6sjXOlX4xph6YIKIhAGfishYY8zeRps0tzf/gzUbjDGLgEUACQkJJjExsQORlVLKejsOH+PKlzfyys2TuGhcvy77viJyqKNf265ZOsaY48Bq4MImN+UAAx1h/IFeQElHQymllLtLPHgMgMmx4RYncZ4zs3SiHHv2iEh34Dxgf5PNlgK3Oy5fA3xvdFU2pZQX23awhNjewUT3CLI6itOcGdLpB7ztGMe3AR8ZY74UkSeARGPMUuAN4B0RSadhz/6GTkuslFIWM8aQeOgYc0dFWx2lXdosfGPMbmBiM9c/3uhyNXCta6MppZR7yiyqpKTyJFM8aDgH9EhbpZRqt8SDDW9RJsRGWJykfbTwlVKqnbYdPEZESCBDI0OsjtIuWvhKKdVOiQdLSBgcjqcdX6qFr5RS7XC0vJqDxVVM8bDhHNDCV0qpdklyzL9P8LA3bEELXyml2mXbwWMEBdgY07+X1VHaTQtfKaXaIfFQCRMGhhHo73n16XmJlVLKIpU1dSTnlnnk+D1o4SullNN2HD5Ovd0webDnjd+DFr5SSjltS1YxfjbxuAOuTtHCV0opJ23OLGbcgF6EdnNqZXm3o4WvlFJOOHGynp3Zx5k+tLfVUTpMC18ppZyw/fAxausN04d65nAOaOErpZRTNmd69vg9aOErpZRTPH38HrTwlVKqTd4wfg9a+Eop1SZvGL8HLXyllGqTN4zfgxa+Ukq1yRvG70ELXymlWuUt4/egha+UUq3ylvF70MJXSqlWecv4PWjhK6VUq7xl/B608JVSqkWVNXVeM34PWvhKKdWirQdLqK03nDk80uooLqGFr5RSLdhwoIhAf5tHnrC8OVr4SinVgvXpRUyJDScowM/qKC7RZuGLyEARWSUiKSKSLCIPN7PNHBEpFZGdjo/HOyeuUkp1jcLyGvbnlzPLS4ZzAJx527kOeNQYs11EegBJIrLSGLOvyXbrjDGXuj6iUkp1vY0ZRQDMHh5lcRLXaXMP3xiTZ4zZ7rhcDqQAAzo7mFJKWWn9gSLCggOI79/T6igu064xfBGJBSYCW5q5eYaI7BKR5SIyxgXZlFLKEsYYNqQXMXNYb/xsYnUcl3G68EUkFFgCPGKMKWty83ZgsDFmPPAC8FkL97FARBJFJLGwsLCjmZVSqlNlFlWSW1rtVeP34GThi0gADWX/njHmk6a3G2PKjDEVjsvLgAAR+cEzZYxZZIxJMMYkREV5z7iYUsq7bEhvGL/3lvn3pzgzS0eAN4AUY8yzLWzT17EdIjLVcb/FrgyqlFJdZf2BImLCuzMoItjqKC7lzCydWcCtwB4R2em47jFgEIAxZiFwDXC/iNQBJ4AbjDGmE/IqpVSnqqu3symzmEvG9cOxH+s12ix8Y8x6oNVHbYx5EXjRVaGUUsoqu4+UUl5d53Xj96BH2iql1H9Ym1aICFr4Sinl7VanFjI+JoyIkECro7icFr5SSjmUVJ5kV85xzh7hnbMItfCVUsph3YFCjIE5I7XwlVLKq61JLSQ8OIAzYsKsjtIptPCVUgqw2w1rDxQyOy7Kq5ZTaEwLXymlgOTcMooqTnrtcA5o4SulFABr0o4CMDtOC18ppbza6tRCxg3oRVSPblZH6TRa+Eopn1daVcv2w8e8djrmKVr4Simftz69CLsXT8c8RQtfKeXzVqcepWeQPxMGeud0zFO08JVSPs1uN6xKPcrZI6Px9/PuSvTuR6eUUm3YmXOcooqTnDc62uoonU4LXynl077dV4CfTZgzQgtfKaW82ncpR5kSG06v4ACro3Q6LXyllM/KLqkitaCc80b3sTpKl9DCV0r5rG9TCgC08JVSytt9m1LA8OhQYiNDrI7SJbTwlVI+qay6li2ZJZzrA7NzTtHCV0r5pDWphdTZjc8M54AWvlLKR32XUkB4cACTBoVbHaXLaOErpXxObb2dVamFnDMq2mtPdtIcLXyllM/ZkllC6YlaLojva3WULqWFr5TyOcv35tE9wM/rl0NuSgtfKeVT6u2GFckFnDMqiu6BflbH6VJa+Eopn5J06BhFFTVcOLaf1VG6nBa+UsqnfL03n0A/G3NH+c78+1PaLHwRGSgiq0QkRUSSReThZrYREXleRNJFZLeITOqcuEop1XHGGFYk5zM7LpLQbv5Wx+lyzuzh1wGPGmNGA9OBB0Qkvsk2FwFxjo8FwCsuTamUUi6wO6eUI8dPcOFY35qdc0qbf+KMMXlAnuNyuYikAAOAfY02mw/8yxhjgM0iEiYi/RxfqyxQUVPH3iOl7M8rI7WgnPzSagorajhWWYsxBgN087cRHhJIRHAgAyOCGRYdyojoUM6ICfO5N7OUb1i+Nx9/m3B+vO8cXdtYu17TiEgsMBHY0uSmAUB2o89zHNf9R+GLyAIaXgEwaNCg9iVVbcoqqmTZnjzWpBWy/dAx6uwGgPDgAGLCg4kK7caI6B7YbIIA1XV2jledJLe0mk2ZxVSdrAfA3yaMi+nFtCG9OT8+mokDw7H50MEpyjsZY/h6bx4zhvUmLDjQ6jiWcLrwRSQUWAI8Yowpa3pzM19ifnCFMYuARQAJCQk/uF21X3VtPZ/tOMLipBySDh0DYOyAnvzorKFMHRLBmH49ierRDZHWC9sYQ15pNfvzy9h28Bjbskp4Y30mC9dkENWjGxeO6cs1k2M4I6ZXm/ellDtKLSjnYHEVPzprqNVRLONU4YtIAA1l/54x5pNmNskBBjb6PAbIPf14qiXl1bW8t+Uwr6/LoqiihrjoUH5z0SiumDiAPj2D2n1/IkL/sO70D+vO3FENL3fLqmtZtf8oK5LzWZyUzTubDxHfryc3ThvE1ZMGEBzoe296Kc/15a48bILPHV3bWJv/Y6Vhd+4NIMUY82wLmy0FHhSRD4BpQKmO33eOervho8Rs/mdFKsWVJ5kdF8n9cyYwY2hvl+959wwKYP6EAcyfMICy6lo+33GEf2/N5vef7eW5lWncOTOW22bE+sSp4ZRnM8awdFcuM4dFEtWjm9VxLOPMLtos4FZgj4jsdFz3GDAIwBizEFgGXAykA1XAna6PqnZlH+exT/eQnFvG1NgI3rhjNBMGhnXJ9+4ZFMCtM2K5ZfpgEg8d45XVGTyzMo2FazK468whLDhrKD2CtPiVe9qVU8rhkioePGe41VEs5cwsnfU0P0bfeBsDPOCqUOo/1dbbeeH7dF5alU5UaDdeuHEil57Rz5KxdBFhSmwEU+6IICWvjBdXpfPC9+n8e8thHj4vjhunDiLAT4/nU+5l6c5cAv1szPPR6Zin6P9MN5ddUsVVL2/k+e8OMH9Cf1b89CwuG9/fLd44Hd2vJy/dNInPH5jF8OhQHv88mXl/X8vGjCKroyn1f+rthi935zJnZBS9uvv2q1AtfDe2/kARl724nkPFlbxy8ySevW6CW/7Cjh8YxgcLpvPG7QnU1Rtuem0LP/twJ0UVNVZHU4otWcUcLa/h8gn9rY5iOZ1m4aZeX5fJU8tSGB4dyqJbE9z+JMsiwrmj+zBreCQvrUpn4ZoMvtt/lN9fGs/Vkwa4xSsS5Zu+2JVLSKAf547yzYOtGtM9fDdjtxueWpbCk1+lcEF8Xz798Sy3L/vGggL8ePSCkSx/eDYj+/Tg54t3ce87Sbq3ryxxss7Osj35nB/fR48eRwvfrdTW2/nFx7tZtDaT22cM5uWbJxHioQs8DY/uwQcLpvO7S0azOq2Qec+t5ZvkfKtjKR+z7kAhpSdqdTjHQQvfTdTV23n4gx0s2Z7Dz84fwR8vH+PxyxnYbMI9s4fyxYNn0qdnEAveSeJ3n+2hurbe6mjKR3y+M5ew4ADOHO5bZ7ZqiRa+G6i3G36+eBfL9uTzu0tG89C5cV415j2ybw8+e2AWC84ayrubD3PNwo0cKq60OpbycmXVtXyzL59LxvUj0F+rDrTwLWe3Gx77ZA+f7czlF/NGcs9s71znI9DfxmMXj+b12xLILjnBpc+v5+u9ejC26jzLdudRXWvn2oSBbW/sI7TwLfbMylQ+TMzmJ3OH84APHAV4XnwfvvzJmQyNDuW+d7fz9PIU6u26jp5yvY+TchgeHcr4mF5WR3EbWvgW+mhbNi+tyuDGqQP52fkjrI7TZQZGBLP43hncPG0Qr67J5J63t1FWXWt1LOVFsooqSTx0jKsnxXjV8Ojp0sK3yPoDRTz26R5mx0XyxPyxPvdLGehv47+vHMeTV4xl3YEirnhpA5mFFVbHUl7ik+052ASunDjA6ihuRQvfAoeKK7n/vSSGR4fy8s2TfHrtmVumD+a9e6ZxvKqW+S9tYE1aodWRlIez2w1LknKYHRdF317tXyrcm/lu01jkxMl67n0nCT+b8NptCbrCJDBtaG+WPjiLmPBg7nprGx9sPWx1JOXBNmUWk1tazTWTY6yO4na08LuQMYbffrqH1IJy/n79BAZGBFsdyW3EhAez+L4ZzBoeya8/2cMz36TSsAirUu3zcVIOPYL8ffa8ta3Rwu9C72/N5pMdR3jk3BHMGRltdRy3E9rNnzduT+D6hIG88H06j360i5N1dqtjKQ9SVl3L13vzuWx8f4ICdCmFpjzzuH0PlH60gie+TGZ2XCQ/mev90y87KsDPxp+vHkdMeHeeWZlGflk1r9wy2S1XCVXu5/MdRzhRW8/1Ove+WbqH3wVO1tl55MMdBAf688y14z1+yYTOJiL85Nw4nr1uPFuzSrj+1U0cLa+2OpZyc8YY3ttymDH9e3KGzr1vlhZ+F3h2ZRp7j5Tx56vGEd2BE4z7qqsmxfDPO6dwuKSKaxduIrukyupIyo3tyD7O/vxybpo2yOemOTtLC7+TJR4s4dW1DQdXXTDGt0+v1hGz46J41zFt85qFG0krKLc6knJT/95ymJBAP+ZP0Ln3LdHC70TVtfX8aslu+vfqzu8uibc6jseaNCicD++djt3Ada9uYmf2casjKTdTeqKWL3fnMn/iAEI9dEnxrqCF34leWpVORmElT181zmPXtXcXo/r2ZMl9M+kR5M/Nr21mY7qeN1f9v0+351Bda+emqYOsjuLWtPA7SUpeGa+szuCqSQM4a4Suxe0Kg3oH8/F9M4kJD+aOf25j5b4CqyMpN2CM4d9bDzM+phdjB+ibta3Rwu8E9XbDr5fsplf3AH6vQzku1adnEB/eO53R/Xty/7tJLN+jSyz7us2ZJaQVVHDztMFWR3F7Wvid4N9bDrErp5THL4snPCTQ6jheJyw4kHfunsr4gWE8+P4OvtiVa3UkZaF/bsgiPDhAT2PoBC18FztWeZL/+SaNmcN6c/l4/QXsLD2DAnj7rqlMHhzOwx/s4JPtOVZHUhbILqliZUoBN00bpEfWOkEL38WeXZlGRU0df7hsjM4F7mSh3fx5684pTB/am0cX7+KjbdlWR1Jd7O2NB/ET4dbpsVZH8Qha+C60L7eM97Yc4tbpgxnZt4fVcXxCcKA/b94xhTOHR/LLJbt5b8shqyOpLlJZU8eHidlcNK6fLoPspDYLX0TeFJGjIrK3hdvniEipiOx0fDzu+pjuzxjDH79Iplf3AH56nu+cvcodBAX48dptCcwdFc1vP93LWxuyrI6kusCS7TmUV9dx56xYq6N4DGf28N8CLmxjm3XGmAmOjydOP5bnWZGcz9asEn4+byS9gnWhr64WFODHwlsmc0F8H/74xT5eX5dpdSTViex2w1sbDzJ+YBiTBoVbHcdjtFn4xpi1QEkXZPFYdfV2/roilbjoUG6Yogd+WCXQ38ZLN0/i4nF9efKrFBatzbA6kuok3+0/SmZhJXfp3n27uGoMf4aI7BKR5SIypqWNRGSBiCSKSGJhofecym5xUg6ZhZX8Yt5I/HQlTEsF+Nn4xw0TueSMfjy1bD+vrNbS9zbGGF5enU5MeHcuGdfP6jgexRXH+28HBhtjKkTkYuAzIK65DY0xi4BFAAkJCV5xOqMTJ+v5+7dpTB4crmfYcRMBfjb+cf0E/ET4y9f7sRvDA+foOQi8xdasEnYcPs5/zR+Dvw+fD7ojTrvwjTFljS4vE5GXRSTSGOMTi528tfEgBWU1vHDjJJ2G6Ub8/Ww8e914bAJ/W5GK3W74ybnN7ocoD/PKmgx6hwRyrZ7kpN1Ou/BFpC9QYIwxIjKVhmGi4tNO5gFKq2p5ZXU6c0dFM3VIhNVxVBP+fjaeuW4CNhGeWZlGvTE8ojOoPNq+3DJWpxby8wtG6IFWHdBm4YvI+8AcIFJEcoA/AAEAxpiFwDXA/SJSB5wAbjA+cvbpNzZkUVZdxy/mjbQ6imqBn034m+MsY3//9gB2Az89L05fjXmohWsyCAn00wOtOqjNwjfG3NjG7S8CL7oskYcoPVHLPzdkcdHYvozu19PqOKoVfjbhr1efgU3g+e8OYIzhZ+eP0NL3MAeLKvlydy53nzlEpz53kC7S3kFvbzxIeXUdD+oJyT2CzSb8+aozsInwwvfp1NsNv5g3Ukvfgzz//QEC/W386KyhVkfxWFr4HVBeXcsb67M4P74PY/rr+tuewmYTnrpyHDab8PLqDOwGfnWhlr4nyCis4LMdR7hn9lCie+gyCh2lhd8B/9p0iNITtTw0V2d9eBqbTXhy/lhs0jAebDeG31w0SkvfzT3/3QGCAvy4V/fuT4sWfjtV1tTx+rpMzhkZxbgY3bv3RDab8F/zx2ITYdHaTOrtht9dMlpL300dKChn6a5c7jt7GL1Du1kdx6Np4bfT+1sPc6yqVud0ezgR4U+Xj8Emwhvrs7Abw+OXxmvpu6G/f3eA4AA/FszWvfvTpYXfDrX1dv654SBTh0Togk1eQET4w2Xx2ER4c0MWxsAfLtPSdyd7j5Ty1e48HjxnuJ49zgW08Nth2Z48jhw/wZ8ub3G5IOVhRITfXzoaPxu8ti6Lertp2PPXNZEsZ4zhqWUphAcHsOBs3bt3BS18JxljWLQ2k2FRIcwdFW11HOVCIsJjF4/GZhNeXZOJ3ZiGMX4tfUutTi1kY0Yxf7wsnp5BOu/eFbTwnbQpo5jk3DL+fNU4LQIvJCL8+sJR2ER4ZXXD7J3/vkJ/1lapq7fz1LIUhkSGcNO0wVbH8Rpa+E56dW0mkaHduGLiAKujqE4iIvxy3kj8RHhxVTp2Ozytf+AtsTgphwNHK1h4yyQC/XVFTFfRwndCan45a9J0wSZfICI8esEIbDbh+e8OUG8Mf7n6DD3PQReqqKnj2ZVpJAwOZ96YvlbH8Spa+E54e9NBuvnbuFlfWvoEEeFn54/AJjgWXDP87ZrxWvpd5B/fplFUUcNrtyXojCkX08JvQ+mJWj7dfoT5E/rrtDAf88h5I7CJ8OzKNIyB/7lWS7+zpeaX8+aGg9wwZSATBoZZHcfraOG3YUlSDidq67ltRqzVUZQFHjo3rmGJ5RWp2I3hmWvH61mWOokxht9/tpeeQf78ct4oq+N4JS38Vtjthnc2H2LSoDDGDtBlFHzVA+cMRwT++nUqdgPPXael3xk+2X6ErQdL+PNV4/TVdCfRwm/F+vQisooqefj6CVZHURb78Zzh+Inw9PL92O2Gv98wgQAtfZcprarl6eUpTBwUxnV66sJOo4Xfin9tOkhkaCAXjdOZAgruPXsYfjbhya9SsBvD8zdO1NJ3kT99mcyxqlrevksPeOtM+tvaguySKr7bf5Qbpgyim79OxVQN7pk9lN9fGs/yvfk8+O/tnKyzWx3J432XUsAn24/wwJxhen6JTqaF34L3thzGJsJN0wZZHUW5mbvPHMIfLotnRXIBC95JpOpkndWRPFZpVS2PfbqHUX178KCeX6LTaeE3o7bezsdJ2cwdFU3/sO5Wx1Fu6M5ZQ3j6qnGsTSvkpte2UFJ50upIHumJL/dRVHGSv10zXo+o7QL6DDfju5SjFFWc5IYp+uaRatmNUwfxyi2TSckr45qFG8kuqbI6kkdZuiuXJdtzeGDOMD2ZUBfRwm/GR4nZRPfoxtkjoqyOotzcvDF9efeeaRSV13D1KxtJySuzOpJHOFxcxW8/2cOkQWE8pCcT6jJa+E3kl1azOvUo1ybE6Fxr5ZQpsRF8fP9M/GzCdQs3sSmj2OpIbq223s5DH+wAgX/cMFH/n3UhfaabWLI9B7tB5wKrdhnRpwdL7p9J315B3P7mVpYk5VgdyW399ev97Mw+zp+vOoOBEcFWx/EpWviN2O2GD7dlM31oBIN7h1gdR3mY/mHd+fi+mUwZEs6ji3fx168bDtJS/+/znUd4bV0Wt80YzCVn9LM6js/Rwm9kc1Yxh0uquF7frFUd1Cs4gLfunMqNUwfx8uoMfvzedp226bD3SCm/WrKbqUMi+P2l8VbH8Ula+I18tC2bHkH+XDRW9zxUxwX42XjqyrH87pLRrNiXz/Wvbia/tNrqWJYqLK/h3neSCA8O5OWbJ+kRyhZp81kXkTdF5KiI7G3hdhGR50UkXUR2i8gk18fsfGXVtSzfm8/8Cf31JCfqtIkI98weyuu3JZBZWMGlL6xnS6ZvvplbWVPHXW9to6TyJK/eOpnI0G5WR/JZzvyZfQu4sJXbLwLiHB8LgFdOP1bX+3pPPjV1dq6ZrMM5ynXOHd2HTx+YRc8gf256fQuvr8vEGN8Z16+tt3P/e9vZl1fGSzdP5IwYXePeSm0WvjFmLVDSyibzgX+ZBpuBMBHxuDGRT3ccYUhkCOP1ABDlYiP69ODzB2dx3uhonvwqhQf/vYOKGu8f17fbDb/6eDdr0wp56sqxzB3Vx+pIPs8VA2kDgOxGn+c4rvsBEVkgIokiklhYWOiCb+0aucdPsDmrmCsmDNBTqqlO0SMogIW3TOY3F41i+d485r+4nuTcUqtjdRq73fCrJbv5ZMcRHj1/BNdP0TWp3IErCr+5hmz2NasxZpExJsEYkxAV5T5HsS7dlYsxMH9Cf6ujKC8mItx79jDevWca5dV1XPnSRl5fl+l1UzdPlf3ipBwePjeOn+iRtG7DFYWfAzQe+I4Bcl1wv13msx1HmDgojNhInXuvOt/MYZF8/chZnD0yiie/SuG2N7dSUOYds3hq6up55MOd/1f2Pz1/hNWRVCOuKPylwG2O2TrTgVJjTJ4L7rdLpOSVsT+/nCsnNjsKpVSniAgJZNGtk3nqynEkHirhgufW8nFSjke/oVt6opbb39zK0l25/PLCkVr2bqjNM16JyPvAHCBSRHKAPwABAMaYhcAy4GIgHagC7uyssJ3hs51H8LcJl4zzuPeZlYcTx/kWpg2N4Jcf7+bni3fx+c4jPHXlOI9bcuBQcSU/+lciWUWV/P36CVyhO1BuSazao0hISDCJiYmWfO9T7HbDzD9/z5j+PXnjjimWZlG+zW43vLvlEH9Zvh+7gUfOi+OOWbEecba1Fcn5/HzxLgRYeMtkZg6PtDqSVxORJGNMQke+1qcPd9ucVUx+WbXujSjL2WzCbTNiWfmzs5k1vDdPL9/PBc+tZUVyvtsO81TX1vPkl/u4950khkSG8NVDs7Xs3ZxPF/6Xu/MIDvTjvNE6P1i5h/5h3Xn99im8fddUAv1s3PtOEje+tpnth49ZHe0/JB4s4eLn1/H6+ixumT6IxffN8LhhKF/U5hi+t6qrt/P13nzOHd2H7oHu/7JZ+ZazR0Qx6+HZvL8tm+dWpnHVyxuZHRfJT+bGMXVIhGW5jpZX89zKA3yw7TD9e3XnnbunMjvOfaZYq9b5bOFvySqhpPIkl4zra3UUpZrl72fj1umDuWriAN7bcohFazO57tVNJAwO57aZsVw4pm+XnQe2tKqWNzdk8dq6TE7W2bljZiw/v2AkId18tkI8ks/+tE4N58wZGW11FKVaFdLNnwVnDePW6bG8v/Uwb208yEPv7yAytBvXTI7hsvH9iO/Xs1OOEs8srODtjQdZnJRD1cl6LhnXj1/MG6nHrHgonyz8uno7K5IbhnN0ZUzlKboH+nHXmUO4Y2Ysaw4U8u6mQ7y2LpOFazIYGhXCvDF9mTUskoTY8A7/XhtjOFRcxbcpBXyxK5ddOaUE+AmXjx/A3WcOIb5/Txc/KtWVfLLwN2eeGs7RuffK89hswjkjozlnZDQllSdZvjePL3fl8draTF5ZnUGgv40zBvRidL+ejOrXg4HhwfTtFURkaDcC/W3424R6u6H0RC2lJ2o5VFxFZlEF+/PK2ZpVQr7jqN+xA3rym4tGceXEAUT3DLL4UStX8MnC/2pPLiGBfswZqW82Kc8WERLIzdMGc/O0wVTU1LE1q5gN6cXszjnOZzuOUL7Z+VU5+/YMIiE2nOlDe3Pm8EgdtvFCPlf4jWfn6HCO8iah3fyZO6rP/y1DbIzhyPET5B6vpqCsmqKKGmrr7dTWG/xsQq/uAfTqHkBMeHeGRIbQIyjA4kegOpvPFf6mzGKOVdXqCZSV1xMRYsKDiQnX+fGqgc8deLVsTx4hgX6cPUKHc5RSvsWnCr/ebli5r4BzRkXrcI5Syuf4VOHvOHyMooqTzBujB1sppXyPTxX+N/sKCPATnZ2jlPJJPlP4xhhyHp7TAAAJnElEQVRWJOczc1ikzkZQSvkknyn8A0crOFRcxQVjdGVMpZRv8pnC/yY5H4DzdSlkpZSP8pnCX7mvgImDwvQQcaWUz/KJws8rPcGunFIuiNfZOUop3+UThf/tvgIAzo/X4RyllO/yicL/Zl8BQ6NCGB4danUUpZSyjNcXfll1LZsyinXvXinl87y+8NelFVFnNzo7Rynl87y+8L/ff5Re3QOYMDDM6ihKKWUpry58u92wJu0oZ4+Iwt/Pqx+qUkq1yatbcM+RUooqTjJ3lJ6oXCmlvLrwv99/FBF07XullMLJwheRC0UkVUTSReTXzdx+h4gUishOx8c9ro/afqtSjzJxYBjhIYFWR1FKKcu1Wfgi4ge8BFwExAM3ikh8M5t+aIyZ4Ph43cU5262wvIbdOaU6nKOUUg7O7OFPBdKNMZnGmJPAB8D8zo11+lanHgXgHC18pZQCnCv8AUB2o89zHNc1dbWI7BaRj0VkoEvSnYZVqUfp07Mb8f16Wh1FKaXcgjOFL81cZ5p8/gUQa4w5A/gWeLvZOxJZICKJIpJYWFjYvqTtUFtvZ11aEeeMjEakufhKKeV7nCn8HKDxHnsMkNt4A2NMsTGmxvHpa8Dk5u7IGLPIGJNgjEmIiuq8mTOJB49RXlOnwzlKKdWIM4W/DYgTkSEiEgjcACxtvIGI9Gv06eVAiusitt+atEIC/IRZwyOtjKGUUm7Fv60NjDF1IvIgsALwA940xiSLyBNAojFmKfCQiFwO1AElwB2dmLlN6w4UMmlQOKHd2nx4SinlM5xqRGPMMmBZk+seb3T5N8BvXButY4oqakjOLeMX80ZaHUUppdyK1x1puyG9CIDZcTqco5RSjXld4a9NKyI8OIAx/XtZHUUppdyKVxW+MYZ1BwqZNTwSP5tOx1RKqca8qvDTCio4Wl7DWXG6WJpSSjXlVYW/7kDDwVxn6vi9Ukr9gFcV/toDRQyPDqV/WHeroyillNvxmsKvrq1nS2axzs5RSqkWeE3hJx48Rk2dXcfvlVKqBV5T+OvSG5ZTmDY0wuooSinllryn8NOKSBgcQXCgLqeglFLN8YrCL6qoYV9emc7OUUqpVnhF4W/KKAbQ1TGVUqoV3lH4mcX06ObP2P56diullGqJdxR+RjFTh0Tg7+cVD0cppTqFxzdkfmk1WUWVzBjW2+ooSinl1jy+8DdlNiyHrIWvlFKt8/jC35heTFhwAKP76vi9Ukq1xuMLf1NmMdOGRGDT5ZCVUqpVHl342SVV5Bw7wcxhOh1TKaXa4tGFf2r+vY7fK6VU2zy78DOLiQwNJC461OooSinl9jy28I0xbMooZvrQ3ojo+L1SSrXFYws/q6iS/LJqHc5RSikneWzhb8psGL/XN2yVUso5nlv4GcX07RlEbO9gq6MopZRH8MjCN8awJauE6UMjdPxeKaWc5JGFf7C4isLyGqYO0fF7pZRylkcW/rasEgCmDgm3OIlSSnkOpwpfRC4UkVQRSReRXzdzezcR+dBx+xYRiXV10Ma2HiwhIiSQYVE6/14ppZzVZuGLiB/wEnAREA/cKCLxTTa7GzhmjBkOPAf8xdVBG9uaVULC4HAdv1dKqXZwZg9/KpBujMk0xpwEPgDmN9lmPvC24/LHwLnSSW1cUFbN4ZIqpg6J6Iy7V0opr+XvxDYDgOxGn+cA01raxhhTJyKlQG+gqPFGIrIAWOD4tEZE9nYkNMCP/gI/6ugXu0YkTR6fh9H81vLk/J6cHTw//8iOfqEzhd/cnrrpwDYYYxYBiwBEJNEYk+DE93dLmt9amt86npwdvCN/R7/WmSGdHGBgo89jgNyWthERf6AXUNLRUEoppVzPmcLfBsSJyBARCQRuAJY22WYpcLvj8jXA98aYH+zhK6WUsk6bQzqOMfkHgRWAH/CmMSZZRJ4AEo0xS4E3gHdEJJ2GPfsbnPjei04jtzvQ/NbS/Nbx5Ozgw/lFd8SVUso3eOSRtkoppdpPC18ppXxElxW+iFwrIskiYheRFqdEtbWMg1VEJEJEVorIAce/zS7kIyL1IrLT8dH0ze0u527LYrSHE9nvEJHCRs/3PVbkbImIvCkiR1s63kQaPO94fLtFZFJXZ2yJE9nniEhpo+f+8a7O2BoRGSgiq0QkxdE7DzezjTs//87kb//PwBjTJR/AaBoOGFgNJLSwjR+QAQwFAoFdQHxXZWwj/1+BXzsu/xr4SwvbVVidtT3PJ/BjYKHj8g3Ah1bnbkf2O4AXrc7aymM4C5gE7G3h9ouB5TQcxzId2GJ15nZknwN8aXXOVvL3AyY5LvcA0pr5/XHn59+Z/O3+GXTZHr4xJsUYk9rGZs4s42CVxstHvA1cYWEWZ7nVshjt5M6/C04xxqyl9eNR5gP/Mg02A2Ei0q9r0rXOiexuzRiTZ4zZ7rhcDqTQsCJAY+78/DuTv93cbQy/uWUcTvtBukgfY0weNPwwgOgWtgsSkUQR2SwiVv9RcOb5/I9lMYBTy2JYzdnfhasdL8c/FpGBzdzuztz5990ZM0Rkl4gsF5ExVodpiWOYciKwpclNHvH8t5If2vkzcGZphfYE+xbo28xNvzXGfO7MXTRzXZfNG20tfzvuZpAxJldEhgLfi8geY0yGaxK2m8uWxbCAM7m+AN43xtSIyH00vFKZ2+nJXMddn3tnbAcGG2MqRORi4DMgzuJMPyAiocAS4BFjTFnTm5v5Erd6/tvI3+6fgUsL3xhz3mnehTPLOHSa1vKLSIGI9DPG5Dle9h1t4T5yHf9mishqGv4yW1X47VkWI8fNlsVoM7sxprjRp6/RyctydwJLf99PR+PyMcYsE5GXRSTSGOM2i5KJSAANZfmeMeaTZjZx6+e/rfwd+Rm425COM8s4WKXx8hG3Az94xSIi4SLSzXE5EpgF7OuyhD/kyctitJm9yXjr5TSMc3qSpcBtjtki04HSU8OG7k5E+p56r0dEptLQJcWtf1XXcWR7A0gxxjzbwmZu+/w7k79DP4MufNf5Shr+otYABcAKx/X9gWWNtruYhnekM2gYCrL8HXNHrt7Ad8ABx78RjusTgNcdl2cCe2iYUbIHuNsNcv/g+QSeAC53XA4CFgPpwFZgqNWZ25H9aSDZ8XyvAkZZnblJ/veBPKDW8bt/N3AfcJ/jdqHh5EIZjt+XZmevuWn2Bxs995uBmVZnbpL/TBqGZ3YDOx0fF3vQ8+9M/nb/DHRpBaWU8hHuNqSjlFKqk2jhK6WUj9DCV0opH6GFr5RSPkILXymlfIQWvlJK+QgtfKWU8hH/C0sr463bwk5KAAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAAlYAAADSCAYAAACIG474AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VGX2wPHvIfQioCAg0lREQQUlIjZAFxsWFFFRFHVVfrbFtvZV194bYkMRULNYABUFF3EBK6KASBGVohRBqvRQQs7vj3MjIaRMyMzcKefzPPNkMnNn5szNzJtz3/u+5xVVxTnnnHPOlV25sANwzjnnnEsVnlg555xzzkWJJ1bOOeecc1HiiZVzzjnnXJR4YuWcc845FyWeWDnnnHPORYknVm4HIvKbiHQOrt8hIq+GFEcnEVkUxms750rmbUX0iMgQETkz7Dh2RUn7X0TWi8g+ETxPJRH5SUT2jG6E8eeJVRIRkR4iMlFENojIsuD61SIisXg9VX1IVS8v6/OISFMRUREpH424wiYig0TkgbDjcK4o3lYkhkjaChE5BGgNfBCfqOJLVaur6rwIttsMvAbcGvuoYssTqyQhIjcBzwKPA/WBesCVwNFAxSIekxG3AJ1zCcHbiqTzf0CWerVugP8AF4tIpbADKRNV9UuCX4CawAbg7BK2GwS8CIwKtu8MnAp8D6wFFgL/LvCYi4D5wErgTuA3oHNw37+BN/Nt2x74GlgN/AB0ynffeOB+4CtgHfAJUCe4bwGgwPrgcmQhsVcJ4v8T+BG4GViU7/69gGHAcuBXoE+++9oBk4L3uBR4Kt99x+SLeSFwSXB7JeCJILalwEtAleC+TsAi4CZgGbAEuDS4rzewFdgSvJcPw/58+MUveRdvK5KvrQDmAcfk+/2HfO9/fbA/OgX3nQHMDGIcDxyY73EHBretDrY5o8Df+wXg4+A5v8KS7meC/fgTcGiE+7DY/V/I+1Ngv3xxPA+MDP72E4F9C2w/G+gY9nepTN/DsAPwSwR/JDgZyAHKl7DdIGANdmRaDqgcfPEPDn4/JGgYzgy2bxl8yToEjcdTwevs1FgCDbEGtUvwXCcEv9cN7h8PzAX2D75444FHgvuaBl+uIuMHHgG+AHYHGgEz8r6swetNBu7Gjrj3CRqjk4L7JwAXBderA+2D642DL+/5QAVgD6BNcN8zwIjg9WoAHwIPB/d1CvbDfcHjugAbgdr59vMDYX8u/OKXghdvK5KrrQCqBe+3bhH398aSnt2C/bUh2J8VgFuAOcH7rBBcvyP4/fjg/bTIF8cKoG3wtx6LJUy9gAzgAWBchPuwyP1fxHsomFitwhLc8kAW8FaB7UeQL5FLxoufCkwOdYAVqpqTd4OIfC0iq0UkW0Q65Nv2A1X9SlVzVXWTqo5X1enB79OAIUDHYNvuwEeq+rna+e27gNwiYrgQGKWqo4LnGoMd+XXJt81AVf1FVbOBd4A2pXiP5wIPquoqVV0I9M133+FYw3Ofqm5RO1//CtAjuH8rsJ+I1FHV9ar6TXB7T+BTVR2iqltVdaWqTg3GmVwB3BC83jrgoXzPl/ec9wWPG4X9U2lRivfjXBi8rUiutqJW8HNdwTtE5Bgs4TlDVdcC5wEjVXWMqm7FetGqAEdhPYTVsQR1i6qOBT7CEsU876nqZFXdBLwHbFLV11V1G/A2cGiE+7C4/R+J4ar6bfAZzWLnv/26fPslKaXEAME0sBKoIyLl8xpMVT0KIJiNkT9BXpj/gSJyBHaEcRB29FEJeDe4e6/826vqBhFZWUQMTYBzROT0fLdVAMbl+/2PfNc3Yl/0SO0QC3bKIf9r7yUiq/PdloEdNQFchh0x/iQivwL3qupH2NHU3EJeqy5QFZicbyyvBM+ZZ2X+f0678H6cC4O3FcnVVuTFWQPY9NcLiDTCEs6LVfWX4Oa9yPdeVTVXRBZiPYQ5wEJVzZ/szg/uy7M03/XsQn7Pi7mkfVjc/o9ESX/7GmzfL0nJE6vkMAHYDHTFznsXp+AAyP8A/YBTVHWTiDyDHdWCjQc4MG9DEamKdYEXZiHwhqpeUcrYC4upMEuwxm1m8HvjAq/9q6o2L/TJVWcD54tIOaAbMFRE9gge166Qh6zAGpJWqvp7ZG9hx5fchcc4Fw/eViRRWxEkqHmnRZcDiEgV4H3gGVX9ON/mi7FTtQTbCbYffge2AY1EpFy+5Kox8AulV+w+pPj9Hw0HAk9G+Tnjyk8FJgFVXQ3cC7wgIt1FpLqIlBORNtg5+uLUAFYFDWU74IJ89w0FThORY0SkInYkV9Rn4k3gdBE5SUQyRKRyUL9k7wjewnLstEFxtUzeAW4XkdrBc/4j333fAmtF5FYRqRK8/kEicjiAiFwoInWDBiXvSGcb1s3cWUTOFZHyIrKHiLQJtnsFeDqvZoqINBSRkyJ4L2BHeiXWZXEu3rytSMq2YhTbT7mClRz4SVUfK+R9nyoifxORCtiA+c3YgPuJ2PirW0Skgoh0Ak4H3oowzvyK3YcUv//LREQaYmO3vilp20TmiVWSCL5kN2IDFpdhX9iXsZofXxfz0KuB+0RkHTYY8Z18zzkTuAY7Ul2CzfIotNBbcC69KzY4cjl2VHMzEXyGVHUj8CDwVTDWo30hm92LdSn/is0SeiPf47dhjUSb4P4VwKvYDCiwAbszRWQ9Ns28RzBmZAE2ruMmbMDkVKxeDNh+mwN8IyJrgU+JfFzEAKBl8F7ej/AxzsWFtxVJ11b0B3rK9nONPYCzxApr5l2OVdWfsfFrzwXv63Tg9GAc1BZsxuApwX0vAL1U9acI4/xLBPuwyP0fBRcAg4NxfElLVP2shnPOORcWEfkP8I6qpu2Bmljtqh+ADqq6LOx4ysITK+ecc865KCmxazY4P/6tiPwgIjNF5N5CtqkkIm+LyByxpROaxiJY55wrLW/DnHPxFMkYq83A8araGjvnenIh570vA/5U1f2Ap4FHoxumc87tMm/DnHNxE8lgQlXV9cGvFYJLwfOHXYHBwfWhwN/yDcRzzrnQeBvmnIuniGYFBtMtp2IzTMao6sQCmzQkKBgWFEpbQ9E1TpxzLq68DXPOxUtEBUKD6ZdtRKQW8J6IHKSqM/JtUtiR3U6j4kWkN7b2EdWqVWt7wAEH7ELIzrlY+uMP+P13aN0ayke5hPDkyZNXqGrd6D5rybwNc86VVaTtV6maTVVdLSLjsVog+RulRVgl1kUiUh6rd7GqkMf3x2p2kJmZqZMmTSrNyzvn4uDUU2G33WDq1Og/t4iUdvmLqPI2zDm3qyJtvyKZFVg3OMrLK7XfGVttO78RwMXB9e7AWPU6Ds4lnW3b4Kuv4Nhjw44kerwNc87FUyQ9Vg2AwSKSgSVi76jqRyJyHzBJVUdg1WXfEJE52FFej6KfzjmXqGbOhDVr4Jhjwo4kqrwNc87FTYmJlapOAw4t5Pa7813fBJwT3dCcc/H2RbB+fSr1WHkb5pyLJ18r0Dn3ly+/hIYNoUmTsCNxzrnk5ImVcw4AVeuxOvZY8ApOzjm3azyxcs4BMH++lVlIsfFVzjkXV55YOeeA1Bxf5Zxz8eaJlXMOsPFVNWtCq1ZhR+Kcc8nLEyvnHACffw5HHw0ZGWFH4pxzycsTK+ccS5bATz9Bp05hR+Kcc8nNEyvnHOPH28/jjgs1DOecS3qeWDnnGDfOxlcdulMZTeecc6XhiZVzjvHjoUMHH1/lnHNl5YmVc2nu999h9mwfX+Wcc9HgiZVzaW7cOPvp46ucc67sPLFyLs2NGwe1a0Pr1mFH4pxzyc8TK+fS3Lhx0LEjlPPWwDnnysybUufS2Pz58OuvfhrQOeeixRMr59KY169yzrno8sTKuTQ2bhzssYevD+icc9HiiZVzaUoVxo61Mgs+vso556KjxOZURBqJyDgRmSUiM0XkukK26SQia0RkanC5OzbhOuei5eefYeFCOPHEsCOJHW+/nHPxVj6CbXKAm1R1iojUACaLyBhV/bHAdl+o6mnRD9E5FwtjxtjPE04IN44Y8/bLORdXJfZYqeoSVZ0SXF8HzAIaxjow51xsffIJ7LcfNGsWdiSx4+2Xcy7eSjWyQkSaAocCEwu5+0gR+UFEPhYRHwrrXALbssUGrqfyacCCvP1yzsVDJKcCARCR6sAw4HpVXVvg7ilAE1VdLyJdgPeB5oU8R2+gN0Djxo13OWjnXNlMmAAbNqRPYhWN9it4Hm/DnHPFiqjHSkQqYI1SlqoOL3i/qq5V1fXB9VFABRGpU8h2/VU1U1Uz69atW8bQnXO7aswYyMhIj4WXo9V+Bfd7G+acK1YkswIFGADMUtWnitimfrAdItIueN6V0QzUORc9n3wC7dtDzZphRxJb3n4556JhxIjIt43kVODRwEXAdBGZGtx2B9AYQFVfAroDV4lIDpAN9FBVLUXMzrk4WbkSJk2Cf/877Ejiwtsv51yZDBwIl18e+fYlJlaq+iUgJWzTD+gX+cs658IydqwVB02H8VXefjnnyuLxx+GWW6wsTV6JmpJ4vWXn0swnn9gpwMzMsCNxzrnEpGoJ1S23wHnnwUcfRf7YiGcFOueSnyqMHg3HHw/l/dvvnHM7ycmB3r3tFODVV0PfvjbZJ1LeY+VcGpkxw5ax6dIl7Eiccy7xZGdD9+6WVN1zD/TrV7qkCrzHyrm0MnKk/fTEyjnndrRiBXTtanX+nnsOrr12157HEyvn0sjIkXDYYbDXXmFH4pxziWPePDj5ZFiwAN55x3qtdpWfCnQuTaxcCV9/DaeeGnYkzjmXOL791ur6rVwJn35atqQKPLFyLm2MHg25uZ5YOedcng8/tBUoqle3A89jjin7c3pi5VyaGDkS6taFww8POxLnnAvfiy/CmWdCq1Y2rqpFi+g8rydWzqWBnBz4739t0Ho5/9Y759JYbi7cdpuVUujSBcaPh3r1ovf8PnjduTTwzTewapWfBnTOpbeNG+Hii2HoULjySpv9F+2afp5YOZcGRo60xiMdlrFxzrnC/P67lVOYMgWeeAJuvBGk2AWvdo0nVs6lgZEjbVBmzZphR+Kcc/E3eTKccQasXQsjRsBpp8XutXy0hXMpbt48mD4dTj897Eiccy7+hg6FY4+1XvuvvoptUgWeWDmX8t57z36edVa4cTjnXDypwoMPwjnnQJs2Vq/qkENi/7p+KtC5FPfee9aoNGsWdiTOORcf2dlwxRWQlQU9e8Krr0LlyvF5be+xci6F/fGHFb3z3irnXLpYuNBO/WVlwQMPwBtvxC+pAu+xci6lffCBdYd7YuWcSweffWan/jZvtkHqYYwt9R4r51LYe+/BvvvCQQeFHYlzzsWOKvTtC3/7G+yxh42nCmvCjidWzqWoNWtg7FjrrYpFrRbnnEsE2dlwySVw3XU242/ixOgtT7MrSkysRKSRiIwTkVkiMlNEritkGxGRviIyR0SmichhsQnXOReRrCxG7teHrVvhrDe62WCDNOTtl3MpKCsLmjaFcuVYuPeRHNtqJa+/DvfeC8OHw267hRteJGOscoCbVHWKiNQAJovIGFX9Md82pwDNg8sRwIvBT+dcvGVlQe/eDN84mPosof3S96H3aLuvZ89wY4s/b7+cSyVB+8bGjYyhMz1/z2Iz5fnwpvGcdnensKMDIkisVHUJsCS4vk5EZgENgfwNU1fgdVVV4BsRqSUiDYLHuhIsXgxTp8LMmTaLa/lyWzQ3IwOqV4e99oJGjWzKfMuWULFi2BG7hHbnnazbWI6RnMplDKAcagtk3Xln2iVW3n45l2LuvJNtGzfxAHdzL/fQkh8Zxtm0GLoZnvgt7OiAUs4KFJGmwKHAxAJ3NQQW5vt9UXDbDg2TiPQGegM0bty4dJGmEFWbufDeezBqFMyZs/2+KlWgbl1LnnJzrfz+ihXb769YEdq1sxW5u3SxYmc+fsbtYMECPqQHm6hCD97a4fZ0Vtb2K3gOb8OcC8vKlSyfv4GefMwYTqQXg3mBq6nGRliQOP8IIx68LiLVgWHA9aq6tuDdhTxEd7pBtb+qZqpqZt26dUsXaQrYsAGefhoOOACOOw7694f997fbPvsMVq2yjoX582H2bJg713qvNm2Cn3+Gt96CPn1smzvusB6sQw6BJ5+EpUvDfncuYey9N2/Rg4Ys4ii+3n57GicC0Wi/wNsw50KxfDncdhtfNerBoXzP53TgFS5nEJdYUgUJ1b5FlFiJSAWsUcpS1eGFbLIIaJTv972BxWUPLzVs3gzPPAP77GOradetC4MHw8qVtjju9ddDhw5Qu3bhj69UyRKw886Dxx+3xSSXLIGXXoIaNeCf/7RThZddBj/9FN/35hLPn4d05L+czHm8bacBAapWtbUd0pC3X84lqaVL4eab0SZNefLRHDpu+i+V96jGN5WP43IGbD8iSrD2LZJZgQIMAGap6lNFbDYC6BXMrmkPrPHxCebLL61n6YYb4OCDrQr2l19Cr172WdhV9evD//2fPd+sWTaW7z//sTFY55xjPV4uDa1YwfufVmcrFTmv/ud2nrhJE+seTbPxVeDtl3NJackS64Vo1oxVTw6k2x6f8U+eoOtZGUyeW5s2r15r7Vqitm+qWuwFOAbrFp8GTA0uXYArgSuDbQR4HpgLTAcyS3retm3bairbtEm1Tx9VUG3aVHXUqNi/5rJlqv/6l2q1aqrly6tee63d5tLITTfpSfxXm+29WXNzww5mZ8AkLaFtiOYlVu2XpkEb5lzcLVpk/zgrV1bNyNDxJz6oe9ffohUqqD79tIbepkXafsWtgSt4SeVGad481cxM27t9+qiuXx/f11+yRPXKK1UzMlRr1VJ99dXwP5AuDn7/XZdV2lszJEdvvz3sYAoX78QqlpdUbsOci6sFC1SvuUa1UiXVjAzdevFletc1K7VcOdXmzVUnTQo7QBNp++WV16Psq68gM9NOxb33Hjz7LFSrFt8Y6teHF1+E6dNtcPvll9tg+Z9/jm8cLs4eeIBhW89gm2Zw3nlhB+OccyWYPx+uusrW3Xr5ZbjoIuaPm0fH2a9y//O706sXTJkCbduGHWjpeGIVRcOG2TpFderYAPMzzww3ngMPhHHj4JVX4IcfoHVrS/S00PlOLqnNmwevvMLrdW+iZUtLqJ1zLiH9+itccQXstx8MGGAzr+bM4d0TX6H16Y2ZMcPGDA8caLUck40nVlEyeLANGj/sMOu12nffsCMy5cpZj9WsWXDCCTYDsUsXL8+Qcu69l9nlWjBh6T5cfLHXNnPOJaC5c+Hvf4fmzeH1120G1ty5rHvsRS6/vwnnnmsdAlOnwvnnhx3srvPEKgrefBMuvRQ6d4ZPP7Ueq0RTvz6MGAHPPw/jx9sMxY8/DjsqFxU//ghvvsnrhz1DuXJw4YVhB+Scc/n88gtcfLGtjDxkCFxzjfWy9+vHF781onVr65264w74/HNo1izsgMvGE6syGjrUPi/HHQfvv1+2EgqxJgJXXw2TJkGDBnDqqXD//Vbh3SWxu+8mt2p1Xl90HCecYEsgOedc6GbNsiO9Aw+Ed9+1Ctfz5sGzz7Jpj4bcfDN07GhnVj7/3EpRVagQdtBl54lVGXz1lX1mjjzSeoMSOanKr1UrmDDBYr/7bhsLtnp12FG5XTJlCgwbxmfdnmXBogwuvjjsgJxzaW/mTDuX16qVzeK68UYbV/XUU9CgAd9/b5O8nnjCzgZOnQpHHx120NHjidUumj0buna1KvoffBD/mX9lVbWqjQt77jk7JXj44XZGySWZf/0Ldt+dwTkXsNtu4U+YcM6lsWnT4NxzbazJRx/BrbfCb7/ZkiH16pGTAw89BEccYUu4jRplM9iTcYB6cTyx2gVr18Lpp9uptVGjYI89wo5o14jAtdfazMF16+Coo2yMmEsSX34JH3/M+hvuYugHFTn3XFvE2znn4mrqVOjWzaae//e/Nljqt9/g4YdtDTesE+voo+HOO23TGTPglFPCDTtWPLEqJVWbZTdnjo2v2m+/sCMqu2OOgYkTbb3BU06x2a8uwalaC1W/Pu/WuYoNG/DTgM65+Jo82U7dHHoojB1rY0t++w0eeOCvHoctW+C++2yTefNs7Ppbb8Huu4cbeiyVDzuAZPPcczYG79FHbdBdqmjSxDpAzj3XEse5c+27Uc5T78Q0ZoyN9uzXj5cHVeLAA1NrjIJzLoF9+61lSyNHQq1acO+9NjC9Vq0dNps0yUpUTZsGPXpA375/dWClNP+3WQrffgs33WQJ+s03hx1N9NWsaafFe/e2HtyePe1owyWYvN6qpk2Z2q43Eyfa38xrVznnYmrCBDutccQRdv2BB6x6+t1375BUZWfb8KojjoDly20c8pAh6ZFUgfdYRWzjRrjoIpvKPmhQ6v4Tq1ABXnrJCpzeeqsNMBw2LPUGFya199+3Q8GBA3n5tQpUrgy9eoUdlHMuZX35pfVQjRljhRofecRq99SosdOmn39uRdV/+cV6q554YqeOrJTnPVYRuv12+6AMHJj6HxIRuOUWeO01G8zeuTOsXBl2VA6AbdvgrrugRQvWn3khWVl2+jaVxys450Ly2Wdw/PFw7LG2Ltrjj1vZhFtv3SmpWrHCiqp37GhnOsaMgVdfTf3/l4XxxCoCY8faueF//MM+Y+ni0kutt2rqVOjQAX7/PeyIHEOG2PSa++5jyLvlWbcOrrwy7KCccylD1f7pdewInTpZkc+nnrKE6p//3On0RW6uHYS3aAFvvGE514wZdkCerjyxKsGGDduXNnrkkbCjib8zz7Q6VwsX2uDo2bPDjiiNbd0K99wDbdpA9+689JKVi2nfPuzAnHNJT9W6mTp0gL/9zaa+P/usTeW74YZCK2DPmGH512WXQcuW8P339n8y2eo6RpsnViXIG5v36qvJU1k92o47zmpdbdhgpRlmzAg7ojQ1cKA1cg88wDfflmPKFKtanKrj/ZxzcaBqtaeOOgpOPNHKJfTrZ1PD+/QptDjehg1w221WQuHHH61Ez2efwUEHxT/8ROSJVTF+/NEG3l1yiSXx6axtW/jiCyhf3nqHp04NO6I0s2mTDR498kjo0oWnn7ZZnF67yjm3S1RtGvgRR9hMv8WLrQz6nDm2SHLlyoU+5K234IADrORQr17w8892VsdL82znu6IIqtsnPTz2WNjRJIYDDrCjkqpVbazZpElhR5RGXnzRBrk99BALFgrDhlmJBZ+t6ZwrFVWrf5CZaUuILF8Or7xi4zyuvBIqVSr0Yd9/b6f9zj/fyiZ88YX1VNWpE+f4k0CJiZWIvCYiy0Sk0BNAItJJRNaIyNTgcnf0w4y/IUMsiXj00fSpvRGJ/faz6bQ1a9pp+AkTwo4oDaxbZ4XFOneGTp147jm7+R//CDesZJGubZhzO8jNheHD7fzdmWfCmjU26vyXX6wqdMWKhT5s+XIbctC2rY1j798fvvvOhoW4wkXSYzUIOLmEbb5Q1TbB5b6yhxWu7Gwrr9C2rQ3Kcztq2tSSqz33tFPyn38edkQp7tlnrXV78EHWr7eDy+7dbQkiF5FBpFkb5txfcnNtuZA2beDss60o4+DB8NNPNvW7QoVCH7Z1qzU9++9v+dd111mn1hVXQEZGnN9DkikxsVLVz4FVcYglYfTtCwsW2PgqP29cuEaNrEdv773t9Pz//hd2RCnqzz/tg9i1K7Rrx8CBdqB5ww1hB5Y80rENc45t22xA1MEHW7G7rVvhzTdt8HCvXjZgthCqVoP4oIPg+uuhXTtbkubpp9OzJtWuiFbacKSI/CAiH4tIq6I2EpHeIjJJRCYtX748Si8dXcuXw0MPwRln2CBtV7S99oLx42GffeC002D06LAjSkGPPw5r18L995OTA888Y+PXjzgi7MBSTsq0YS7N5eRAVpZlRuefb7cNGWLTuXv2LDKhAvjmG5uoddZZ1is1YoRNGDzwwDjFniKikVhNAZqoamvgOeD9ojZU1f6qmqmqmXUTdODSfffZVNJHHw07kuRQr56VYjjgAOtUGTUq7IhSyB9/WF98jx5w8MG8/bZVW7jllrADSzkp1Ya5NJWTY6f4WraECy+0U3zvvgvTp1sbUsz5u9mz4Zxz7KBt9mx4+WXrpTr9dC/nsivKnFip6lpVXR9cHwVUEJGknCfw22+2Tt4VV1ii4CJTp46dCmzVyo50Pvww7IhSxMMPw+bNcO+95OZaT+rBB1tvqoueVGrDXBraunV76fNLLrHqnMOHW02c7t2LHc+ybJmVqmrZ0gpB//vfVm2hd+9iO7ZcCcqcWIlIfRHLaUWkXfCcSbmy3IMPWlL/r3+FHUny2X13W1ewdWsbH/nee2FHlOQWLLAs/9JLoXlz3n/fhkbccYeP+4u2VGrDXBrZssVmsuy/v82yql3byihMmWJHuMU0FKtWWVuyzz7wwgv28DlzbGEHL+FSdiXmpCIyBOgE1BGRRcA9QAUAVX0J6A5cJSI5QDbQQ1U1ZhHHyK+/wqBBcNVV0LBh2NEkp9q1bUWEk0+2sZJDhtgBk9sF999vP++6C1VbAaB5c+uud6WTLm2YSxObN9sqDA8/bAdghx9uldK7dCnxvN3atTZO88knrYrLeedZL1WLFvEJPV2UmFip6vkl3N8P6Be1iEKS11t1661hR5Lcata0Qexduthp/aws+/K6Upg92xrOa66Bxo0ZNdKK8w0c6NOcd0W6tGEuxW3aZBU5H3kEFi2yRUJffhlOOqnEhGrDBnj+eRs7vGqVdWjde68NLXDR5ycVsAHBgwbZeWXvrSq73XazmSRHHw0XXGDJlSuFe+6x6sd33EFuLtx5p3XZ9+wZdmDOubjLzrYaQPvuC9dea4UEx4yBr7+20wPFJFXr11vv1L77WqfBEUdYcc/hwz2piiUfnob1VpUvb4tKuuioXt1mCJ5+Olx0kU1Y8XXtIjBtmtWeue02qFePd96CH36w8jNF1PFzzqWijRutR+qxx2yGcMeO1hB06lRiD9Xq1dZD9fTTsHKlrZIxbJgd7LrYS/vEavFieOMNmwm4115hR5NaqlWzNT67drUx2Dk5Xsm+RHfdZV1+N9/M1q02keKQQ7aXo3HOpbj1621t0CeesGl7xx9vB1sdO5b40BUrbAzVc88ZzX8OAAAbLklEQVTZeKrTTrMe7/bt4xC3+0vaJ1Z9+1qB2htvDDuS1FS1qhWZ69bNlqPKybF1p1whJk60nfXAA1C7NgNegrlzLTn1mYDOpbh166yb6cknLUM64QS4++6IFuVbtMh6p156yc4cnn22zfo79NA4xO12ktaJ1bp19kHs1s3OQbvYqFLFyi90726Lp+fk2LhsV8C//mUrfl93HevW2eDSo4+2iQDOuRS1Zo3N6nvqKRtZfsop1nN95JElPvT77y0Pe/ttW4rmggtsnVuvlB6utE6sBgywz/Q//xl2JKmvcmU7x3/eeTb+MifHFvV0gXHjrBDYU09B9eo8dLsNq3jvPa987FxKWr3aTpk8/bRdP+00S6jatSv2Yao2OeiJJ2DsWBvP+o9/WHvapEmcYnfFStvEautW+zwfe6yvuxYvlSrBO+/YeKHrr7fk6qabwo4qAajaQIi994arrmLOHMuvevXysRHOpZxVq2ypqmeftSP7rl0toWrbttiHZWfDf/5jbcOPP9oM9sces/HBvjhyYknbxGroUKut1s+r18RVxYo2DrNnT+spzMnx2mGMGgUTJtgMoMqVufFG20+PPBJ2YM65qFm50o7m+/a1cSjdullC1aZNsQ/79Vcbyz5ggOVkrVvbhKtzz7V2wiWetE2s+va1lQBOPTXsSNJPhQp25JVX4iJv9ltayitUte++cOmlfPyxrbX42GPQoEHYwTnnymz5chsI9fzzVqnznHOswSumkFRurpWq6tcPRo60yStnnmljUyOotuBClpaJ1ZQp8M03Ni3VZ1uFo3x5O+oqX94O2nJyrC5m2jUYQ4f+Vahq/eYKXHWVLS/h48+cS3JLl9pAqBdesPN4PXrYQVSrVkU+ZMUKaxdfeMHW7ttzT3vI//2fjRRwySEtE6sXX7QyAF6wMlwZGduXabn3Xkuu7r8/jZKrnBybTt2qFfTowZ032unpL77wLn7nktaSJfD44zblfPNmm6p3551wwAGFbp6ba/NWBgyA99+3tZWPPNLaxLPPtrGpLrmkXWL155+2xMqFF/qAv0SQkWENSvnyVgF/zRrrSUyLNfHeeAN+/hmGD2fCtxk895x19Xt1ZOeS0O+/2zn8/v1tfMOFF1oxqf33L3Tz+fPtwHLgQDug2n13K0dz2WVWFNglr7RLrAYPtl7Zq68OOxKXp1w5G7dds6YNRVi61HKOlD5S27zZDkkzM8k+6Uwuy4RGjeChh8IOzDlXKgsX2urGr75q1aZ79bKEqpDiiNnZVgP4tddsDBVA586Wj3XtamVpXPJLq8QqN9fOXR95ZIkTMVyclStnwxEaNLDZgsuXW7d4zZphRxYjr75qh6z9+3PzLcKsWTB6NNSoEXZgzrmIzJ9vU3cHDLDfL7nEqnM2a7bDZjk5Vm8qK8vq0q1bZwdRd91lS301bRr3yF2MpVViNXYszJ5tg6RdYrrpJqhXzxqcjh3h449TcHbcxo22bE2HDnyw8QSef96WVDrxxLADc86V6NdfrWt50CA7Irz8cqsZk686pyp8950lU2+/bb3wNWvahMCePa1tS4vhDmkqrRKrAQOgdm0bEOgS14UX2souZ58NRx1lyVUR4z6TU79+8McfLHr+A/5+mXDooX4K0LmEN3euDQR9/XUbFHrllXDLLdb9hCVT06bZChNDhtisvooVraB6z562NJWf6ksPaZNY/fmndcNecYV/uJPBSSfB+PFWZ6x9ezvqO+mksKOKgjVr4NFH2XTiGXR7pB1btlgjnNLjyZxLZr/8YglVVpYV4bv2Wkuo9tqL3FyYOAGGD7fLvHnWidWpk50V7NbNJ0mlo7RJrIYMsfHCf/972JG4SGVmwrffwhln2NHe00/bmlhJXY7hqafQVau4suJrfPedJfstWoQdlHNuJ7NmWUKVd+Rz3XVw883k1KnPZ5/B8Aft+7tkieVbnTtbMnXGGVZ/yqWvEstjishrIrJMRGYUcb+ISF8RmSMi00TksOiHWXavvWZLARx6aNiRuNJo0gS++gpOP93atSuvtJnMSWnFCnjqKZ45ZCCDP9qDe+6xasoutlKlDXNxMnOmLWjaqpVlTjfdxNJv5zP4kCc577r61K1rSdSgQTZUISvLJtuMGmXDrTypcpH0WA0C+gGvF3H/KUDz4HIE8GLwM2FMmwaTJ9ualy75VK9u3ex33mmTcH780dYbbNgw7MhK6dFHeWvD6dw47RK6dbPaoC4uBpHkbZiLg2nTbFLJ0KFsq1qDSRf1Y1Tdixk1rhqTHrdN6te303unnWZDE6pWDTdkl5hKTKxU9XMRaVrMJl2B11VVgW9EpJaINFDVJVGKscwGDrSu2gsuCDsSt6vKlYOHH7bCeVdcYeUysrKSZCZdVhbceitjfj+QXoykwwFLycqq58spxUkqtGEuirKy7ChtwQJo3Bh690a/m8Sv709lbOVTGdtyKmP+OIgVr5ejXDkb4/nAAzYcoU2bJB+K4OIiGmOsGgIL8/2+KLhtp0ZJRHoDvQEaN24chZcu2ZYt8OabVnytTp24vKSLofPPt9O53bvDySfbWqb33JPAU5ezsqB3b/63sT1n8j4HMosP5p9C5WGP21QhlwgSug1zURR8H9m4kcU0YOz8Yxh7Z33GytPMpwlsgvor4eRTLJE68UTYY4+wg3bJJhqJVWH5uxa2oar2B/oDZGZmFrpNtI0ebUNbLrkkHq/m4uGAA2DiRFv+5f774X//s/EOzZuHHVkh7riDURs70o3h7M8vfMKJ1MpeZkfMnlglioRuw1wZbduGzp7DnE/m8fVtX/JV9tN8Tgd+xmq47M5Kjqs8kVueaMLxx9tkEu+VcmURjcRqEdAo3+97A4uj8LxR8Z//2BFHUpwychGrVs2Sqc6dbaZg69a2qsQ115A4p9iWL2fwgk5cwSsczHQ+4UT2YJXdt2BBqKG5HSR0G+ZKYdUqmD6d7O9mMHn8Or76oTpfL27K17lHsIIWwCnU4k+O4muu4BWOZyyt+YFym4Crc8OO3qWIaCRWI4BrReQtbMDnmkQZm7B+PXzwgfVWVagQdjQuFi68EI47zsZd9eljxfmef94m9IRp20cfc8d5c3mMwfyNTxlKd2qxZvsGfhopkSRsG+aKkJNj9aWmTWPzlJlM/3odk2dWZvLqfZhMW6ZzBVupCEDzWss47aDVHNUpm6Ne/TsH/jGWcgU7JBs3KeRFnNs1JSZWIjIE6ATUEZFFwD1ABQBVfQkYBXQB5gAbgUtjFWxpffCBLXp5/vlhR+JiqWFDGDnSSmrcfLP1Xv3jHzb2Ku7F+TZu5PerH+TiwcfxP67lynZT6Dv9bCpkr92+TdWqVh/HxUUyt2EOG8sxbRpMm8bqb39h5uRNTJ9Xlck5rZlMW2bQ7a8kqnaVbA5rsZEbj97CUSdU4MijhLp19wSCGggHXAq9J9iyUnn8++iiTVVDubRt21ZjrUsX1UaNVLdti/lLuQSxfLlq796qIqp77qn6zDOqGzfG6MXefFO1SRN7sSZNNPe++/Wthjfo7qzQquU36SsvbNHc3J230zffjFFAiQ+YpCG1OdG+xKMNSwmRfv63bFGdPl01K0vXXf8vndi+j75W63q9kSf0JD7WhixUWzjGLrtX26QntF+rt92co+++qzpvntr3LVrxOFdApO2X2Lbxl5mZqZMmTYrZ869YYYv33nijjb1x6WXyZFvQ+bPP7HNw++1Wdb9atSi9QL7ZRQAzaEUf+jKO48ncfw1ZH9Zk//2j9FopREQmq2pm2HFEQ6zbsJRQ4HsCWA/R44+zsVEL5oxfxJzJa5j9Sy6z/9iNOboPs2nOYrYXqatcfistm2XTqnUFDjq8CgcdBAcdZEv0+SBzF0+Rtl8pu6TN0KF2Gt5rV6Wntm1trcHx4+2UYJ8+cNddcPHFVr39wANLeIKrr4aXX4bcAgNaMzJg27a/fk7mMB7mdobTjdr8yQu17uCKmQ9RPmW/WS7h5dVpmj9/++c1r2bAqlU2vi/v1Ff+ek4PPrjzTNWCNZ8K26YgVXTNWlb+tJwFfQaxcGNnFtB4+2VjY+Zf04Ql7LXDw/asuo7me2dz4oEV2O+wbbQ6JIODDoJmzSqQkeGDZF3ySNkeqw4dYOVKmDHDj2qcLYvzwgvw7ru2JM4hh1gtrC5dbEzWDolQRsbOCVU+v9GE9zmTN7mQyWRSk9Vcw/PcyFPsIX8W+9h0l0o9ViKZWr/+JPbck50udetC7dp2qVVr+6VmzRjVXKtdG1avjnz7ihXtjFr+9aGqVoX+/bcnTvl6mxTYQDVWV6rP6l59WF7/YJYu2MyyJTksXSosXVWBpWursDS7Bsu21GYpe7KJKju8ZCU2/ZVeNWIh+/29A/u1r0vztrux336w225l3w3OxVKk7VdKJlaLF9uA5vvus14K5/IsXWolOIYNg6+/tv8tNWrAEUdYL9Z+L91Eg60LqM56KrLF/plQi3nsw8+0YAJHMp+mALRlEhfyJpcykJoEg9ObNIHffgvt/SW6VEqsGjTI1NNPn8SyZexwWbeu+MfVqLFjolWtmuU0xV0qV7bkv3x5S8x2uH5uN8pvWI2g5FKOXMqxjYxCr2+hItlUKfqSUYN1dZqxelNlVq8V1uhurKYWa6jJtiJOcGSQw57lV1Gvylr2rJ5NvdpbqFc3l70bQeMPn6fxmuk0ZgF1WLG9YJh/T1wSSutTge+9Zz/PPjvcOFziqVcPbrjBLkuW2BisL76Ab76xuljrtj5Z6OOEXJown7ZM5p88wQmMoQW/7LiRzy5KKw0bWgdPQdnZtijv6tXw55/2s+Al7/Y1a+z677/bMKS8S3Z2aRYbH16m91GerdtTq23Z1FiTQ62KG9lLF9CSmdRiNTVZQy1W/3W97oBH2bN5Teq1qMXudcpTrly+mXf5ZZ0Avd/1WXguraRkYjVsmFXnbtky7EhcImvQAHr0sAtY79WycvVYTl3WU50tVKQ666nBOhqxkMps3vlJ8k4bRjr+xKW8KlXs41DWUmVbt1qClZdobdtm40Zzcna8ntP+6L96pTLYFvRN5RZ6vUL+JCq4lGfb9hfN35PU9Awbp1VQkybw92aRvYm870Npx2k5l8RSLrFascJ6IW6/PexIXLIRgXosox7LIntAwTEpzkVRhQp2KXns0delf/K/xljlS6oK9iQ9+GDhM/pK29vUs6d/R1xaSZTFP6Lmgw+sA8FPA7pdUqVK8ffnjTxu0sSTKpcYIqmCu8cedhGxz+5rr8HAgXY977aCn+eePe224rZxzu0k5Xqshg2DZs2gTZuwI3FJaeNGOyrPzt7x9quusmmFziWaP/8sfFZgtWpWMqS4RKikJMl7m5wrtZRKrFavhk8/heuu8xILrgzyn/pwLhn8+WfYETjnAil1KvCjj2zAp58GdM4551wYUiqxGjbMpkC3axd2JM4555xLRymTWG3cCKNHw1lnQbmUeVfOOeecSyYpk4L873823rhr17Ajcc4551y6SpnE6qOPbKmIDh3CjsQ555xz6SolEitVS6xOOsnq3jnnnHPOhSElEqvvv7eFl087LexInHPOOZfOUiKx+vBDq1vVpUvYkTjnnHMunUWUWInIySLys4jMEZHbCrn/EhFZLiJTg8vl0Q+1aB99BO3bQ9268XxV51wySPT2yzmXWkqsvC4iGcDzwAnAIuA7ERmhqj8W2PRtVb02BjEWa/FimDQJHnoo3q/snEt0id5+OedSTyQ9Vu2AOao6T1W3AG8BCVPUYORI++njq5xzhUjo9ss5l3oiSawaAgvz/b4ouK2gs0VkmogMFZFGUYkuAh99ZIuuH3RQvF7ROZdEErr9cs6lnkgSq8KWM9YCv38INFXVQ4BPgcGFPpFIbxGZJCKTli9fXrpIC7Fpky26fNppvuiyc65QUWu/IPptmHMu9USSWC0C8h/B7Q0szr+Bqq5U1c3Br68AbQt7IlXtr6qZqppZNwojzb/4wpayOeWUMj+Vcy41Ra39CraNahvmnEs9kSRW3wHNRaSZiFQEegAj8m8gIg3y/XoGMCt6IRZt9GgrCNqpUzxezTmXhBK2/XLOpaYSZwWqao6IXAuMBjKA11R1pojcB0xS1RFAHxE5A8gBVgGXxDDmv3zyCRxzDFSrFo9Xc84lm0Ruv5xzqUlUCw43iI/MzEydNGnSLj9+8WJo2BAefRRuuSWKgTnnYkZEJqtqZthxRENZ2zDnXHKJtP1K2srrn3xiP086Kdw4nHPOOefyJG1iNXo01KsHhxwSdiTOOeeccyYpE6tt22DMGDjxRC+z4JxzzrnEkZSJ1ZQpsHKlnwZ0zjnnXGJJysQqb3zVCSeEG4dzzjnnXH5JmViNHg2HHQZ77hl2JM4555xz2yVdYrV2LUyY4KcBnXPOOZd4ki6x+uILyMmBzp3DjsQ555xzbkdJl1iNGweVKsGRR4YdiXPOOefcjpIusRo7Fo46CqpUCTsS55xzzrkdJVVitXIlTJ0Kxx0XdiTOOeeccztLqsTqs89AFY4/PuxInHPOOed2llSJ1bhxUK0aHH542JE455xzzu0sqRKrsWPh2GOhYsWwI3HOOeec21nSJFZ//AE//uinAZ1zzjmXuJImsRo3zn56YuWcc865RJVUiVWtWtCmTdiROOecc84VLmkSq7FjoWNHyMgIOxLnnHPOucJFlFiJyMki8rOIzBGR2wq5v5KIvB3cP1FEmkYzyIULYe5cr1/lnCu9sNsv51x6KTGxEpEM4HngFKAlcL6ItCyw2WXAn6q6H/A08Gg0g/z8c/vZsWM0n9U5l+oSof1yzqWXSHqs2gFzVHWeqm4B3gK6FtimKzA4uD4U+JuISLSC/OIL2G03OPjgaD2jcy5NhN5+OefSSySJVUNgYb7fFwW3FbqNquYAa4A9ohEgWGJ19NE+vso5V2qht1/OufRSPoJtCjty013YBhHpDfQOft0sIjMieH3AaljF+BiyDrAipq9QeokWk8dTvESLBxIvphZxfr2otV9QtjYsCSXaZyfaUv39Qeq/x3i/vyaRbBRJYrUIaJTv972BxUVss0hEygM1gVUFn0hV+wP9AURkkqpmRhJkPCRaPJB4MXk8xUu0eCDxYhKRSXF+yai1X5DYbVi0+ftLfqn+HhP1/UVyKvA7oLmINBORikAPYESBbUYAFwfXuwNjVbXQIz7nnIsjb7+cc3FVYo+VquaIyLXAaCADeE1VZ4rIfcAkVR0BDADeEJE52JFej1gG7ZxzkfD2yzkXb5GcCkRVRwGjCtx2d77rm4BzSvna/Uu5fawlWjyQeDF5PMVLtHgg8WKKezwxar8g8fZttPn7S36p/h4T8v2J93g755xzzkVH0ixp45xzzjmX6OKWWInIOSIyU0RyRaTIUfwlLT8RxXh2F5ExIjI7+Fm7iO22icjU4FJw0Gs04ki45TYiiOkSEVmeb79cHsNYXhORZUVNaxfTN4h1mogcFqtYIoynk4isybdv7i5suyjG00hExonIrOD7dV0h28R7H0USU1z3066K4O/dM9in00TkaxFpHe8Yy6Kk95dvu8ODtrB7vGKLhkjeX/BZnBp8Vj+LZ3zREMFntKaIfCgiPwTv8dJ4x1gWidjGlUhV43IBDsRq2IwHMovYJgOYC+wDVAR+AFrGKJ7HgNuC67cBjxax3foY7pMS3y9wNfBScL0H8HaM/06RxHQJ0C9On5sOwGHAjCLu7wJ8jNUiag9MDDmeTsBH8dg3wes1AA4LrtcAfink7xXvfRRJTHHdTzH8ex8F1A6unxLrfRvv9xdskwGMxcapdQ875ij//WoBPwKNg9/3DDvmGLzHO/L+vwF1sQkaFcOOuxTvL+HauJIuceuxUtVZqvpzCZtFsvxEtORfxmIwcGaMXqc4ibjcRjz/BiVS1c8poqZQoCvwuppvgFoi0iDEeOJKVZeo6pTg+jpgFjtXFo/3PookpqRQ0t9bVb9W1T+DX7/B6mQljQg/z/8AhgHLYh9RdEXw/i4AhqvqgmD7VHyPCtQI/m9UD7bNiUds0ZCIbVxJEm2MVSTLT0RLPVVdAvaHA/YsYrvKIjJJRL4RkWgnX4m43Eakf4Ozgy7XoSLSqJD74yWen5lIHRl0u38sIq3i9aLBaeJDgYkF7gptHxUTE4S0n2LoMuyoOWWISEPgLOClsGOJkf2B2iIyXkQmi0ivsAOKgX7YGaPFwHTgOlXNDTekXZOIbVxhIiq3ECkR+RSoX8hdd6rqB5E8RSG37fK0xeLiKcXTNFbVxSKyDzBWRKar6txdjamAqC63ESWRvN6HwBBV3SwiV2I9asfHMKbixHv/lGQK0ERV14tIF+B9oHmsX1REqmO9Cter6tqCdxfykJjvoxJiCmU/xYqIHIclVseEHUuUPQPcqqrbYttRHpryQFvgb0AVYIKIfKOqv4QbVlSdBEzF2uh9gTEi8kUh38mElohtXFGimlipaucyPkUky09EJR4RWSoiDVR1SdBlWGgXsKouDn7OE5HxWLYcrcQqqsttxCsmVV2Z79dXgEdjGE9JovqZKav8X3hVHSUiL4hIHVWN2XpWIlIBa3CyVHV4IZvEfR+VFFMY+ylWROQQ4FXglALfjVSQCbwVJFV1gC4ikqOq74cbVtQsAlao6gZgg4h8DrTGxvGkikuBR9QGI80RkV+BA4Bvww0rconYxhUn0U4FRrL8RLTkX8biYmCnHjURqS0ilYLrdYCjsYGO0ZKIy22UGFOBc9dnYOe8wzIC6BXMCmkPrMk7xRsGEamfNwZORNph37GY/bMNXmsAMEtVnypis7juo0hiivd+ihURaQwMBy5KsV4OAFS1mao2VdWm2BjPq1MoqQJr948VkfIiUhU4gnDbs1hYgPXIISL1sElk80KNqBQSsY0rSVR7rIojImcBz2GzEkaKyFRVPUlE9gJeVdUuWsTyEzEK6RHgHRG5DPvgnRPEmQlcqaqXY+elXxaRXKzhf0RVo5ZYFfV+JcTlNiKMqY+InIENgFyFzRKMCREZgs0gqyMii4B7gApBrC9hM5W6AHOAjdjRWcxEEE934CoRyQGygR4xToSPBi4CpovI1OC2O4DG+WKK6z6KMKZ476ddEsHf+25szOMLQZ6Yowm4KGxRInh/Sa2k96eqs0Tkv8A0IBf7X1Rs6YlEE8Hf8H5gkIhMx06Z3ZpkPcOJ2MYVyyuvO+ecc85FSaKdCnTOOeecS1qeWDnnnHPORYknVs4555xzUeKJlXPOOedclHhi5ZxzzjkXJZ5YOeecc85FiSdWzjnnnHNR4omVc84551yU/D8bG2sBKoqLogAAAABJRU5ErkJggg==\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": [
    "## Let's implement SLIM BPR "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In order to implement a SLIM BPR we need to:\n",
    "#### Randomly sample the triplets (user, positive_item, negative_item)\n",
    "#### Compute the score of each triplet\n",
    "#### Update the similarity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from urllib.request import urlretrieve\n",
    "import zipfile, os\n",
    "\n",
    "# If file exists, skip the download\n",
    "data_file_path = \"data/Movielens_10M/\"\n",
    "data_file_name = data_file_path + \"movielens_10m.zip\"\n",
    "\n",
    "# If directory does not exist, create\n",
    "if not os.path.exists(data_file_path):\n",
    "    os.makedirs(data_file_path)\n",
    "\n",
    "if not os.path.exists(data_file_name):\n",
    "    urlretrieve (\"http://files.grouplens.org/datasets/movielens/ml-10m.zip\", data_file_name)\n",
    "\n",
    "dataFile = zipfile.ZipFile(data_file_name)\n",
    "URM_path = dataFile.extract(\"ml-10M100K/ratings.dat\", path=\"data/Movielens_10M\")\n",
    "URM_file = open(URM_path, 'r')\n",
    "\n",
    "\n",
    "def rowSplit (rowString):\n",
    "    \n",
    "    split = rowString.split(\"::\")\n",
    "    split[3] = split[3].replace(\"\\n\",\"\")\n",
    "    \n",
    "    split[0] = int(split[0])\n",
    "    split[1] = int(split[1])\n",
    "    split[2] = float(split[2])\n",
    "    split[3] = int(split[3])\n",
    "    \n",
    "    result = tuple(split)\n",
    "    \n",
    "    return result\n",
    "\n",
    "\n",
    "URM_file.seek(0)\n",
    "URM_tuples = []\n",
    "\n",
    "for line in URM_file:\n",
    "   URM_tuples.append(rowSplit (line))\n",
    "\n",
    "userList, itemList, ratingList, timestampList = zip(*URM_tuples)\n",
    "\n",
    "userList = list(userList)\n",
    "itemList = list(itemList)\n",
    "ratingList = list(ratingList)\n",
    "timestampList = list(timestampList)\n",
    "\n",
    "import scipy.sparse as sps\n",
    "\n",
    "URM_all = sps.coo_matrix((ratingList, (userList, itemList)))\n",
    "URM_all = URM_all.tocsr()\n",
    "\n",
    "\n",
    "\n",
    "from Data_manager.split_functions.split_train_validation_random_holdout import split_train_in_two_percentage_global_sample\n",
    "\n",
    "\n",
    "URM_train, URM_test = split_train_in_two_percentage_global_sample(URM_all, train_percentage = 0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<71568x5000 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 6769221 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "URM_train = URM_train[:,0:5000]\n",
    "URM_test = URM_test[:,0:5000]\n",
    "URM_train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 - Sampling\n",
    "\n",
    "#### Create a mask of positive interactions. How to build it depends on the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<71568x5000 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 3916039 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 9,
     "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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_users = URM_mask.shape[0]\n",
    "n_items = URM_mask.shape[1]\n",
    "\n",
    "\n",
    "# Extract users having at least one interaction to choose from\n",
    "eligibleUsers = []\n",
    "\n",
    "for user_id in range(n_users):\n",
    "\n",
    "    start_pos = URM_mask.indptr[user_id]\n",
    "    end_pos = URM_mask.indptr[user_id+1]\n",
    "\n",
    "    if len(URM_mask.indices[start_pos:end_pos]) > 0:\n",
    "        eligibleUsers.append(user_id)\n",
    "                \n",
    "                \n",
    "\n",
    "def sampleTriplet():\n",
    "    \n",
    "    # By randomly selecting a user in this way we could end up \n",
    "    # with a user with no interactions\n",
    "    #user_id = np.random.randint(0, n_users)\n",
    "    \n",
    "    user_id = np.random.choice(eligibleUsers)\n",
    "    \n",
    "    # Get user seen items and choose one\n",
    "    userSeenItems = URM_mask[user_id,:].indices\n",
    "    pos_item_id = np.random.choice(userSeenItems)\n",
    "\n",
    "    negItemSelected = False\n",
    "\n",
    "    # It's faster to just try again then to build a mapping of the non-seen items\n",
    "    while (not negItemSelected):\n",
    "        neg_item_id = np.random.randint(0, n_items)\n",
    "\n",
    "        if (neg_item_id not in userSeenItems):\n",
    "            \n",
    "            negItemSelected = True\n",
    "\n",
    "    return user_id, pos_item_id, neg_item_id\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8288, 105, 1630)\n",
      "(2982, 1094, 3317)\n",
      "(26453, 653, 358)\n",
      "(11989, 1148, 1244)\n",
      "(55371, 3358, 3572)\n",
      "(33429, 3996, 1101)\n",
      "(55792, 50, 3366)\n",
      "(57923, 2081, 2453)\n",
      "(26150, 498, 1001)\n",
      "(38236, 110, 1682)\n"
     ]
    }
   ],
   "source": [
    "for _ in range(10):\n",
    "    print(sampleTriplet())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Computing prediction\n",
    "\n",
    "#### The prediction depends on the model: SLIM, Matrix Factorization... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We have to initialize our model. In case of SLIM it works best to initialize S as zero, in case of MF you cannot because of how the gradient is computed and you have to initialize at random. Here we initialize SLIM at random just so that we have some numbers to show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#similarity_matrix = np.zeros((n_items,n_items))\n",
    "\n",
    "similarity_matrix = np.random.random((n_items,n_items))\n",
    "similarity_matrix[np.arange(n_items),np.arange(n_items)] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_id, positive_item_id, negative_item_id = sampleTriplet()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "372"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "positive_item_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "933"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "negative_item_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  22,   32,  110,  372,  380,  441,  593,  608,  858, 1042, 1193,\n",
       "       1513, 1721, 1747, 2144, 2268, 2762, 2858, 2959, 3863, 4865, 4963])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "userSeenItems = URM_mask[user_id,:].indices\n",
    "userSeenItems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_i is 10.36, x_j is 9.61\n"
     ]
    }
   ],
   "source": [
    "x_i = similarity_matrix[positive_item_id, userSeenItems].sum()\n",
    "x_j = similarity_matrix[negative_item_id, userSeenItems].sum()\n",
    "\n",
    "print(\"x_i is {:.2f}, x_j is {:.2f}\".format(x_i, x_j))"
   ]
  },
  {
   "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": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7458571825185416"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_ij = x_i - x_j\n",
    "x_ij"
   ]
  },
  {
   "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": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.32172466929646293"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient = 1 / (1 + np.exp(x_ij))\n",
    "gradient"
   ]
  },
  {
   "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": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = 1e-3\n",
    "\n",
    "similarity_matrix[positive_item_id, userSeenItems] += learning_rate * gradient\n",
    "similarity_matrix[positive_item_id, positive_item_id] = 0\n",
    "\n",
    "similarity_matrix[negative_item_id, userSeenItems] -= learning_rate * gradient\n",
    "similarity_matrix[negative_item_id, negative_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": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_i is 10.37, x_j is 9.61\n"
     ]
    }
   ],
   "source": [
    "x_i = similarity_matrix[positive_item_id, userSeenItems].sum()\n",
    "x_j = similarity_matrix[negative_item_id, userSeenItems].sum()\n",
    "\n",
    "print(\"x_i is {:.2f}, x_j is {:.2f}\".format(x_i, x_j))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5 - Write the iterative epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "def epochIteration():\n",
    "\n",
    "    # Get number of available interactions\n",
    "    numPositiveIteractions = int(URM_mask.nnz*0.01)\n",
    "\n",
    "    start_time_epoch = time.time()\n",
    "    start_time_batch = time.time()\n",
    "\n",
    "    # Uniform user sampling without replacement\n",
    "    for num_sample in range(numPositiveIteractions):\n",
    "\n",
    "        # Sample\n",
    "        user_id, positive_item_id, negative_item_id = sampleTriplet()\n",
    "        \n",
    "        userSeenItems = URM_mask[user_id,:].indices\n",
    "        \n",
    "        # Prediction\n",
    "        x_i = similarity_matrix[positive_item_id, userSeenItems].sum()\n",
    "        x_j = similarity_matrix[negative_item_id, userSeenItems].sum()\n",
    "        \n",
    "        # Gradient\n",
    "        x_ij = x_i - x_j\n",
    "\n",
    "        gradient = 1 / (1 + np.exp(x_ij))\n",
    "        \n",
    "        # Update\n",
    "        similarity_matrix[positive_item_id, userSeenItems] += learning_rate * gradient\n",
    "        similarity_matrix[positive_item_id, positive_item_id] = 0\n",
    "\n",
    "        similarity_matrix[negative_item_id, userSeenItems] -= learning_rate * gradient\n",
    "        similarity_matrix[negative_item_id, negative_item_id] = 0\n",
    "        \n",
    "\n",
    "        if(time.time() - start_time_batch >= 30 or num_sample == numPositiveIteractions-1):\n",
    "            print(\"Processed {} ( {:.2f}% ) in {:.2f} seconds. Sample per second: {:.0f}\".format(\n",
    "                num_sample,\n",
    "                100.0* float(num_sample)/numPositiveIteractions,\n",
    "                time.time() - start_time_batch,\n",
    "                float(num_sample) / (time.time() - start_time_epoch)))\n",
    "\n",
    "\n",
    "            start_time_batch = time.time()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processed 5378 ( 13.73% ) in 30.01 seconds. Sample per second: 179\n",
      "Processed 10458 ( 26.71% ) in 30.01 seconds. Sample per second: 174\n",
      "Processed 16640 ( 42.49% ) in 30.00 seconds. Sample per second: 185\n",
      "Processed 22990 ( 58.71% ) in 30.00 seconds. Sample per second: 192\n",
      "Processed 28602 ( 73.04% ) in 30.01 seconds. Sample per second: 191\n",
      "Processed 33604 ( 85.81% ) in 30.01 seconds. Sample per second: 187\n",
      "Processed 38721 ( 98.88% ) in 30.01 seconds. Sample per second: 184\n",
      "Processed 39159 ( 100.00% ) in 2.33 seconds. Sample per second: 184\n"
     ]
    }
   ],
   "source": [
    "epochIteration()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sps\n",
    "\n",
    "\n",
    "def similarityMatrixTopK(item_weights, forceSparseOutput = True, k=100, verbose = False, inplace=True):\n",
    "    \"\"\"\n",
    "    The function selects the TopK most similar elements, column-wise\n",
    "\n",
    "    :param item_weights:\n",
    "    :param forceSparseOutput:\n",
    "    :param k:\n",
    "    :param verbose:\n",
    "    :param inplace: Default True, WARNING matrix will be modified\n",
    "    :return:\n",
    "    \"\"\"\n",
    "\n",
    "    assert (item_weights.shape[0] == item_weights.shape[1]), \"selectTopK: ItemWeights is not a square matrix\"\n",
    "\n",
    "    start_time = time.time()\n",
    "\n",
    "    if verbose:\n",
    "        print(\"Generating topK matrix\")\n",
    "\n",
    "    nitems = item_weights.shape[1]\n",
    "    k = min(k, nitems)\n",
    "\n",
    "    # for each column, keep only the top-k scored items\n",
    "    sparse_weights = not isinstance(item_weights, np.ndarray)\n",
    "\n",
    "    if not sparse_weights:\n",
    "\n",
    "        idx_sorted = np.argsort(item_weights, axis=0)  # sort data inside each column\n",
    "\n",
    "        if inplace:\n",
    "            W = item_weights\n",
    "        else:\n",
    "            W = item_weights.copy()\n",
    "\n",
    "        # index of the items that don't belong to the top-k similar items of each column\n",
    "        not_top_k = idx_sorted[:-k, :]\n",
    "        # use numpy fancy indexing to zero-out the values in sim without using a for loop\n",
    "        W[not_top_k, np.arange(nitems)] = 0.0\n",
    "\n",
    "        if forceSparseOutput:\n",
    "            W_sparse = sps.csr_matrix(W, shape=(nitems, nitems))\n",
    "\n",
    "            if verbose:\n",
    "                print(\"Sparse TopK matrix generated in {:.2f} seconds\".format(time.time() - start_time))\n",
    "\n",
    "            return W_sparse\n",
    "\n",
    "        if verbose:\n",
    "            print(\"Dense TopK matrix generated in {:.2f} seconds\".format(time.time()-start_time))\n",
    "\n",
    "        return W\n",
    "\n",
    "    else:\n",
    "        # iterate over each column and keep only the top-k similar items\n",
    "        data, rows_indices, cols_indptr = [], [], []\n",
    "\n",
    "        item_weights = check_matrix(item_weights, format='csc', dtype=np.float32)\n",
    "\n",
    "        for item_idx in range(nitems):\n",
    "\n",
    "            cols_indptr.append(len(data))\n",
    "\n",
    "            start_position = item_weights.indptr[item_idx]\n",
    "            end_position = item_weights.indptr[item_idx+1]\n",
    "\n",
    "            column_data = item_weights.data[start_position:end_position]\n",
    "            column_row_index = item_weights.indices[start_position:end_position]\n",
    "\n",
    "            non_zero_data = column_data!=0\n",
    "\n",
    "            idx_sorted = np.argsort(column_data[non_zero_data])  # sort by column\n",
    "            top_k_idx = idx_sorted[-k:]\n",
    "\n",
    "            data.extend(column_data[non_zero_data][top_k_idx])\n",
    "            rows_indices.extend(column_row_index[non_zero_data][top_k_idx])\n",
    "\n",
    "\n",
    "        cols_indptr.append(len(data))\n",
    "\n",
    "        # During testing CSR is faster\n",
    "        W_sparse = sps.csc_matrix((data, rows_indices, cols_indptr), shape=(nitems, nitems), dtype=np.float32)\n",
    "        W_sparse = W_sparse.tocsr()\n",
    "\n",
    "        if verbose:\n",
    "            print(\"Sparse TopK matrix generated in {:.2f} seconds\".format(time.time() - start_time))\n",
    "\n",
    "        return W_sparse\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SLIM_BPR_Recommender(object):\n",
    "    \"\"\" SLIM_BPR recommender with cosine similarity and no shrinkage\"\"\"\n",
    "\n",
    "    def __init__(self, URM):\n",
    "        self.URM = URM\n",
    "        \n",
    "        self.similarity_matrix = np.zeros((n_items,n_items))\n",
    "        \n",
    "        self.URM_mask = self.URM.copy()\n",
    "        self.URM_mask.data[self.URM_mask.data <= 3] = 0\n",
    "        self.URM_mask.eliminate_zeros()\n",
    "        \n",
    "        self.n_users = URM_mask.shape[0]\n",
    "        self.n_items = URM_mask.shape[1]\n",
    "\n",
    "\n",
    "        # Extract users having at least one interaction to choose from\n",
    "        self.eligibleUsers = []\n",
    "\n",
    "        for user_id in range(n_users):\n",
    "\n",
    "            start_pos = self.URM_mask.indptr[user_id]\n",
    "            end_pos = self.URM_mask.indptr[user_id+1]\n",
    "\n",
    "            if len(self.URM_mask.indices[start_pos:end_pos]) > 0:\n",
    "                self.eligibleUsers.append(user_id)\n",
    "\n",
    "\n",
    "\n",
    "    def sampleTriplet(self):\n",
    "\n",
    "        # By randomly selecting a user in this way we could end up \n",
    "        # with a user with no interactions\n",
    "        #user_id = np.random.randint(0, n_users)\n",
    "\n",
    "        user_id = np.random.choice(self.eligibleUsers)\n",
    "\n",
    "        # Get user seen items and choose one\n",
    "        userSeenItems = URM_mask[user_id,:].indices\n",
    "        pos_item_id = np.random.choice(userSeenItems)\n",
    "\n",
    "        negItemSelected = False\n",
    "\n",
    "        # It's faster to just try again then to build a mapping of the non-seen items\n",
    "        while (not negItemSelected):\n",
    "            neg_item_id = np.random.randint(0, n_items)\n",
    "\n",
    "            if (neg_item_id not in userSeenItems):\n",
    "\n",
    "                negItemSelected = True\n",
    "\n",
    "        return user_id, pos_item_id, neg_item_id\n",
    "\n",
    "\n",
    "        \n",
    "    def epochIteration(self):\n",
    "\n",
    "        # Get number of available interactions\n",
    "        numPositiveIteractions = int(self.URM_mask.nnz*0.01)\n",
    "\n",
    "        start_time_epoch = time.time()\n",
    "        start_time_batch = time.time()\n",
    "\n",
    "        # Uniform user sampling without replacement\n",
    "        for num_sample in range(numPositiveIteractions):\n",
    "\n",
    "            # Sample\n",
    "            user_id, positive_item_id, negative_item_id = self.sampleTriplet()\n",
    "\n",
    "            userSeenItems = self.URM_mask[user_id,:].indices\n",
    "\n",
    "            # Prediction\n",
    "            x_i = self.similarity_matrix[positive_item_id, userSeenItems].sum()\n",
    "            x_j = self.similarity_matrix[negative_item_id, userSeenItems].sum()\n",
    "\n",
    "            # Gradient\n",
    "            x_ij = x_i - x_j\n",
    "\n",
    "            gradient = 1 / (1 + np.exp(x_ij))\n",
    "\n",
    "            # Update\n",
    "            self.similarity_matrix[positive_item_id, userSeenItems] += learning_rate * gradient\n",
    "            self.similarity_matrix[positive_item_id, positive_item_id] = 0\n",
    "\n",
    "            self.similarity_matrix[negative_item_id, userSeenItems] -= learning_rate * gradient\n",
    "            self.similarity_matrix[negative_item_id, negative_item_id] = 0\n",
    "\n",
    "\n",
    "            if(time.time() - start_time_batch >= 30 or num_sample == numPositiveIteractions-1):\n",
    "                print(\"Processed {} ( {:.2f}% ) in {:.2f} seconds. Sample per second: {:.0f}\".format(\n",
    "                    num_sample,\n",
    "                    100.0* float(num_sample)/numPositiveIteractions,\n",
    "                    time.time() - start_time_batch,\n",
    "                    float(num_sample) / (time.time() - start_time_epoch)))\n",
    "\n",
    "                start_time_batch = time.time()\n",
    "\n",
    "                \n",
    "    def fit(self, learning_rate = 0.01, epochs = 10):\n",
    " \n",
    "        self.learning_rate = learning_rate\n",
    "        self.epochs = epochs\n",
    "\n",
    "        for numEpoch in range(self.epochs):\n",
    "            self.epochIteration()\n",
    "            \n",
    "        self.similarity_matrix = self.similarity_matrix.T\n",
    "        \n",
    "        self.similarity_matrix = similarityMatrixTopK(self.similarity_matrix, k=100)\n",
    "        \n",
    "        \n",
    "        \n",
    "    def recommend(self, user_id, at=None, exclude_seen=True):\n",
    "        # compute the scores using the dot product\n",
    "        user_profile = self.URM[user_id]\n",
    "        scores = user_profile.dot(self.similarity_matrix).toarray().ravel()\n",
    "\n",
    "        if exclude_seen:\n",
    "            scores = self.filter_seen(user_id, scores)\n",
    "\n",
    "        # rank items\n",
    "        ranking = scores.argsort()[::-1]\n",
    "            \n",
    "        return ranking[:at]\n",
    "    \n",
    "    \n",
    "    def filter_seen(self, user_id, scores):\n",
    "\n",
    "        start_pos = self.URM.indptr[user_id]\n",
    "        end_pos = self.URM.indptr[user_id+1]\n",
    "\n",
    "        user_profile = self.URM.indices[start_pos:end_pos]\n",
    "        \n",
    "        scores[user_profile] = -np.inf\n",
    "\n",
    "        return scores  \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processed 5494 ( 14.03% ) in 30.01 seconds. Sample per second: 183\n",
      "Processed 11473 ( 29.30% ) in 30.01 seconds. Sample per second: 191\n",
      "Processed 17290 ( 44.15% ) in 30.02 seconds. Sample per second: 192\n",
      "Processed 23027 ( 58.80% ) in 30.01 seconds. Sample per second: 192\n",
      "Processed 28604 ( 73.04% ) in 30.01 seconds. Sample per second: 191\n",
      "Processed 33575 ( 85.74% ) in 30.01 seconds. Sample per second: 186\n",
      "Processed 38543 ( 98.42% ) in 30.01 seconds. Sample per second: 183\n",
      "Processed 39159 ( 100.00% ) in 3.91 seconds. Sample per second: 183\n"
     ]
    }
   ],
   "source": [
    "recommender = SLIM_BPR_Recommender(URM_train)\n",
    "recommender.fit(epochs=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluated user 0 of 71568\n",
      "Evaluated user 10000 of 71568\n",
      "Evaluated user 20000 of 71568\n",
      "Evaluated user 30000 of 71568\n",
      "Evaluated user 40000 of 71568\n",
      "Evaluated user 50000 of 71568\n",
      "Evaluated user 60000 of 71568\n",
      "Evaluated user 70000 of 71568\n",
      "Recommender performance is: Precision = 0.2621, Recall = 0.0896, MAP = 0.1989\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'precision': 0.26210298398320586,\n",
       " 'recall': 0.08955017906833086,\n",
       " 'MAP': 0.19894397921779672}"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Notebooks_utils.evaluation_function import evaluate_algorithm\n",
    "\n",
    "evaluate_algorithm(URM_test, recommender)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reasonable performance\n",
    "#### Note that we ran just one epoch on a micro-subsample of the data. Being a machine learning approach, more time is needed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's load the cython impementation and run some serious learning. Here I show just two epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLIM_BPR_Recommender: URM Detected 1692 (2.36 %) cold users.\n",
      "SLIM_BPR_Recommender: URM Detected 96 (1.92 %) cold items.\n",
      "Unable to read memory status: list index out of range\n",
      "SLIM_BPR_Recommender: Automatic selection of fastest train mode. Unable to get current RAM status, you may be using a non-Linux operating system. Using dense matrix.\n",
      "Processed 71568 ( 100.00% ) in 0.68 seconds. BPR loss is 2.16E-03. Sample per second: 105931\n",
      "SLIM_BPR_Recommender: Epoch 1 of 2. Elapsed time 0.45 sec\n",
      "Processed 71568 ( 100.00% ) in 1.13 seconds. BPR loss is 1.37E-02. Sample per second: 63589\n",
      "SLIM_BPR_Recommender: Epoch 2 of 2. Elapsed time 0.90 sec\n",
      "SLIM_BPR_Recommender: Terminating at epoch 2. Elapsed time 2.05 sec\n",
      "Deallocating Cython objects\n"
     ]
    }
   ],
   "source": [
    "from SLIM_BPR.Cython.SLIM_BPR_Cython import SLIM_BPR_Cython\n",
    "\n",
    "\n",
    "\n",
    "recommender = SLIM_BPR_Cython(URM_train, recompile_cython=False)\n",
    "\n",
    "recommender.fit(epochs=2, batch_size=1, sgd_mode='sgd', learning_rate=1e-4, positive_threshold_BPR=4)\n",
    " \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EvaluatorHoldout: Processed 59510 ( 85.33% ) in 30.02 sec. Users per second: 1983\n",
      "EvaluatorHoldout: Processed 69739 ( 100.00% ) in 36.05 sec. Users per second: 1934\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "({5: {'ROC_AUC': 0.37906933948961963,\n",
       "   'PRECISION': 0.2629948809131884,\n",
       "   'PRECISION_RECALL_MIN_DEN': 0.2686786924581338,\n",
       "   'RECALL': 0.08816346989303116,\n",
       "   'MAP': 0.19634895825869142,\n",
       "   'MRR': 0.43688753781957185,\n",
       "   'NDCG': 0.139975255418048,\n",
       "   'F1': 0.1320574675907752,\n",
       "   'HIT_RATE': 1.3149744045655947,\n",
       "   'ARHR': 0.6473783201173414,\n",
       "   'RMSE': 3.351348491785974,\n",
       "   'NOVELTY': 0.008235469902706695,\n",
       "   'AVERAGE_POPULARITY': 0.8221555650858704,\n",
       "   'DIVERSITY_MEAN_INTER_LIST': 0.6993780325105289,\n",
       "   'DIVERSITY_HERFINDAHL': 0.939873600800755,\n",
       "   'COVERAGE_ITEM': 0.0786,\n",
       "   'COVERAGE_USER': 0.9744438855354348,\n",
       "   'DIVERSITY_GINI': 0.0483961992295712,\n",
       "   'SHANNON_ENTROPY': 4.566429947440342}},\n",
       " 'CUTOFF: 5 - ROC_AUC: 0.3790693, PRECISION: 0.2629949, PRECISION_RECALL_MIN_DEN: 0.2686787, RECALL: 0.0881635, MAP: 0.1963490, MRR: 0.4368875, NDCG: 0.1399753, F1: 0.1320575, HIT_RATE: 1.3149744, ARHR: 0.6473783, RMSE: 3.3513485, NOVELTY: 0.0082355, AVERAGE_POPULARITY: 0.8221556, DIVERSITY_MEAN_INTER_LIST: 0.6993780, DIVERSITY_HERFINDAHL: 0.9398736, COVERAGE_ITEM: 0.0786000, COVERAGE_USER: 0.9744439, DIVERSITY_GINI: 0.0483962, SHANNON_ENTROPY: 4.5664299, \\n')"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from Base.Evaluation.Evaluator import EvaluatorHoldout\n",
    "\n",
    "evaluator_validation = EvaluatorHoldout(URM_test, cutoff_list=[5])\n",
    "\n",
    "\n",
    "evaluator_validation.evaluateRecommender(recommender)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compared to content and collaborative k-nearest neighbor SLIM improves map significantly (on this dataset)"
   ]
  },
  {
   "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": {
    "collapsed": true
   },
   "source": [
    "## Let's now move to SLIM ElasticNet \n",
    "\n",
    "#### ElasticNet optimizes prediction error and l1, l2 norms over the similarity matrix to ensure sparsity and counter overfitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Key idea:\n",
    "##### Learn every item similarity column independently, in such a way that those similarities will allow the model to predict the interaction score (rating) for that item using all other item iteractions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import scipy.sparse as sps\n",
    "\n",
    "from sklearn.linear_model import ElasticNet\n",
    "\n",
    "import time, sys\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "class SLIMElasticNetRecommender(object):\n",
    "    \"\"\"\n",
    "    Train a Sparse Linear Methods (SLIM) item similarity model.\n",
    "    NOTE: ElasticNet solver is parallel, a single intance of SLIM_ElasticNet will\n",
    "          make use of half the cores available\n",
    "\n",
    "    See:\n",
    "        Efficient Top-N Recommendation by Linear Regression,\n",
    "        M. Levy and K. Jack, LSRS workshop at RecSys 2013.\n",
    "        https://www.slideshare.net/MarkLevy/efficient-slides\n",
    "\n",
    "        SLIM: Sparse linear methods for top-n recommender systems,\n",
    "        X. Ning and G. Karypis, ICDM 2011.\n",
    "        http://glaros.dtc.umn.edu/gkhome/fetch/papers/SLIM2011icdm.pdf\n",
    "    \"\"\"\n",
    "\n",
    "    RECOMMENDER_NAME = \"SLIMElasticNetRecommender\"\n",
    "\n",
    "    def __init__(self, URM_train):\n",
    "\n",
    "        super(SLIMElasticNetRecommender, self).__init__()\n",
    "\n",
    "        self.URM_train = URM_train\n",
    "\n",
    "\n",
    "    def fit(self, l1_penalty=0.1, l2_penalty=0.1, positive_only=True, topK = 100):\n",
    "\n",
    "        self.l1_penalty = l1_penalty\n",
    "        self.l2_penalty = l2_penalty\n",
    "        self.positive_only = positive_only\n",
    "        self.topK = topK\n",
    "\n",
    "        if self.l1_penalty + self.l2_penalty != 0:\n",
    "            self.l1_ratio = self.l1_penalty / (self.l1_penalty + self.l2_penalty)\n",
    "        else:\n",
    "            print(\"SLIM_ElasticNet: l1_penalty+l2_penalty cannot be equal to zero, setting the ratio l1/(l1+l2) to 1.0\")\n",
    "            self.l1_ratio = 1.0\n",
    "\n",
    "\n",
    "\n",
    "        # initialize the ElasticNet model\n",
    "        self.model = ElasticNet(alpha=1.0,\n",
    "                                l1_ratio=self.l1_ratio,\n",
    "                                positive=self.positive_only,\n",
    "                                fit_intercept=False,\n",
    "                                copy_X=False,\n",
    "                                precompute=True,\n",
    "                                selection='random',\n",
    "                                max_iter=100,\n",
    "                                tol=1e-4)\n",
    "\n",
    "\n",
    "        URM_train = sps.csc_matrix(self.URM_train)\n",
    "\n",
    "        n_items = URM_train.shape[1]\n",
    "\n",
    "\n",
    "        # Use array as it reduces memory requirements compared to lists\n",
    "        dataBlock = 10000000\n",
    "\n",
    "        rows = np.zeros(dataBlock, dtype=np.int32)\n",
    "        cols = np.zeros(dataBlock, dtype=np.int32)\n",
    "        values = np.zeros(dataBlock, dtype=np.float32)\n",
    "\n",
    "        numCells = 0\n",
    "\n",
    "\n",
    "        start_time = time.time()\n",
    "        start_time_printBatch = start_time\n",
    "\n",
    "        # fit each item's factors sequentially (not in parallel)\n",
    "        for currentItem in range(n_items):\n",
    "\n",
    "            # get the target column\n",
    "            y = URM_train[:, currentItem].toarray()\n",
    "\n",
    "            # set the j-th column of X to zero\n",
    "            start_pos = URM_train.indptr[currentItem]\n",
    "            end_pos = URM_train.indptr[currentItem + 1]\n",
    "\n",
    "            current_item_data_backup = URM_train.data[start_pos: end_pos].copy()\n",
    "            URM_train.data[start_pos: end_pos] = 0.0\n",
    "\n",
    "\n",
    "\n",
    "            # fit one ElasticNet model per column\n",
    "            self.model.fit(URM_train, y)\n",
    "\n",
    "            # self.model.coef_ contains the coefficient of the ElasticNet model\n",
    "            # let's keep only the non-zero values\n",
    "\n",
    "            # Select topK values\n",
    "            # Sorting is done in three steps. Faster then plain np.argsort for higher number of items\n",
    "            # - Partition the data to extract the set of relevant items\n",
    "            # - Sort only the relevant items\n",
    "            # - Get the original item index\n",
    "\n",
    "            nonzero_model_coef_index = self.model.sparse_coef_.indices\n",
    "            nonzero_model_coef_value = self.model.sparse_coef_.data\n",
    "\n",
    "\n",
    "            local_topK = min(len(nonzero_model_coef_value)-1, self.topK)\n",
    "\n",
    "            relevant_items_partition = (-nonzero_model_coef_value).argpartition(local_topK)[0:local_topK]\n",
    "            relevant_items_partition_sorting = np.argsort(-nonzero_model_coef_value[relevant_items_partition])\n",
    "            ranking = relevant_items_partition[relevant_items_partition_sorting]\n",
    "\n",
    "\n",
    "            for index in range(len(ranking)):\n",
    "\n",
    "                if numCells == len(rows):\n",
    "                    rows = np.concatenate((rows, np.zeros(dataBlock, dtype=np.int32)))\n",
    "                    cols = np.concatenate((cols, np.zeros(dataBlock, dtype=np.int32)))\n",
    "                    values = np.concatenate((values, np.zeros(dataBlock, dtype=np.float32)))\n",
    "\n",
    "\n",
    "                rows[numCells] = nonzero_model_coef_index[ranking[index]]\n",
    "                cols[numCells] = currentItem\n",
    "                values[numCells] = nonzero_model_coef_value[ranking[index]]\n",
    "\n",
    "                numCells += 1\n",
    "\n",
    "\n",
    "            # finally, replace the original values of the j-th column\n",
    "            URM_train.data[start_pos:end_pos] = current_item_data_backup\n",
    "\n",
    "\n",
    "            if time.time() - start_time_printBatch > 300 or currentItem == n_items-1:\n",
    "                print(\"Processed {} ( {:.2f}% ) in {:.2f} minutes. Items per second: {:.0f}\".format(\n",
    "                                  currentItem+1,\n",
    "                                  100.0* float(currentItem+1)/n_items,\n",
    "                                  (time.time()-start_time)/60,\n",
    "                                  float(currentItem)/(time.time()-start_time)))\n",
    "                sys.stdout.flush()\n",
    "                sys.stderr.flush()\n",
    "\n",
    "                start_time_printBatch = time.time()\n",
    "\n",
    "\n",
    "        # generate the sparse weight matrix\n",
    "        self.W_sparse = sps.csr_matrix((values[:numCells], (rows[:numCells], cols[:numCells])),\n",
    "                                       shape=(n_items, n_items), dtype=np.float32)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "        \n",
    "    def recommend(self, user_id, at=None, exclude_seen=True):\n",
    "        # compute the scores using the dot product\n",
    "        user_profile = self.URM_train[user_id]\n",
    "        scores = user_profile.dot(self.W_sparse).toarray().ravel()\n",
    "\n",
    "        if exclude_seen:\n",
    "            scores = self.filter_seen(user_id, scores)\n",
    "\n",
    "        # rank items\n",
    "        ranking = scores.argsort()[::-1]\n",
    "            \n",
    "        return ranking[:at]\n",
    "    \n",
    "    \n",
    "    def filter_seen(self, user_id, scores):\n",
    "\n",
    "        start_pos = self.URM_train.indptr[user_id]\n",
    "        end_pos = self.URM_train.indptr[user_id+1]\n",
    "\n",
    "        user_profile = self.URM_train.indices[start_pos:end_pos]\n",
    "        \n",
    "        scores[user_profile] = -np.inf\n",
    "\n",
    "        return scores  \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Deallocating Cython objects\n",
      "Processed 599 ( 11.98% ) in 5.00 minutes. Items per second: 2\n",
      "Processed 1141 ( 22.82% ) in 10.00 minutes. Items per second: 2\n",
      "Processed 1537 ( 30.74% ) in 15.01 minutes. Items per second: 2\n",
      "Processed 1842 ( 36.84% ) in 20.01 minutes. Items per second: 2\n",
      "Processed 2198 ( 43.96% ) in 25.02 minutes. Items per second: 1\n",
      "Processed 2588 ( 51.76% ) in 30.03 minutes. Items per second: 1\n",
      "Processed 2960 ( 59.20% ) in 35.04 minutes. Items per second: 1\n",
      "Processed 3340 ( 66.80% ) in 40.05 minutes. Items per second: 1\n",
      "Processed 3726 ( 74.52% ) in 45.05 minutes. Items per second: 1\n",
      "Processed 4074 ( 81.48% ) in 50.06 minutes. Items per second: 1\n",
      "Processed 4445 ( 88.90% ) in 55.06 minutes. Items per second: 1\n",
      "Processed 4913 ( 98.26% ) in 60.07 minutes. Items per second: 1\n",
      "Processed 5000 ( 100.00% ) in 60.64 minutes. Items per second: 1\n"
     ]
    }
   ],
   "source": [
    "recommender = SLIMElasticNetRecommender(URM_train)\n",
    "\n",
    "recommender.fit()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluated user 0 of 71568\n",
      "Evaluated user 10000 of 71568\n",
      "Evaluated user 20000 of 71568\n",
      "Evaluated user 30000 of 71568\n",
      "Evaluated user 40000 of 71568\n",
      "Evaluated user 50000 of 71568\n",
      "Evaluated user 60000 of 71568\n",
      "Evaluated user 70000 of 71568\n",
      "Recommender performance is: Precision = 0.3303, Recall = 0.1133, MAP = 0.2641\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'precision': 0.3302657049858151,\n",
       " 'recall': 0.11332678222135233,\n",
       " 'MAP': 0.2640531561329721}"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate_algorithm(URM_test, recommender)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As opposed to SLIM BPR, ElasticNet has some advantages, namely it relies upon a parallel solver and works on all OS, whyle Cython is sequetial and needs to be compiled or adapted. It has however disadvantages like a VERY slow training on some datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}