{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "import dateutil\n",
    "import gzip\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import scipy\n",
    "import sklearn\n",
    "import tensorflow as tf\n",
    "import time\n",
    "from collections import defaultdict\n",
    "from fastFM import als\n",
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is available at http://cseweb.ucsd.edu/~jmcauley/pml/data/. Download and save to your own directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataDir = \"/home/jmcauley/pml_data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Autoregression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example based on Bay-Area bike-share data. Extract the time information from the events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "events = {}\n",
    "\n",
    "for f in [dataDir + \"201408_trip_data.csv\", dataDir + \"201402_trip_data.csv\"]:\n",
    "    f = open(f, 'r')\n",
    "    l = f.readline()\n",
    "    for l in f:\n",
    "        l = l.split(',')\n",
    "        tripID = l[0]\n",
    "        timeString = l[2]\n",
    "        timeUnix =\\\n",
    "            time.mktime(datetime.datetime.strptime(timeString, \"%m/%d/%Y %H:%M\").timetuple())\n",
    "        events[tripID] = [timeUnix, timeString]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the earliest event (so that we can sort events from the first to the last hour)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "earliest = None\n",
    "for event in events:\n",
    "    if earliest == None or events[event][0] < earliest[0]:\n",
    "        earliest = events[event]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "earliestTime = earliest[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count events by hour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hourly = defaultdict(int)\n",
    "\n",
    "for event in events:\n",
    "    t = events[event][0]\n",
    "    hour = int(t - earliestTime) // (60*60)\n",
    "    hourly[hour] += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Autoregressive feature representation. Here we don't include a bias term, though could easily include one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature(hour):\n",
    "    previousHours = []\n",
    "    # Features for last 5 hours, one day ago, one week ago, and one year ago\n",
    "    for i in [1,2,3,4,5,24,24*7,24*7*365]:\n",
    "        previousHour = hour - i\n",
    "        previousHourExists = previousHour in hourly\n",
    "        if previousHourExists:\n",
    "            previousHours += [hourly[previousHour]]\n",
    "        else:\n",
    "            previousHours += [0]\n",
    "    return previousHours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feature(x) for x in hourly]\n",
    "y = [hourly[x] for x in hourly]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observation one week ago is the most predictive, followed by the observation from the previous hour:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.47105067, -0.28432565,  0.10555845,  0.01445134, -0.02131945,\n",
       "        0.17501989,  0.53969635,  0.        ])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sliding window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse ratings and timestamps from (a small fraction of) Goodreads fantasy novel data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratingsTime = []\n",
    "z = gzip.open(dataDir + \"goodreads_reviews_fantasy_paranormal.json.gz\")\n",
    "\n",
    "for l in z:\n",
    "    d = eval(l)\n",
    "    t = dateutil.parser.parse(d['date_updated'])\n",
    "    ratingsTime.append((t,d['rating']))\n",
    "    if len(ratingsTime) >= 50000:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sort observations by time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50000"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratingsTime.sort()\n",
    "len(ratingsTime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep track of a window (wSize) of ratings and timestamps (the raw time is just for printing the plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "wSize = 1000\n",
    "x = [r[0] for r in ratingsTime] # as raw times\n",
    "y = [r[1] for r in ratingsTime] # ratings\n",
    "xu = [time.mktime(d.timetuple()) for d in x] # as unix times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use a dynamic-programming approach to build the sliding window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "xSum = sum(xu[:wSize])\n",
    "ySum = sum(y[:wSize])\n",
    "sliding = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(wSize,len(x)-1):\n",
    "    xSum += xu[i] - xu[i-wSize]\n",
    "    ySum += y[i] - y[i-wSize]\n",
    "    sliding.append((xSum*1.0/wSize,ySum*1.0/wSize))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X and Y coordinates for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [a[0] for a in sliding]\n",
    "Y = [a[1] for a in sliding]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEZCAYAAABB4IgrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4FFXW+PHvSYgBgSBLwhIg7CA7CIiCEERWRRRHR3QARUVxG8RxNl8FfqOOjo7jq/iqoI7IsLiDOiCLEtmEECDsOyRA2An7EiC5vz+qOnY63Ukn6U53us/nefJYXXWr6nQ31um71C0xxqCUUkr5U0SgA1BKKRX6NNkopZTyO002Siml/E6TjVJKKb/TZKOUUsrvNNkopZTyO002YURE7hORHwIdR6CISD0ROS0iUkrne1RE3vSi3CIRGWkvF/gdFaVsaRCRMyLSoJj75r4XfxORlSJybWmcS7mnySaIiUiaiJy3L5AHROTfInK1l/smiEiOiOR+x8aY6caY/v6LOLiIyB4Rudnx2hizzxgTY0rh5jIRiQKeB/5RlP2K8h0Fw/dpjKlsjEkLZAxeeh34W6CDCGeabIKbAW41xsQA7YEOwF+83Ffs/UvlV3yglVZtpQgGA1uMMYcCHYgC4Dugl4jEBTqQcKXJJvgJgDHmCDAPK+lYG0QGisgaETklIukiMs5pv5/t/560a0bXi8gIEVnitH+O3dSzXUQyRWSi07YIEfmniBwVkV0i8oRzTUlEHrDXn7b/O9Rt8CJXichbIpIhIvtF5F/2r35EZLOIDHQqGykiR0Skvf26q4gsE5ETIrJWRHo6lV0kIi+JyFIROQc0dDnvp0B94Ds7xj+41vbsY/zNPscZEZktItVE5D/2Z7pSROo7HbOFiMwXkeMiskVE7i7gexvg9B0gItEiMlVEjtnvZ6WIxLr5vFy/oz72uU6IyDs4/Xjw5ffpEsMDIvKt0+sdIvKZ0+u9ItLW6ZyN7OV/i8hEEfne/sx/EZGGTvsV9F5ERP5HrNr8IRH5REQq29s+EZFn7OU69jlH268bi8hxe7m6iHxnH/+4iOR+/saYLGA10M/jN6b8yxijf0H6B+wBbraX6wLrgTedtvcAWtnLrYGDwO326wQgGxCn8iOAxU6vc4BvgcpAPeAI0Nfe9hiwEagNVAEW2MeLAK4GTgFN7LI1gWs9vIf/BywHqtt/y4AJ9rYXgP84lb0V2GQvxwPHgH7269726+r260VAGtDCjinSw+fXy+m14zOJcDrGdqCB/RlsArYCvexjTgE+ssteDewFhmNdJNvZn1cLD+87GbjL6fUoYDYQbe/fAajkFMdI1+8IqAGcBu4EIoExwGV3ZUvyfbqJvSGQaS/Xtj/nvfbrRsBxp7LZQCN7+d/AUeA6+/P7DzDd3la9kPcy0v4uEuzP+itgir3tQWC2vTwU2AHMcNr2jb38CvB/jn8PQDeX9/W/wBuB/v86XP+0ZhP8ZonIaawL3WFgvGODMWaxMWaTvbwRmAn0dNm/sOalvxtjzhhj9mFd9Bw1p7uB/zXGHDTGnAJeddkvG2gjIuWNMYeNMVs8HP8+rORy3BhzHJiAdcEGmA7cLiLl7ddDgRn28v3Af40x8+z39yOQAuTWhIBPjDFbjTE5xphsD+cv7P3/2xiTZow5A8wFdhljFhljcoAvsJICwG3AHmPMp8ayDvga63Ny5xrgjNPry1gX3Gb2/muNMWcLiW0AsNEY840xJtsY8xZQWLNccb/PXMaYPcAZu4bZA6tGfUBEmtmvlzgVd/18vzHGrLY/v2lO5x9YyHu5D+uHVLox5jxWc/FQu+b1M9DdLtcDqx+sm/26J7/WIC9jJceG9jmWucR2But7UQGgySb4DTZWn01PrF/xNRwbRKSLiPxkNz2dBB513u6lw07L54FK9nIdYJ/Tttxl+2LwW2A0cNBuumju4fh1sBKlQzrWBQFjzC5gMzBIRCoAt2NdoMD6hXuP3RyUKSInsC4wtdzFVALO7/+Cm9eOzyMB6OoSz30u8Tg7gVXDcPgU66I9025OfE1EIguJzfU7wM1rV0X+Pj34GauG1wNIsv8SyXtxd8c5gRR0ftcY6mD923BIB8oBNY0xu4FzItIBuAn4nl+Tn3M8/wB2AfNFZKeI/MnlfJWBkwXErvxIk03wc/TZLMFq1vmn07bpwCwg3hhzDfABv/7SLOmIq4NYTXcO9Z03GmMWGGP6Yl1stwGTPRwnA+tC7ZAAHHB6PRProj0Yqwltj71+H/CpMaaa/VfVWCOfXncOo5D34MtRZ/uAJJd4YowxT3govx5olhuI9Uv7b8aYVsCNWDWl4R72dTiIy+eO1TxWHAV+n24sxkou3bEu5ouxLuw9KDjZFHT+gt7LAfL/O7nMr8nzZ+A3QJQx5qAdzwismkoqgDHmnDHmD8aYxlg/XMaKSC+nY14LrCtG7MoHNNmULW8BfUSkjf26EnDCGHNZRLpgXbQdjmK14Tcu5rk+B35vd8heA/zRsUFE4kTkdrGGYV8GzmI1q7kzE/gfEakhIjWw+mmmumzvi1VLmu60/j9YNZ6+dud2eRHpKSJ1ivAeDmH1MTgr7qi174FmIvI7ESknIlEi0klEWngoPwfrYm2dVCRRRFrbzUJnsT43T5+Zw3+BliJyh1iDJ36P55pUYTx+nx44ajYVjDEHsJrO+mM1Ba4txvkLey8zgGdEpIGIVAJeBmbazXFgJZcn7f+CVdN6ElhqjNUhIyK3iojj3/sZ4ArW/wOISDRWX9KCYsSufECTTXDL88vcGHMMq3bzor3qCeBvInIK+B/gM6eyF7D+h11mN/t0Kez4Lq8nA/OxfqGvxrpYXLH/548AxmLVWo5h/dod7eE9vITV17Ie61dlih2XI85DwC9AV5f492PVdv6KlTjTgT/w679Zb2otrwIv2O9/rJv9vK752P0rfYF7sX6FH7CPf5WHXb4DmouI44JaC/gSa2DFJqz+lP8UFIfdx3U38BrW59wYWFpQmAW8Luj7dHfuHVgX7MX26zNYTVS5F/eCYi/Ge/kY60fIYvs854Gnnbb/jPXjylGrWgpUIG8tqymwUETOYA1EedcY49h+O7DI6FD0gJG8/278dBLr11wKsN8Yc7vLtquw2rOvw/pH+FtjzN78R1GBJCL9gfeMMQ0LLawAEJGHgZbGmLGFFi5l4fZ9isgvwEPGmM2BjiVclVbN5vdYHcHuPIQ1zLIpVjNRke64Vv5hN1sNsJs84oFxWKOvlJeMMR8GS6IJ9+/TGHODJprA8nuyEZG6WMMeP/RQZDBW0xBYzQy9/R2T8opgDVPOxGp22YR1gVJlk36fKqDKlcI5/gU8h3UjmTvx2EMgjTHZInJSRKoZYzJLITblgd3n466fR5VB+n2qQPNrzUZEbgUOG2NSsX5ZeTMSKNjmuFJKKVVC/q7ZdMO6Q3wg1siRyiLyqTHG+f6C/Vjj7Q/YN7nFuKvViIj/RzIopVQIMsYE/Ed8qYxGAxBrEsVn3YxGexxobYx5XETuBe4wxtzrZn9T1FjHjx/P+PHjSxC1UkqVTY7rn4gERbIJyH02IjJBRG6zX34E1BCRHViT8/05EDEppZTyn9IYIACAfXPVz/byOKf1WcA9pRWHUkqp0ldqzWgl5dqM1qBBA9LT0wvYQwWThIQE0tLSAh2GUmEjKSmJxMTEoGlGK7PJxv4AAxiRKgr9vpQKjGBJNjo3mlJKKb/TZFOI9PR0atWqxc0330z//v1z17/xxhvcdNNNDBs2jOzs7CKtK8q5Fy1aBMDhw4f5+9//7qN39auxY8d6VeN48MEH2btXp6xTShWPJhsv9O3bl59++okffvgBgKNHj5KUlMSSJUto06YNs2bNKnRd27ZtmTVrltvje7rYp6Wl8dNPPwFQs2ZN/vKXv/j8vb355puIBLyGrZQKcZpsvPDTTz/Rs2dP3nrrLQBSUlJITEwE4JZbbmH58uWFruvduze//PJLnuM++OCDPPXUUwwYMICDBw9y880306NHD5588kkAJk2axNSpU+nTpw/p6ekMGzYMgBtuuIFRo0bRsWNH5s+fD8D777/PjTfeyJ///Gd69eqV5zx//OMf2bRpEwsXLqRDB+spxw888ABHjx6lV69e5OTkMGHCBEaMGEGfPn0YNWoUYCW7rl27MnjwYPbssZ5pdvr0aQYNGkRiYiJjxowBYOjQoZw5c4YPP/yQIUOGADBo0CDto1FK5dJkU4g6deqwY8cOFi1axI8//sjGjRs5efIkMTExAFSpUoWTJ09y6tQpr9a56t69Oz/88AOxsbEsXLiQxYsXc+rUKXbt2sWoUaMYPnw4CxZYz3ty1EAyMzP5+9//zvfff8/7779PdnY2U6ZMYfny5QwZMiRfTaVbt24sW7aMZcuWER8fz9mzZzly5AixsbF5ynbs2JEFCxawd+9eTp8+zeuvv85bb73F119/zfHjxwErAd57770kJSVx/vx5kpOT6dq1K7/88gurVq0iOjqaK1euEBkZGfY1pqysLJ566im337tS4Sbkko2IFOmvMFFRUVSoUIGIiAhuvfVWNm7cSJUqVTh16hRg/dKvWrWqV+uuueaafMe/7rrrADh27Bh33XUXvXr1YtmyZRw4cCBfWYe4uDiqV69OnTp1OHXqFMePH6d+feuJu+3bt89X/sYbb2T58uXs2bOH+++/n9mzZ1OzZs185Vq3bg1YCfbkyZPs3r2b9u3bExkZSdu2bQHYtWsXHTt2zI19165ddOvWjeXLl3Px4kXatm3LZ599lluDCmfPPPMMEydOZMmSJYEORamAC7lkY4wp0l9hzp49m7u8bNkyGjduTOfOnVm82Ho67cKFC+natSudOnXyap2riAjrK5g+fTp33nknixYt4sYbb8QYQ1RUFFeuXCn0/VavXp39+/djjGHduvyPWI+NjeXgwYNERkZy44038sYbb3DjjTfm7u/gSL6OdY0aNSI1NZXs7GzWr18PQJMmTUhJSQGs5sTGjRvTvn175s+fT82aNfMdP1xNnTqVhQsXMmbMmHzNp0qFo5BLNr62ZMkSOnXqRPfu3albty6dO3cmNjaWm266iZtuuol169Zxxx13eL3OmXPN6uabb+aNN95gyJAhnD9/HrBqGsuWLWPo0KEe9wOIjIxk+PDhdOvWjZkzZxIVFZXvfdSpU4d27dqRkJDAsWPH6NatW55jOR/Tsfzss88yZswY7rrrLmrVsp5u/PDDDzNz5kx69uxJ+fLl6dKlC+XKlSM6Opru3bvTpUsXtm3b5jaxhov169czduxYvvrqK/r166fJRin0ps6QkZ2dTWRkJMnJyfz73//mvffeC3RIeYTL93Xy5Ek6d+7M+PHjuf/++zlx4gT169fnxIkTlCtXarNDKZVLb+pUPvXOO++QmJjIM888w7PPPhvocMKSMYYHH3yQvn37cv/99wNQtWpV6taty8aNGwMcXWBdunSJ1157LdBhqADSmo0qFeHwfa1atYqhQ4eyadMmoqOjc9ePHDmSzp07M3r06ABGF1irVq3ihhtu4NKlS7n9lKp0aM1GqRDz9ddfc8899+RJNGDdFxXu/TYpKSlkZ2eTmalPew9XmmyU8gFjDF999VXuTa3ONNmQO4LxyJEjAY5EBUqZ7bFMSEgI+5sGy5KEhIRAh+BXmzdv5uLFi7n3TTm79tprOXLkCEePHiU2NjYA0QVeSkoKMTExHDlyhJYtWwY6HBUAZbZmk5aWVuR7avQvcH+h/iybr7/+2u3sDWANTe/SpQsrV64MQGSBd/78eXbs2EHPnj21ZuMju3bt4pFHHgl0GEVSZpONUsHEkWw8CeemtHXr1tGyZUvq1aunycZH1qxZw8cff8yxY8cCHYrXNNkoVUK7d+/mwIEDuTfKuhPOySYlJYVOnToRFxenycZH0tLSyMnJYc6cOYEOxWuabJQqoW+++YbBgwcTGRnpscz111/PqlWrCp1+KBRpsvG9tLQ02rdvz3fffRfoULymyUapEiqsCQ2gWrVqxMfHh+XNnatWrcpNNocPHw50OCEhPT2dJ554ggULFpCVlRXocLyiyUapEjh48CCbN2/m5ptvLrTsDTfcwIoVK0ohquBx5swZ0tPTadWqldZsfCgtLY0uXbrQsmVLfv7550CH4xVNNkqVwKxZs7j11lu56qqrCi0bjv02a9eupW3btkRFRVGzZk1NNj7gGN2ZkJDAoEGDykxTml+TjYhEi8hKEVkrIhtEZJybMvVFZKGIrBORn0Skjj9jUsqXvGlCcwjHZOPorwG0ZuMlxzOwPDl+/DhRUVFUqVIlN9mUhamg/JpsjDFZQC9jTAegPTBARLq4FHsD+MQY0w74f8Cr/oxJKV/JzMwkOTmZfv36eVW+ZcuWHD58OCiHq/orCTgnmypVqnDhwgUuXrzol3OFghUrVtCoUaMCk0daWhoNGjQAoFWrVkRERLBhw4ZSirD4/N6MZow5by9GY81Y4PoptgQW2WWTgMH+jkkpX5g9eza9e/emYsWKXpV33NwZbP02O3bsoF69emRkZBRr/+3bt3vc5pxsRIS4uDiOHj1arPOA1YRU0PnKunfeeYfMzEz27t3rsYxzshGRMtOU5vdkIyIRIrIWOAQsMMascimSCgyxyw4BKolIVX/HpVRJzZgxI9+D7QrTvXt3FixY4KeIimfp0qVcvnyZ//u//yvyvmlpaTRv3pyJEyfm23by5EkOHjxIixYtcteVtClt6dKlJCYmFnm/tWvX8v333xf7vKXh0KFDzJkzh86dO7Np0yaP5ZyTDVBmko3f50YzxuQAHUQkBpglIi2NMZudijwHTBSRB4DFQAaQ7e5Y48ePz11OTEws1j86pXzhyJEjJCcnM2vWrCLtN3LkSNq1a8f48eOpWjU4flMtXbqUp59+mkmTJvH8889z9dVXe73v1KlTGTx4MC+//DJNmzbN06S4Zs0a2rdvn+f+o5Imm3nz5nHw4EEuXrxI+fLlvdrHGMMTTzxBbGwst912W7HP7W+TJk3it7/9LVdffTUbN25k4MCBbsulpaXRtGnT3Nc9evRg27ZtHDp0iFq1apGUlERSUlIpRV0EpTk/FvACMLaA7RWBvR62GaWCxcSJE839999frH2HDx9uXn75ZR9HVHzNmzc3a9euNYMGDTLvv/++1/vl5OSYJk2amJUrV5olS5aY2NhYs2nTptztr732mhkzZkyefYYNG2Y++eSTYsfauXNnIyJm+/btXu/zww8/mJiYGNOuXbtin9ffLl26ZOrUqWPWr19vPvzwQzN8+HCPZW+77TYza9asPOvuuece8+GHH7otb187Az4/or9Ho9UQkSr2cgWgD7DVpUx1+XX2wr8AH/szJqV8YebMmdx7773F2ve5557j7bffDoqO8mPHjnHw4EHatGnDM888w1tvvUVOTo5X+y5fvpxy5crRuXNnunfvzhtvvMGgQYNy+2Sc+2scSnJj57Fjx9i6dSs33HAD6enpXu1jjGHcuHG88sorXu8TCF9//TVNmzalTZs2tG7dusCbf12b0aBsNKX5u8+mNrBIRFKBlcA8Y8wcEZkgIo76bCKwTUS2AnHAy36OSSkA/vnPf7Js2bIi77d37162bNlC3759i3Xe1q1b06lTJ6ZMmVKs/X1p+fLlXH/99URGRpKYmEh0dDTz5s3zat9PP/2U4cOH5850PXz4cO655x6GDBlCVlaW22RTknttfvzxR3r27EnTpk29Thw//PADZ8+e5bHHHiMrK4vTp08X69z+NnHiRJ566inAGrW4detWsrPz9yYYp3tsnA0cOJCffvqJCxculEq8xeHvoc8bjDEdjTHtjTFtjTEv2+vHGWO+t5e/MsY0M8a0MMaMMsZc9mdMSjksWLCgWPe9fP7559x5551e3cjpyZ/+9CfeeOMNtxeU0rRs2TK6d+8OWCObHLWbwly4cIEvvviCYcOG5Vn/8ssvExsby3333cfx48fz9C1Ayfps5s2bR9++fUlISPAq2RhjGD9+POPGjSMyMpKEhIQCR3kFSmpqKmlpaQwebA3ErVy5MrGxsezZsydf2czMTMqVK8c111yTZ321atXo0KEDP/30U6nEXBw6g4AKW0ePHi3Wxac4o9Bcde/enRo1avDNN9+U6DgltXTp0jyzVd97772sX7++wNFQAN999x3XXXcddevWzbM+IiKCqVOnsmfPHjp27EhERN5LTHGTjTGG+fPn069fP6+Tzdy5czl//jx33XUXQNAmm4kTJ/LYY49Rrtyv47VatWrl9jtw14TmEOxNaZpsVNgqTrLZvn07Bw4coGfPniU6t4jwpz/9iddee83vd3///e9/d9u8cvHiRVJTU7n++utz10VHRzN69OhCazdTpkxh+PDhbrdVrFiRH374gXfeeSfftuImm82bNxMVFUXTpk29SjaOvppx48blJrz69esHXb9NZmYmX331Vb4HoXnqtyko2dx+++3MmzcvaGcT0GSjwpIxpljJZubMmdxzzz0FPk7AW7fffjtnzpzx6zDVM2fO8Ne//pUvvvgi37bVq1fTokULKlWqlGf9Y489xpdffunx5stDhw6xfPnyAqfpiYuLo3Xr1m7XFyfZzJ8/n759+yIiXiWbOXPmcOnSpTwxBlvN5uLFizz66KPceeedxMXF5dlWUM3G0yPWmzVrxrp169w+LTYYaLJRYencuXNkZWUV6eJjjGHGjBnFHoXmKiIigueee4433njDJ8dzZ+vWrURFRTFp0qR825YtW+b2gW9xcXHcddddvPvuu26POX36dO644w6vZ05wFhsby5EjR4r863vevHm59/DUq1ePAwcOFNjf9e677/LnP/85TzNeMNVsjh49ys0330xkZKTbm2lbtWpV5JoNQExMjC/D9ClNNiosHT16lHr16nHu3DnOnTvn1T7r16/nwoULdO3a1Wdx3HfffSQnJ7Nz506fHdPZli1buPPOO9mzZ0++X8rOgwNc/eUvf+GDDz5g+vTp+bYV1IRWmPLly1OhQoVCJ5t0dvHiRZYtW5b7GIfo6GiqV6/OgQMHPO6zYcMGbrjhhjzrgqVms2XLFq6//np69+7N9OnT3d6ceu2117Jjxw4uX847Xio9Pb3AZBPM/D6DgFLB6OjRo8TFxREdHc2+ffvyTKniiePeGl82U1SoUIGRI0fy3nvv8c9//tNnx3XYvHkzbdq0oXnz5nzwwQe8/fbbgFVLW7Zsmccpaho3bszChQvp27cvly9fZsSIEQCsW7eOU6dOlajPytGU5jqiypMlS5bQtm3bPOUdTWn16tXLV/7MmTOcOHGC+vXr51kfiJrNxo0bOXz4MCdOnODEiRMcPnyYd955h3/84x+5n6k7V199NfHx8ezcuZNrr702d31hNZtgpjUbFZaOHDlCbGws9evX9+rXrjGmRDdyFmT06NFMmTKl0BrWRx99xOOPP16kY2/evJmWLVvy8MMPM23aNM6ft+bF3bZtG5UqVSI+Pt7jvq1ateLHH3/k+eefZ/LkyYBVqxk2bFi+UWZFUbNmzSLd2Onor3FWUL/N1q1badasWb4Y4+PjOXz4cL7agr/MmTOHxMREXn75ZWbMmMHKlSs5deoU33zzTYGJxqF169Z5aqOOe2zKarLRmo0KS0ePHiU2NpbIyEivks26deuIioqiXbt2Po+lQYMG3HjjjUyfPj3fqCSArKwsnnrqKZYsWcK+fft48803vZ4XzJFs6tevT9euXfn888954IEHPPbXuGrRogWLFi3illtu4cKFC0yfPp3FixcX+T06K+oggXnz5uXrcyos2birqToe4JaRkVEqF+zZs2fz17/+lbFjxxZrf0e/zW9+8xsATpw4QUREhNc1wmCjNRsVlhzJxtt2/LS0NFq0aOG3kT5PPvkkEydOzNdxvm/fPm666abcZ+e0atXK60cUXLhwgYyMDBo3bgzAo48+mnvR9jbZADRt2pSkpCTefPNNGjVqRLNmzYrwzvIrSrI5ePAg+/fvzzcTQUHJZsuWLXmanlz381W/zdatWz1O7WOM4fvvvy/RxJ+uNZuyXKsBTTYqTDmSjbft+IcPH6ZmzZp+i+eWW27h4sWLLF26FICcnBymTZtGly5duPvuu/niiy+oXLkyiYmJXg+V3rZtG40bNyYqKgqwpjTZt28fGzZsYOnSpR4HB7jTsGFDfvnlFz755JOivrV8ipJsfvzxR3r16pXnhkcoPNl46oPzZb9N//79PU7tk5qaSsWKFUuUmF1HpGmyUaoMck423vzS9XeyiYiI4IknnmDixInMnTuXjh078s477/Dll1/y3HPP5daoipJsHE1oDuXKleOhhx7ipZde4siRI7Rq1apIMdauXbvEtRooWrJZtWqV29F/hTWjFVSz8UWyuXLlCvv37/f4iImS1moAmjdvTlpaGllZWYAmG6XKpGBLNgAjRoxgzpw5jB07lnHjxvHLL7/ka+rq1q0bKSkpXs0Y7ZpsAB566CG+/PJLunbt6pMbU4ujKMlm9erVXHfddfnWO5rDXJsdL1++zJ49e/LNyebg7fddmAMHDhAVFcW3337rtinNF8kmOjqaBg0asG3bNkCTjVJlkiPZ1K1bl/379xc6rf6hQ4f8nmyqVKnC+vXr2bBhA3feeafb/qGYmBhatWrFypUrCz2eu2RTr149brvtthJPt1MS3iab7OxsUlNT6dixY75tlStXpnz58hw7dizP+p07d1KvXj2io6PdHtNXNZv09HQ6dOhAtWrVSE5OzrPt8OHDbNu2rUjNlJ4499uU5XtsQJONClOOZFOhQgWuueaaQofiHj58mFq1avk9roYNG+brn3DlbVOau2QD1gwAf/jDH4obYol5m2y2bt1KnTp1PI6+cpc4CmpCA9/VbPbu3UtCQgJ33HFHvqa0uXPn0qdPnxLNCu7g3G+jNRulyiBHsgHvLkCl0YzmLW+SzaVLl0hPT3fbnFSxYsXcQQOB4G2ySUlJcduE5uAu2RQ0Eg1+HSBQ0skq09PTPSYbXzShOThqNp6eY1OWaLJRYefChQtcvnyZypUrA2Uv2XTr1o1Vq1YV2G+zY8cO6tev77E5KZDtBiHKAAAgAElEQVSqVavG6dOnC725cvXq1fmGPDvzlGwKmg0iJiaG6Ohojh8/XrSgXaSnp1O/fn2uu+46zp49y9at1gOIL126xMKFCxkwYECJju/gqNmcPHkSoMzeYwOabFQYctRqHH0ihSWb8+fPc+nSJapUqVJaIRbIm34bT01owSAiIoIaNWp4nFXawZuaTVpaWp51hTWjOfYrab+NoxktIiKCwYMH59ZuFi9eTIsWLfLN4lxcTZo0ISMjg82bN9OgQYOgndHZG5psVNhxbkKDwpONo1YTTP+jF9aUFszJBgpvSrty5Qrr1693OzjAwTVpGGM8zh7gzBf9No6aDZCnKc2XTWhA7jN85s6dW6b7a0CTjQpDxU02wSTUk82WLVuIj48vcMp812Szf/9+KlWqVGhTU0lrNsaY3JoNQM+ePdm+fTsHDx70ebIBq9/mv//9ryYbpcqaUEg2hfXblPVkU1h/DeRPGt40oUHJazaZmZlERUXlJsKrrrqKAQMG8Prrr3Px4kWfz5/XqlUrUlNTNdkoVdaEQrJx9Nu43uMBVhPUzp07ad68eQAi805hyaaw/hqA6tWrc+nSJU6fPg0UPhLNoaQ1G+cmNIfBgwfz9ttvc9ttt/m8udXxxFNNNkqVMa7JJi4ujrNnz+ZOv+8qGJMNeG5K27VrF3Xq1OHqq68u/aC85Iuajesjor3pr4GS12ycm9Ac+vfvT2RkJLfeemuxj+uJY1ohTTZKlTGuyUZEqFevHvv27XNb/tChQ6VyQ2dR9ezZ022yCfYmNCg42TgGB3To0KHQ4zgnm0DWbGJiYpg/fz79+/cv9nE9adiwIbVq1aJhw4Y+P3Zp8muyEZFoEVkpImtFZIOIjHNTpp6I/CQia0QkVUR8M0BdKQ9ckw0U/Gs3WGs23bt3Jzk5OV+/TVlINjVr1vSYbDZv3kz9+vVz74MqiGuy8aZmU7NmTU6dOsWFCxeKFrTNcUOnq549e/rlZlnHM5fK8j024OdkY4zJAnoZYzoA7YEBItLFpdj/AJ8ZYzoCQwH3z6lVykdCJdnExMTQsmXLfP02mzdv9uoXfiDFxcV5nCLIm/4aB0eyOXHiBOfPny/wyaMOERER1KtXr9hNae6a0fwtkDM++Irfn9RpjHE0hEfb53OdJyIHcIxvvAbI8HdMKryFSrIBq9/mq6++4sKFC+zYsYOdO3eyYMECnnnmmUCHVqCCmtFSUlIK7a9xSEhIIDU1Nbe/xtvOeccNocUZROGuGU0Vzu99NiISISJrgUPAAmPMKpciE4BhIrIP+B54yt8xqfDmLtkU9Eu3tCbhLI5Bgwbx/fff8/rrr7Np0ybq1avH1KlTva4ZBEpsbCxHjhxxO0eZp8cKuOOo2XjbhObQvHlztmzZ4nV5Z4Go2YSC0qjZ5AAdRCQGmCUiLY0xm52KDAX+bYz5l4h0Bf4DuH2q0/jx43OXExMTSUxM9FvcKjRdunSJc+fO5Wv/rlu3LhkZ+SvVFy5cICsrK2imqnF10003sWvXrkCHUWQVK1YkMjKSs2fP5umbuXz5Mhs3bvRqcAD8mmy8vcfGoU2bNqxevbrIcV+4cIFTp04FbU0XICkpyesH7JUmvycbB2PMaRFZBPQHnJPNQ0A/u8wKESkvIjWMMcdcj+GcbJQqjsOHD1OjRg0iIvJW6uPj49m/f7/b8nFxcUE1VU2oqFu3LrNmzWLYsGG56zZt2kSDBg2oVKmSV8eoXbs2mZmZpKam8vjjj3t97jZt2hTrEdd79+6lbt26+f79BBPXH+ITJkwIXDBO/D0arYaIVLGXKwB9gK0uxdKBW+wy1wLR7hKNUr6Qmprq9g5vTzWbYO6vKeumT5/Oiy++yB//+EeuXLkCFG1wAFgjteLj41myZEmRmtEcU/cX9tA8V9qEVnz+Ts+1gUUikgqsBOYZY+aIyAQRcUwg9AfgEbvMNGCEn2NSYSw5OZkuXVwHRFpPyczOzs69G90hmPtryrqOHTuSkpJCamoqffv25ciRI17dzOkqISGBK1eu0LhxY6/3ueaaa6hWrRp79uwp0rl0cEDx+bUZzRizAcg3basxZpzT8hag5M9PVcoLycnJPPVU/jEoIkJ8fDwZGRl5Jn8sjcdBh7Pq1aszd+5cXnzxRTp16kRkZCS/+93vinSMhIQEDh8+XOThwW3atGHDhg1FSlJasym+4G14VMrHjDGsWrWKzp07u93urilNm9H8LzIykpdffpm3336bChUq0L59+yLtn5CQUKQmNIfWrVuzYcOGIu3j6YZOVThNNips7Ny5k5iYGI/Jw90gAU02peeOO+5g8+bNVKxYsUj7DRkyhNGjRxf5fI6aTVFoM1rxabJRYcNTf41D3bp1NdmUQe3ataNPnz5F3q9NmzZs3LixSPtoM1rxabJRYcObZOOuGU0HCISmFi1asGfPHrKysvKsz8nJYebMmfluOM3OziYjI4O6deuWZpghQ5ONChuFJRt3zWg6QCB0RUdH06hRo3wzCaSkpDB06FB++eWXPOsPHTpEtWrVKF++fGmGGTI02aiwcOnSpUKfaa/NaOHHXb/N7NmzqVWrFpMmTcqzXgcHlIwmGxUWHENcC7oz3bUZ7eLFi1y8eLHMT+2uPHOXbGbNmsXkyZOZPXs2J06cyF2vgwNKRpONCguFNaGBNTnkqVOncp8P46jV6FQ1oct1+POOHTvIzMxk4MCBDBgwgKlTp+Zu08EBJaPJRoUFb5JNREQEtWvX5sCBA4A2oYUD1xFps2fP5vbbbyciIoJRo0bxwQcf5A4U0JpNyZSpZONuOnKlvOFNsoG8TWk6OCD0NWjQgJMnT+Y2l82aNYs77rgDsJ68mZ2dzbJlywDtsympMpVsjh3T+TlV0Z0+fZr09HRatXL75Io8nEekac0m9EVERNCqVSs2btzIkSNH2LhxIzfffDNgTWE0atSo3IEC2oxWMmUq2Wzbti3QIagyaPXq1bRv396rubOcR6RpsgkPjkEC3333HX379iU6Ojp324gRI/juu+/IzMzUZrQSKlPJZvv27YEOQZVB3jahQd5mNL2hMzw4ks3s2bNzm9Acqlevzq233srbb78NoCMTS6BMJRut2ajiKEqy0Wa08NOmTRtWrFhBUlISAwcOzLf90Ucf5c0336R+/fo6MrEEylSy0ZqNKo7k5GSPMz27cm5G0wEC4aF169akpqZy/fXXu625dO/enbp162p/TQmV2mOhfUFrNqqoDhw4wPnz52nUqJFX5V2b0TTZhL7Y2Fhq1qyZrwnNQUQYO3as28eGK++VqWSze/durly5QrlyZSpsFUCrVq2iS5cuXjd/1K5dm8OHD3PlyhXtswkjr7zyCrfffrvH7Q8//HApRhOaylQzWq1atUhPTw90GKoMKUp/DUBUVBTVq1dn7969XLhwQTuEw8TIkSOpUaNGoMMIaWUq2TRr1kyb0lSRFDXZgNWUtmbNGuLi4rRDWCkfKVPJpnnz5jpIQHktJyenwMdAexIfH09KSor21yjlQ2Uq2WjNRhXFjh07qFq1KnFxcUXaz1Gz0f4apXynTCWb5s2ba7JRXitOExpYNZvVq1drzUYpHypTyaZZs2bajKa8VtxkU7duXTIzMzXZKOVDfk02IhItIitFZK2IbBCRcW7KvGlvXyMi20Qk09Px6tevT2ZmJmfPnvVn2CpElCTZAJpslPIhv96wYozJEpFexpjzIhIJLBORucaYZKcyYx3LIvIk0N7T8SIiImjSpAk7duygQ4cO/gxdlXFZWVls2LChwMdAexIfHw9oslHKl/zejGaMOW8vRmMlt4IeSjMUmFHQ8XSQgPLG+vXradq0KRUrVizyvo5kowMElPIdvycbEYkQkbXAIWCBMWaVh3L1gQbATwUdT4c/K28UtwkNoGLFilStWlVrNkr5kN/nfTHG5AAdRCQGmCUiLY0xm90UvRf40hTwOM7x48ezadMmdu3aRY8ePUhMTPRT1KqsS05Opnv37sXe/7PPPqNZs2Y+jEip0pGUlERSUlKgw8hHSvNRyyLyAnDOGPOmm21rgMeNMSs87GuMMaxYsYInn3ySlJQUf4eryrCWLVsyY8YM2rVrF+hQlAooEcEYE/CpMPw9Gq2GiFSxlysAfYCtbsq1AK7xlGicOYY/l2aSVGXLqVOn2Lt3r1ePgVZKlQ6vmtFE5G03q08BKcaY2QXsWhuYIiIRWIntM2PMHBGZAKwyxnxvl/stMNObWKpVq0Z0dDSHDh2idu3a3uyiwszq1avp0KGDzg6uVBDx9v/G8kAL4Av79V3AHqCdPbR5jLudjDEbgHxjT40x41xeT/A6Yn4dJKDJRrlTksEBSin/8LYZrS3QyxjzjjHmHeAWrORzJ9DXX8F5osOfVUE02SgVfLxNNlWBSk6vKwLVjDHZQJbPoyqEDn9WBdFko1Tw8bYZ7R9AqogkAQL0AF4RkYrAQj/F5lGzZs1YunRpaZ9WlQEZGRlkZWXRoEGDQIeilHLiVbIxxnwkInMAx8/FvxpjDtjLz/klsgLo7M/Kk6I+BlopVTqKMvQ5AjgKnACaiEgP/4RUuMaNG7N3714uX74cqBBUkNImNKWCk7dDn1/DGp68CcixVxtgsZ/iKlB0dDTx8fHs3r2b5s2bByIEFaSSk5N59tlnAx2GUsqFt302dwDNjTGlPhjAE8cgAU02yiEnJ4eUlJQiPwZaKeV/3jaj7Qai/BlIUenwZ+Vq+/btVK9enRo1agQ6FKWUC29rNuexRqP9iNNQZ2PM036JygvNmzdn7dq1gTq9CkLaX6NU8PI22Xxr/wWNZs2aMXOmVzPcqDChyUap4OXt0Ocp/g6kqHT4s3KVnJzM0KFDAx2GUsqNAh8xICKfG2PuEZENuHnCpjGmrT+Dc4klz6NujDFUrlyZjIwMqlSpUlphqCCVlZVFtWrVOHr0KFdffXWgw1EqaATLIwYKq9n83v7vbf4OpKhEhKZNm7J9+3YdfaRYt24dzZo100SjVJAqcDSaMeagvfi4MSbd+Q943P/hFUznSFMO2l+jVHDzduhzHzfrBvgykOLQ4c/KQZONUsGtwGQjIqPt/prmIrLe6W8PsL50QvRMazbKQZONUsGtsD6b6cBc4O/An53WnzHGZPotKi9pzUYBnDx5koyMDK699tpAh6KU8qDAZGOMOYX1+OehACISh/XUzkoiUskYs9f/IXrmqNnk5OQQEVGUOUVVKElJSdHHQCsV5Ly6QovIIBHZgfUo6J+BNKwaT0DFxMQQExPDgQMHCi+sQpY2oSkV/LytDrwEdAW2G2MaAr2BFX6Lqgi0KU1pslEq+HmbbC4bY44DESISYYxZBHTyY1xe00EC4c0Yw8qVKzXZKBXkvG3kPikilbCeXzNNRI4A5/wXlve0ZhPeMjIyyM7OJiEhIdChKKUK4G3NZjDWzM/PAD8Au4BB/gqqKHSOtPDmaELTx0ArFdy8SjbGmHPGmBxjzBV7Us6JQP/C9hORaBFZKSJrRWSDiIzzUO4eEdlkl/lPUd6ANqOFN+2vUapsKOymzhgR+YuITBSRvmJ5EuthavcUdnD7yZ69jDEdgPbAABHJc2UQkSbAn4AbjDFtgDFFeQMNGzYkIyODrKygeYioKkWabJQqGwqr2UwFmgMbgIeBRcDdwB3GmMHenMAYc95ejMbqI3KdPfoR4F1jzGm7/DHvQrdERUWRkJDArl27irKbCgHZ2dn6GGilyojCBgg0smsbiMiHwEGgvjHmorcnEJEIYDXQGCuprHIp0swutxQr+U0wxszz9vjw6yCBli1bFmU3VcZt27aNuLg4qlevHuhQlFKFKCzZXHYsGGOyRWR/URKNvV8O0EFEYoBZItLSGLPZJYYmQA+gPrBYRFo7ajrOxo8fn7ucmJhIYmIioP024Uqb0JTKLykpiaSkpECHkU9hyaadiDgu+gJUsF8LYIwxMd6eyBhzWkQWYQ0scE42+4EVdlJKE5HtQFOs2lAezsnGWbNmzVixIijuMVWlaNWqVZpslHLh/EMcYMKECYELxklhz7OJNMbE2H+VjTHlnJYLTTQiUkNEqtjLFbAeVbDVpdgsoJejPFai2V2UN6HDn8OT1myUKjv8PXNhbWCK3W8TAXxmjJkjIhOAVcaY740x8+yRbpuAK8AfjDEninISbUYLPxcvXmTz5s106NAh0KEopbwgxrgODgtOImI8xWqMoUqVKqSlpVGtWrVSjkwFwsqVKxk9ejRr1qwJdChKBTURwRgT8LueQ2JefhHR2k2Y0SY0pcqWkEg2oHOkhRtNNkqVLSGTbLRmE1402ShVtoRMstGaTX7GGC5duhToMHzuxIkTHDx4UB8DrVQZEjLJRoc/57Vr1y4GDBhAtWrVeOutt8jOzg50SD6RmZnJyJEj6d27N5GRkYEORynlpZBJNk2bNmXnzp3k5OQEOpSAysrK4qWXXuL666+nd+/erFixglmzZtG1a1fWrl0b6PBK5Oeff6Z9+/Y0aNCAGTNmBDocpVQRhEyyqVSpEtWrV2fv3r2BDiVgkpKSaN++PcnJyaxevZrnnnuO1q1bs2jRIkaPHk2/fv34wx/+wLlzQfHcO69duXKFF154gXvvvZf333+ff/3rX5QvXz7QYSmliiBkkg2E7yCBI0eOMHz4cEaMGMGrr77Kt99+m+fJlSLCyJEj2bhxIwcPHqR169bMnTs3gBF7b8+ePfTo0YPk5GTWrl3LwIEDAx2SUqoYQirZhNsggZycHCZPnkzr1q2pWbMmmzZtYvBgz09+iIuLY9q0abz//vs88cQTDB06lMOHD5dixEUzY8YMunTpwl133cXcuXOpVatWoENSShWTv6erKVXhVLNZv349jz32GAALFy6kbdu2Xu/br18/Nm7cyIQJE2jTpg2vvPIKI0eOJCIiOH57nD17lqeeeoply5bxww8/cN111wU6JKVUCQXH1cVHwmFE2tmzZ3nuuee45ZZbePDBB1m6dGmREo3D1VdfzWuvvcb8+fOZNGkSiYmJbNmyxQ8RF01KSgodO3ZERFizZo0mGqVCREglm1BvRps9ezatWrXiyJEjbNy4kUceeaTEtZH27dvzyy+/8Jvf/IabbrqJ8ePHB+QR2zk5Obz++usMHDiQv/3tb3z88cdUqlSp1ONQSvlHSEzE6ZCdnU3FihU5ceIEFSpUKKXI/C89PZ2nn36abdu28d5779GrVy+/nGffvn08+eSTbNu2jQ8++ICePXv65TyuDh48yIgRIzh37hzTpk2jQYMGpXJepcKBTsTpB5GRkTRq1IgdO3YEOhSfuHz5Mq+//jrXXXcdXbp0Yd26dX5LNAD16tVj1qxZvPLKK9x///08/PDDZGZm+u18AP/973/p2LEjN9xwAz///LMmGqVCVEglGwidQQLLli2jY8eO/Pjjj6xcuZLnn3+e6Ohov59XRBgyZAibNm0iOjqaVq1aMWPGDHxdA7548SK///3vefzxx/n888+ZMGEC5cqF1HgVpZSTkEs2Zb3f5vjx4zzyyCP89re/5cUXX2Tu3Lk0bty41OOoUqUK7777Ll9//TWvvPIKAwcOZM+ePT459pYtW+jatSsZGRmkpqZy0003+eS4SqngFXLJpqzWbIwxTJkyhVatWlGhQgU2b97M3XffjUhgm1pvuOEG1qxZQ48ePejcuTNvvPEGV65cKdaxjDFMmjSJHj168MQTT/DFF19QtWpVH0eslApGITVAAKzmp2effZYVK1aUQlS+sWXLFkaPHs3Zs2d5//336dSpU6BDcmvnzp089thjHD9+nEmTJtG5c2ev983MzOSRRx5h165dzJgxQ2dsVqqU6AABP3E0o5WFJHr+/Hmef/55evTowd13383KlSuDNtEANGnShAULFjB27FgGDRrEmDFjOHPmTKH7OSbQrF+/PitWrNBEo1QYCrlkU6NGDUSEY8eOBTqUAs2dO5fWrVuze/du1q9fzxNPPFEmpswXEYYNG8bGjRs5efIkrVq14rvvvnNbVifQVEo5hNzwHxHJrd3ExsYGOpx8MjIyGDNmDGvXruW9996jX79+gQ6pWGrUqMEnn3zCjz/+yGOPPcaUKVN4++23qVOnDgBpaWncd999VK5cmbVr1+q8ZkqFuZCr2UBwDhK4cuUKb7/9Nu3bt+faa69lw4YNZTbROOvduzfr16+nRYsWtGvXjvfee08n0FRK5RNyNRsIvjnSUlNTeeihh6hSpQpLliyhRYsWgQ7JpypUqMBLL73Evffey6OPPkpmZiZz587Vec2UUrn8OhpNRKKBxcBVWIntS2PMBJcyI4DXgf32qonGmI/dHMur0WgAX375JdOmTeObb74pSfg+YYyhRYsWjB07llGjRgV8KLNSKrwEy2g0v9ZsjDFZItLLGHNeRCKBZSIy1xiT7FJ0pjHmaV+dN5hqNosXL6ZcuXKaaJRSYc3vfTbGmPP2YjRWcnNXPfHpVbhJkybs3r272Dcf+tLkyZN55JFHNNEopcKa35ONiESIyFrgELDAGLPKTbEhIpIqIp+LSN2SnrNChQrUqlWL9PT0kh6qRDIzM/n+++8ZNmxYQONQSqlA8/sAAWNMDtBBRGKAWSLS0hiz2anIt8B0Y8xlERkFTAF6uzvW+PHjc5cTExNJTEz0eF7H8OdAzCvm8J///Idbb72V6tWrBywGpVR4SUpKIikpKdBh5FOq09WIyAvAOWPMmx62RwCZxphr3GzzeoAAwFNPPUXjxo0ZM2ZMseMtCWMMbdu25Z133ikwKSqllD8FywABvzajiUgNEaliL1cA+gBbXco434QxGHCu9RRboAcJrFy5kqysrFJ7AJlSSgUzf/fZ1AYWiUgqsBKYZ4yZIyITROQ2u8zTIrLR7td5EnjAFycO9KMGJk+ezMMPP6wDA5RSihCc9dkhPT2dbt26sX///sIL+9jp06dJSEhg69at1KxZs9TPr5RSDmHRjBZI9erVIzMzk7Nnz5b6uadPn07v3r010SillC1kk01ERARNmjRhx44dpX5ux701SimlLCGbbCAw/TZr1qzh+PHj9OnTp1TPq5RSwSykk00gZn92DAyIiAjpj1YppYokJGd9dmjevDnz5s0rtfOdO3eOzz77jA0bNpTaOZVSqiwI6Z/fpd2M9vnnn9O9e3fi4+NL7ZxKKVUWhHyy2b59O6U1vHvSpEk6MEAppdwI6WRTrVo1oqOjOXTokN/PtXHjRvbt28eAAQP8fi6llCprQjrZQOkNEpg8eTIPPvgg5cqFdDeYUkoVS1gkG3/321y8eJFp06bx0EMP+fU8SilVVoV8snH02/jTV199RadOnWjQoIFfz6OUUmVVyCeb0qjZ6IwBSilVsJBPNv4e/rx9+3a2bt3KoEGD/HYOpZQq60I+2TRu3Ji9e/dy+fJlvxx/8uTJjBgxgquuusovx1dKqVAQ8kOnoqOjiY+PZ/fu3TRv3tynx7506RKffvopS5cu9elxlVIq1IR8zQb8N/x59uzZtGzZkqZNm/r82EopFUrCJtn4o99GBwYopZR3wiLZ+GP48549e1i7di1Dhgzx6XGVUioUhUWy8UfN5qOPPuJ3v/sd5cuX9+lxlVIqFIX8AAHw/fDnK1eu8O9//5sFCxb47JhKKRXKwqJmEx8fz9mzZzl16pRPjjdnzhwaNGhAy5YtfXI8pZQKdWGRbESEpk2b+qzfRh8loJRSRRMWyQZ8N/x5//79LF++nLvvvtsHUSmlVHjwa7IRkWgRWSkia0Vkg4iMK6DsXSKSIyId/RGLrwYJfPzxx9x7771UrFjRB1EppVR48GuyMcZkAb2MMR2A9sAAEeniWk5EKgFPAyv8FYsvhj9nZ2fz0UcfMWrUKB9FpZRS4cHvzWjGmPP2YjTW6Dd3z2j+G/AqkOWvOHxRs1mwYAFxcXG0b9/eR1EppVR48HuyEZEIEVkLHAIWGGNWuWzvANQ1xsz1ZxyOmk1OTk6xj6EzBiilVPH4/T4bY0wO0EFEYoBZItLSGLMZQEQEeBMY4bSLeDrW+PHjc5cTExNJTEz0Oo6YmBhiYmI4cOAAdevWLdJ7ADh06BA//fQTn3zySZH3VUqp0pKUlERSUlKgw8hHjHHXquWnk4m8AJwzxrxpv44BdgJnsZJMLeA4cLsxZo3LvqaksSYmJvLCCy/Qu3fvIu/76quvsnPnTj788MMSxaCUUqVJRDDGePwRX1r8PRqthohUsZcrAH2ArY7txpjTxpg4Y0wjY0xDrAECg1wTja8Ud5BATk4OH374oTahKaVUMfm7z6Y2sEhEUoGVwDxjzBwRmSAit7kpbyigGa2kijtIICkpiauvvpouXfINpFNKKeUFv/bZGGM2APnumzHGuL3fxhhzsz/jadasGQsXLizyfo6BAVYXk1JKqaIq1T6bkvBFn8327dsZMGAAu3bt8nqfY8eO0bRpU3bv3k3VqlVLdH6llCptYdFnE2waNmxIRkYGWVne384zdepUBg0apIlGKaVKIKySTVRUFAkJCV7XbIwxem+NUkr5QFglGyjaIIFly5ZhjKF79+5+jkoppUJb2CWbogx/njx5Mg8//LAODFBKqRIKu2Tjbc3m5MmTzJ49m+HDh5dCVEopFdrCLtl4+4joadOm0a9fP2JjY0shKqWUCm1hl2y8eYiaY2CAPkpAKaV8I+ySTc2aNbl06RKZmZkey6SkpHDmzBl69epVipEppVToCrtkIyKFDhJwDAyIiAi7j0cppfwiLK+mBQ0SOHPmDF988QUPPPBA6QallFIhLCyTTUE1m5kzZ5KYmEjt2rVLOSqllApdYZlsCqrZ6IwBSinle2GZbDwNf163bh2HDh2iX79+AYhKKaVCV9gmm127dpGTk5Nn/eTJkxk5ciSRkZEBikwppUKTX59nE6wqVqxI9erV2bt3Lw0aNADg/PnzzJgxg7Vr1wY2OKWUCkFhWZrvb6cAAAikSURBVLOB/IMEvvzyS7p27Ur9+vUDGJVSSoWmsE02roMEdGCAUkr5T9gmG+eazebNm9m1axe33nprgKNSSqnQFLbJxrlm8+GHH/LAAw8QFRUV4KiUUio0heUAAfh1+HNWVhZTp05lxYoVgQ5JKaVCVtjWbBo0aMCRI0eYPn067dq1o3HjxoEOSSmlQpZfk42IRIvIShFZKyIbRGScmzKPish6u8xiEWnhz5gcIiMjadSoEePGjdNHCSillJ/5NdkYY7KAXsaYDkB7YICIdHEpNs0Y09Yu8zrwL1+dPykpqcDtzZo148KFCwwePNhXp1RKqaBQ2PWvtPm9Gc0Yc95ejMbqIzIu2886vawE5L2tvwQK+7DbtWvHww8/THR0tK9OqZRSQSHYko3fBwiISASwGmgMvGuMWeWmzOPAWCAKuNnfMTmMGzcOY0zhBZVSSpVIadRscuwmsrrA9SLS0k2Z/zPGNAH+BLzg75gcREQfkKaUUqVASvOXvYi8AJwzxrzpYbsAJ4wx17jZplUQpZQqBmOMBDoGvzajiUgN4LIx5pSIVAD6AK+6lGlijNlpv7wNcPtUs2D4sJRSShWPv/tsagNT7H6bCOAzY8wcEZkArDLGfA88KSK3AJeAE8AIP8eklFKqlJVqM5pSSqnwVOq94yLykYgcFpH1LuufEpEt9s2frk1tKSISJSIvicheETntsv0ZEdkkIqkiskBE6rlsnyMidUTkPyKy1b6J9EMRiXQq87aI7LCP0d5p/QgR2S4i20RkuNP6ofZxUu3jV/PVZ6SUCk0i0t++Bm0XkT/a6xaLyBr7xvYMEfnaqXw5EVltL3u6dlYVkfn2NWqeiFRx3V9E4kXkJ/s6uUFEni5sfxFpLiLLReSiiIx1Kt/MjtUR8ynn43lkjCnVP6A71g2e653WJQLzgXL26xpO2xoAs+zlLkBN4LTLMXsC5e3lx4CZTtvKAyvs5f5O66cDj9rLA4D/2svXO5WvCuwCqgDXOC1HAoeBqna514AXS/uz1D/907+y84f1434nkIB1m0cq0MKlzJfA75xeJwL/ay/nu3ba618D/mgv/wl41XV/+7rZ3l5XCdjmOLen/YFY4Drgb8DYAt7TAaBeYe+/1Gs2xpilWH0zzkZjvcErdpljTtv6Az/Y65ONMYfdHPNnY8xF++UKIN5pcyKQZJf7wWl9MtZwbIDBwKd2mZVAFRGpCfQD5htjThljTmIlxP6AY7BCZXsEXQzWB66UUp50AXYYY9KNMZeBmVjXHgBEJAbrPsNZTvv0B+aCx2sn9jGm2MtTgDtc9zfGHDbGpNrHOQts4dfrpNv9jTFHjTGrgSsFvKdbgF3GmH0FlAGCZyLOZkAPEVkhIotEpJPTttxk46WHsL8c2wDX/UWkHDDMqVw84Pxh7bfXua7PAOLtpPg4sMEuey3wURFiVEqFH0/XGYfBwEKTd1aVXtg/lgsQ5/gRbow5hFWL8bi/iDTAqiE5prp33T+u0Hfyq98CM7wpGCzJphxWk1RX4I/A5wAiEoV1cU/z5iAi8jusat/rTqu7AUtdiv4f8LMxZrmnQxVynnJYtbF2xph4rKTzV29iVEopD4bidOEWkTrAcadWG2/leNpfRCphNdX93hhzzsP+Xo0as6/PtwNfeFM+WJLNPuBrAGNNZ5MjItWBm8ifKNyyh0//BRhkV1ERkYbAXkfznL3uRaw+obFOu2cAzoMK6trrMoD6bta3t2NNs9d/DtzgTZxKqbDl6XqCfb3rDPzXaXt/YJ4Xxz1sN/sjIrWAI+72t38kfwlMNcbM9mL/wgwAVhtjjnpTOFDJRshbe5iFPSeaiDQDoowxx3Fqr3Sz/68vRDoA7wO32/s55GlCE5GHsfphhroc71tguF2mK3DSrlbOA/qISBURqYp1U+o8rH8g19r/QLDXb/HurSulwtQqoImIJIjIVcC9WNcegLuB740xl5zKu7v+uV47sY/xgL38AOBIJK77fwxsNsb8bwH7j3Da3/W8rvLUxAoVgBEZ07E607OAvcCDWM1oU7Gao1KAnnbZZCDaZdTFPqwOq73YI8CABcBBYA2wll9Hr30L1Hfa/zKwwy6zBvgfp20TsUaKrAM6Oq1/wN5nOzDcaf0oYDPWiJLZ2CPT9E//9E//PP1hJYBt9jXlz07rfwL6Or2OwKo1OO+b79ppr68GLLSPOx9r5Gye/bG6E7Lt65Xj+tff0/7/v707eLEpDOM4/v0ZTciM1WRhITLZYKGsKDslZcWGshuLiRR/wOwoG4koK9lY+QNkJxZKo+QfsJOpSSYpjMfivjMxdZnUce+Y72d17znnvfc9m/vr3POe52nbt7ff2w/AfPvOrW3fFmAOGFvtuQ/tQ51JdgD3qurEX44fBZ5V1cr+OZI01JIcBs5W1fQgxndhaMNGkvT/GJYFApKk/5hhI0nqnGEjSeqcYSNJ6pxhI0nqnGEj9ZFksZVRf9NKqV9uhVd/N2ZnkpUPDUvrnmEj9fepqg5W1T56VSKOAzN/GLMLONP5zKQ1xrCRVqF6bS/OAxdg+QrmaWvs97KVOQK4BhxpV0SXkmxIcj3Ji9Zob2pQ5yANkg91Sn0k+VhV4yu2zQN7gQXge1V9SbIHeFhVh5IcBa5U1cl2/BQwUVVXW1WL58Cpqnr7b89GGqyNg56AtMYs3bMZBW63FuKLwGSf448B+5Ocbu/H27GGjdYVw0ZapSS7gW9VNZdkBnhXVQeSjACf+w0DLlbVk382UWkIec9G6m955VmSCeAucKtt2kav0jj02lOMtNcLwNhPn/EYmG69REgymWRzl5OWhpFXNlJ/m5LM0vvL7CvwoKputH13gEdJztHrmbTU9fA1veZ/r4D7VXWzteGdbcum3/Nrj3hpXXCBgCSpc/6NJknqnGEjSeqcYSNJ6pxhI0nqnGEjSeqcYSNJ6pxhI0nqnGEjSercDxASt2Ek21ljAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f612144c240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X[::1000],Y[::1000], label=\"5000 rating window\", color='k')\n",
    "plt.xticks([X[600], X[-350]], [x[wSize+600].strftime(\"%d/%m/%Y\"), x[-350].strftime(\"%d/%m/%Y\")])\n",
    "plt.xlim(X[0], X[-1])\n",
    "plt.ylabel(\"Rating\")\n",
    "plt.xlabel(\"Date\")\n",
    "plt.legend(loc=\"best\",fontsize=8)\n",
    "plt.title(\"Ratings over time (sliding windows)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FPMC in Tensorflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse(path):\n",
    "    g = gzip.open(path, 'r')\n",
    "    for l in g:\n",
    "        yield eval(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the interaction data, including the timestamps associated with each interaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "userIDs = {}\n",
    "itemIDs = {}\n",
    "interactions = []\n",
    "interactionsPerUser = defaultdict(list)\n",
    "\n",
    "for d in parse(dataDir + \"goodreads_reviews_comics_graphic.json.gz\"):\n",
    "    u = d['user_id']\n",
    "    i = d['book_id']\n",
    "    t = d['date_added']\n",
    "    r = d['rating']\n",
    "    dt = dateutil.parser.parse(t)\n",
    "    t = int(dt.timestamp())\n",
    "    if not u in userIDs: userIDs[u] = len(userIDs)\n",
    "    if not i in itemIDs: itemIDs[i] = len(itemIDs)\n",
    "    interactions.append((t,u,i,r))\n",
    "    interactionsPerUser[u].append((t,i,r))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interaction with timestamp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1386269065, 'dc3763cdb9b2cae805882878eebb6a32', '18471619', 3)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "interactions[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "542338"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(interactions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sort interactions by time (including interaction sequences for each user). Useful when building data structures that include adjacent pairs of interactions (but consider whether this is desirable if making train/test splits!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactions.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "itemIDs['dummy'] = len(itemIDs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a data structure including users, items, and their previous items"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactionsWithPrevious = []\n",
    "\n",
    "for u in interactionsPerUser:\n",
    "    interactionsPerUser[u].sort()\n",
    "    lastItem = 'dummy'\n",
    "    for (t,i,r) in interactionsPerUser[u]:\n",
    "        interactionsWithPrevious.append((t,u,i,lastItem,r))\n",
    "        lastItem = i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "itemsPerUser = defaultdict(set)\n",
    "for _,u,i,_ in interactions:\n",
    "    itemsPerUser[u].add(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "items = list(itemIDs.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the tensorflow model. Similar to models from Chapter 5, with the addition of the term associated with the previous interaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = tf.keras.optimizers.Adam(0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FMPC class. UI and IJ are given as initialization options, allowing us to exclude certain terms (for exercises later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FPMC(tf.keras.Model):\n",
    "    def __init__(self, K, lamb, UI = 1, IJ = 1):\n",
    "        super(FPMC, self).__init__()\n",
    "        # Initialize variables\n",
    "        self.betaI = tf.Variable(tf.random.normal([len(itemIDs)],stddev=0.001))\n",
    "        self.gammaUI = tf.Variable(tf.random.normal([len(userIDs),K],stddev=0.001))\n",
    "        self.gammaIU = tf.Variable(tf.random.normal([len(itemIDs),K],stddev=0.001))\n",
    "        self.gammaIJ = tf.Variable(tf.random.normal([len(itemIDs),K],stddev=0.001))\n",
    "        self.gammaJI = tf.Variable(tf.random.normal([len(itemIDs),K],stddev=0.001))\n",
    "        # Regularization coefficient\n",
    "        self.lamb = lamb\n",
    "        # Which terms to include\n",
    "        self.UI = UI\n",
    "        self.IJ = IJ\n",
    "\n",
    "    # Prediction for a single instance\n",
    "    def predict(self, u, i, j):\n",
    "        p = self.betaI[i] + self.UI * tf.tensordot(self.gammaUI[u], self.gammaIU[i], 1) +\\\n",
    "                            self.IJ * tf.tensordot(self.gammaIJ[i], self.gammaJI[j], 1)\n",
    "        return p\n",
    "\n",
    "    # Regularizer\n",
    "    def reg(self):\n",
    "        return self.lamb * (tf.nn.l2_loss(self.betaI) +\\\n",
    "                            tf.nn.l2_loss(self.gammaUI) +\\\n",
    "                            tf.nn.l2_loss(self.gammaIU) +\\\n",
    "                            tf.nn.l2_loss(self.gammaIJ) +\\\n",
    "                            tf.nn.l2_loss(self.gammaJI))\n",
    "\n",
    "    def call(self, sampleU, # user\n",
    "                   sampleI, # item\n",
    "                   sampleJ, # previous item\n",
    "                   sampleK): # negative item\n",
    "        u = tf.convert_to_tensor(sampleU, dtype=tf.int32)\n",
    "        i = tf.convert_to_tensor(sampleI, dtype=tf.int32)\n",
    "        j = tf.convert_to_tensor(sampleJ, dtype=tf.int32)\n",
    "        k = tf.convert_to_tensor(sampleK, dtype=tf.int32)\n",
    "        gamma_ui = tf.nn.embedding_lookup(self.gammaUI, u)\n",
    "        gamma_iu = tf.nn.embedding_lookup(self.gammaIU, i)\n",
    "        gamma_ij = tf.nn.embedding_lookup(self.gammaIJ, i)\n",
    "        gamma_ji = tf.nn.embedding_lookup(self.gammaJI, j)\n",
    "        beta_i = tf.nn.embedding_lookup(self.betaI, i)\n",
    "        x_uij = beta_i + self.UI * tf.reduce_sum(tf.multiply(gamma_ui, gamma_iu), 1) +\\\n",
    "                         self.IJ * tf.reduce_sum(tf.multiply(gamma_ij, gamma_ji), 1)\n",
    "        gamma_uk = tf.nn.embedding_lookup(self.gammaUI, u)\n",
    "        gamma_ku = tf.nn.embedding_lookup(self.gammaIU, k)\n",
    "        gamma_kj = tf.nn.embedding_lookup(self.gammaIJ, k)\n",
    "        gamma_jk = tf.nn.embedding_lookup(self.gammaJI, j)\n",
    "        beta_k = tf.nn.embedding_lookup(self.betaI, k)\n",
    "        x_ukj = beta_k + self.UI * tf.reduce_sum(tf.multiply(gamma_uk, gamma_ku), 1) +\\\n",
    "                         self.IJ * tf.reduce_sum(tf.multiply(gamma_kj, gamma_jk), 1)\n",
    "        return -tf.reduce_mean(tf.math.log(tf.math.sigmoid(x_uij - x_ukj)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelFPMC = FPMC(5, 0.00001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trainingStep(model, interactions):\n",
    "    with tf.GradientTape() as tape:\n",
    "        sampleU, sampleI, sampleJ, sampleK = [], [], [], []\n",
    "        for _ in range(100000):\n",
    "            _,u,i,j,_ = random.choice(interactions) # positive sample\n",
    "            k = random.choice(items) # negative sample\n",
    "            while k in itemsPerUser[u]:\n",
    "                k = random.choice(items)\n",
    "            sampleU.append(userIDs[u])\n",
    "            sampleI.append(itemIDs[i])\n",
    "            sampleJ.append(itemIDs[j])\n",
    "            sampleK.append(itemIDs[k])\n",
    "\n",
    "        loss = model(sampleU,sampleI,sampleJ,sampleK)\n",
    "        loss += model.reg()\n",
    "    gradients = tape.gradient(loss, model.trainable_variables)\n",
    "    optimizer.apply_gradients((grad, var) for\n",
    "                              (grad, var) in zip(gradients, model.trainable_variables)\n",
    "                              if grad is not None)\n",
    "    return loss.numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run 100 batches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration 10, objective = 0.48359895\n",
      "iteration 20, objective = 0.45931146\n",
      "iteration 30, objective = 0.45187306\n",
      "iteration 40, objective = 0.45339736\n",
      "iteration 50, objective = 0.4563751\n",
      "iteration 60, objective = 0.45513394\n",
      "iteration 70, objective = 0.4522531\n",
      "iteration 80, objective = 0.45192516\n",
      "iteration 90, objective = 0.45243126\n",
      "iteration 100, objective = 0.4488452\n"
     ]
    }
   ],
   "source": [
    "for i in range(100):\n",
    "    obj = trainingStep(modelFPMC, interactionsWithPrevious)\n",
    "    if (i % 10 == 9): print(\"iteration \" + str(i+1) + \", objective = \" + str(obj))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(still using the hourly bikeshare interaction data from examples above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "xsort = list(hourly.keys())\n",
    "xsort.sort()\n",
    "X = [feature(x) for x in xsort]\n",
    "y = [hourly[x] for x in xsort]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X, y)\n",
    "theta = model.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "ypred = model.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = y - ypred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEZCAYAAACq1zMoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXn8FlP7xz9Xi9ImtEh7UZZQKUpUtFA8ylPZKYkk28+aPJ7Kg+KxRb4iFT2WQiRERb7taSdatWkv2tPe9ftjZu7mO/fs98w9c9/f6/16fV/fuc+cOeeamTPnOst1rkPMDEEQBEEImgJRCyAIgiBkJ6JgBEEQhFAQBSMIgiCEgigYQRAEIRREwQiCIAihIApGEARBCAVRMIJviOhXImpqca4ZEa0LKJ8fiairz2unEdEFQciR7RDRpUS0JGo5vOClnBFRHyL6n3pcjogWE1HhcCXM34iCyQcQ0Roi+puIdhPRRiIaTkTFUk2Xmesw8xS7KKnmkQpEdA2A3cz8cwBpdSaiqQGIFVuYeRoznx21HD7wUs4YAJh5K4BJALqHIpEAQBRMfoEBXM3MpQDUBVAPwJPRipQW7gHwv4DSIqSgMImoYEByhJJePuUjiIIJFVEw+QcCEi238VAUjXKC6AQieomI1hLRJiLKIaIi6rlTiegrItpBRH8R0WTddauJ6Ar1uCgRvUdE24noVwAN82ROdIyIauh+DyeiZ9Tj0moeW9U8viKiiqY3QVSTiHKJaKca/2OLeIUBXAFAL29DIpqh3ssGInqDiAqp56qqMhbQxf+RiLoS0VkA3gLQmIj2ENF29XwpIhqhyrGaiJ7SXdtZHZ57hYj+BNBHDe+qDs38RUTfElEV3TWtiWipKt+b6n12TSG9V4loCxHtIqKfiegcNbwtEf2m9mjXEdHDanhiuImIHieiTw3PdCARvaa793fVHvE6IvoPEZHFu+hDRJ8Q0f/UPH8mojOJqJcq31oiaqmLX4GIvlTvaTkRddOdcypnFYjoM/WdrCSi+81kUvkJQA0iqmwTR0gBUTD5DCKqBKANgBW64BcAnAHgfPV/RQD/Vs89AmAdgFMBlAPQ2yLpvgCqq39XAuhsOG/X+i8AYBiAygCqAPgbwCCLuP8BMJ6ZSwOoBOANi3hnAjjKzBt1YUcBPATgFACNoSige51kZOalUHpDM5m5JDOfop4aBKAkgGoAmgO4nYju0F16MYDfoTy354ioHYBeANoDKAtgKoCPAYCIygD4FMATUJ71MlVG+EyvNYBLAZzBzCcBuB7AX2o67wK4S+3R1oEyVGR8BiMBtCGi4mp6BQB0AvChev59AIcA1IDSI24FoBusuUa9pjSAhVAaOQTgdCjv9B1d3FEA/gBwmprn80TUXD3XFxblTFVwXwFYAKACgBYAHiSiVmYCMfNRKM9T5ujCgpnlL8v/AKwGsFv9OwZgIoBSuvN7AVTX/W4MYJV63A/AFwBqWqR7hXq8EkAr3bm7APyh+30MQA3d7+EAnrGQty6Av3S/fwTQVT1+H8BgABUd7vkSABsd4jwIYLR6XBWKAipgkW9nAFN05woAOAigti7sbgCTdPHXGPIbB+AOQxr7oCjW2wBMN8T/w5C/l/QuB7AUilIiw3Vr1PdT0hDezPDOpgC4VT1uBWCFelwewAEARXRxb9Tu3eQ594HSKNB+X6OWRVJ/l1CffSlV9sMAiuniPw9gmFM5U+/V+Ix6ARiqk2OE4fw07R7lL/g/6cHkH9qx0mJtBuAsAGUAgIjKAigGYJ467LAdwLdQWtEA8F8oH/UEIvqdiJ6wSP90AOt1v9e6FYyITiSit0kxRtgJZVirtMWQy2NQKtLZRLTI0GPQswNK70Kfz5nq8NsmNZ/noD4HH5QBUAiKEtBYC6X3p2G0bqoKYKDuOf8FpcdQEcrzM8Zfb/jtOj1m/hFKD+tNAFuIaDARlVCv6wDgagBr1WHARhb3+DGAm9Tjm6DMWQBKL7MwgE1q3jugKH27Z7lFd7wfwJ+s1vDqb4KiaCoA2M7Mf+vi65+rXTmrAqCi9jxUuZ6E0uOzoiSAnTbnhRQQBZN/0OZgpkLpBbyshv8JZUjqXGY+Rf0rzcqwCph5LzM/ysw1AVwL4GEiutwk/U1QWp8aVQ3n/4aiyDRO0x0/CmVIqyErQ1+a6XOSgmHmrcx8NzNXhDJslUO6uR0dv0MZNamgC3sLwBIovbHSAJ7S5bFP/W8lo3H47E8oLW39fVYFsMHmmj8AdNc955OZuQQzz0Ly8wOUIUA9XtIDMw9i5gYAzgFQG4pyBjPPY2ZtWO1LAJ/AnE8BNCdlPuw6HFcw66D0YE7V5Vuamc+3SMcLGwGcog3NqVTB8edqV87WQel565/HScz8D7OMSDGUOANAylaGgjmiYPInrwFoRUTnqa3IIQBeU3szIKKK6hg+iOhqIqqpXrcHwBEowxlGPgHwJCkT9pUA3Gc4vwDAzURUgIiugtKT0igBpRW7m4hOgTLObgoRdaTjBgA7oQy9HTPGY+bDAL435FMSitny36RM3PfQxf8TSiV2qypjVwA1ddduAVCJ1HUTzHxMvefniKgEEVUF8H+wt1p7G0Bv3WT7SUTUUT33DYA6RHQtERUkovugDEXZYZkeETUgootIMWLYD0UhHCOiwkR0MxGVYmUOYg/M36f2TCZDGc5cxczL1PDNACYAeJWISpJCDbJYE+UFZl4PYAaA/kRUhIjOB3Anjj9Xu3I2G8AeUgwUiqrP8VwiamCR3UUAVjNzIOu1hGREweQP8rR81YrjfRyfyO8FpcU/Sx06mgCglnruTADfE9EeANMBvMnH177o0+0HpUW9GsB3AEYYZHgISg9oB5Thli90516D0nP4E0rlMs5G/oYAfiKi3QDGAHiAmddY3Pc7AG7X/X4UwC3qtW9DmcjWcxeAx1U5zlbvV2MSgN8AbCairWrYA1B6ZqugzFd8wMzDLWQBM48BMADASPU5/wLgKvXcX1AmtP+r5n8WgLlQ5nk8pwdlPmMIgO1Q3smfatqAMt+zWr3mbgA3W+UBpdfSAscn9zVuB3ACgMVqHp8ib4/PK/p3fBOUSfyNAEYDeFod8gNsypmq9K+BMoe3GsBWKM+glEWet0AZ2hNCQptki04AoqFQCsUWrYtNRH2gfOzah9ybmb9Tzz0JoCuUlvSDzDwh/VILmQIpiyPv4wAWW6YTdf5pPYCbmXmyU3zBG2pvPRdAPWY+FLE4WUscFMylUKyYRhgUzB5mfsUQ92woLaqGUManvwdwJkd9E4IQAOqw5E9QhrMegzKEV4OZLXsxghBnIh8iY+ZpUIZNjJhZELUDMJKZj6jDIiugjKMKQjbQGIrF3lYoVl7tRLkImUzkCsaGnkS0kJTVwiepYRWR11RzA/KahQpCxsLM/Zi5jGr51JiZ50YtkyCkQlwVTA4UU9K6ADbjuEmtIAiCkCEUiloAM5h5m+7nECjuHwClx6K3ga+EvOsOAABEJHMygiAIPmBmU59yfohLD4agm3MhIr254z8B/KoejwVwIynOGatDWSQ12yzBqF0kmP316dMnchlEJpEpP8olMrn7C5rIezBE9BEUR4GnEtEfUPwFXU5EdaEsoFsD1aU2My8mok+g2N4fBnAvh/FUBEEQhJSJXMEws9kiL7vFav0B9A9PIkEQBCEI4jJEli9o3rx51CIkITK5Q2RyTxzlEpmiIfKFlmFARDJyJgiC4BEiAmfhJL8gCIIjRIQpU6Y4RxRigSgYQRAyigULFkQtguASUTCCIAhCKIiCEQRBEEJBFIwgCIIQCqJgBEEQhFAQBSMIQkah7MUmZAKiYARBEIRQEAUjCIIghIIoGEEQBCEURMEIgpBRiBuozEEUjCAIghAKomAEQcgoxIoscxAFIwiCIISCKBhBEAQhFETBCIIgCKEgCkYQBEEIBVEwgiBkFDLJnzmIghEEQRBCQRSMIAgZhSy0zBxEwQiCIAihIApGEARBCAVRMIIgZBQyyZ85iIIRBEEQQkEUjCAIghAKomAEQRCEUBAFIwiCIISCKBhBEDIKmeTPHETBCIIgCKEgCkYQhIxCVvJnDqJgBEEQhFAQBSMIgiCEgigYQRAEIRREwQiCkFGIFVnmELmCIaKhRLSFiH7RhZ1MRBOIaBkRjSeik3TnXieiFUS0kIjqRiO1IAiC4ETkCgbAcABXGsJ6AfiemWsDmATgSQAgojYAajLzmQC6AxicTkEFQRAE90SuYJh5GoAdhuB2AN5Xj99Xf2vhI9TrfgJwEhGVT4ecgiAIgjciVzAWlGPmLQDAzJsBaEqkIoB1ungb1DBBEAQhZhSKWgCXeF5Z1bdv38Rx8+bN0bx58wDFEQQhKmSSPzhyc3ORm5sbWvpxVTBbiKg8M28hotMAbFXDNwCorItXSQ1LQq9gBEHIHmQlf3AYG9/9+vULNP24DJGR+qcxFkAX9bgLgC914bcDABE1ArBTG0oTBEEQ4kXkPRgi+ghAcwCnEtEfAPoAGADgUyLqCmAtgOsBgJnHEVFbIvodwD4Ad0QjtSAIguBE5AqGmW+2ONXSIv59IYojCIIgBERchsgEQRBcIZP8mYMoGEEQMoqtW7c6RxJiAWWjRQYRcTbelyDkd7Tei3zf4UBEYObAuojSgxEEQRBCIWsVzP79+6MWQciHrFy5Eo888kjUYghCLMhaBSPjtEIUjBw5Eq+88krUYghCLMhaBSMIgiBEiygYQRAEIRREwQhCgIh1kyAcRxSMIAiCEAqiYAQhQGSVuSAcRxSMIASIDJEJcWHWrFmYMGFCpDJE7uwyLKQlKQhCfqZdu3bYunVrpI2erO3BSEtSEAQhWrJWwQiCIAjRIgpGEPIBubm5GD16dNRiuObAgQMYMmRI1GIIKZK1CkbmYAThODfffDM6duwYtRiumTx5Mu6+++6oxcho4jBNkLUKRhCE48ShshHyH1mrYOSDEqJAyl0wyHPMDrJWwQiCkSlTpqB169ZRiyEI+YasVTAyByMY+eqrrzBx4sRQ84hruQurR/D555/j3XffDTxd6cFkB1mrYITs4tChQ5gxY0ZKaRQsWDAgaazJbxVj9+7dcdddd0UthmBCHMqiKBghI/jggw/QpEmTlNJIh4KJK1plM3z48IglcUccKkchdUTBZDhr167Fk08+GbUYoXP48OGU04jr8FU6ee2111zF++KLL1ChQoWQpclsmjRpEki5DIs4lPesVTDffPNN1CKkhZEjR2LAgAFpz3P58uVpzTPTWrTpfj5OaM9v6dKlruJPmTIFmzdvdp1ufmTGjBnYvXt31GLEmqxVMPPnz49ahKzlpptuwlNPPRW1GLFm3759ruOOGDEC77zzTojSHFcEhw4dCjWfoMjPiiso4vAMs1bBCEKm0L17d3Tv3j1qMfIQh8opE4jDMFScyVoFk19efLorghUrVqQ1PyEYwionmZaukF6yVsEI4TB27NhI8p0+fXrKaaSj0SEVYzJ//fUXdu7cGbUYQgRkrYJZvXo12rRpE7UYWUe6K9BDhw6BiPDBBx+kNd90ErTiW7duHV544YU8YV7fW5DvuWbNmmjatGlg6QmZQ9YqmBkzZuC7776LWgwhRQ4ePBi1CJ6Iw9Dse++9h169euUJi7JntWvXLqxbty6y/MMkDu/bijj0prNWwZx44olRi5DVRFV4V69eHUm+bonDRx0E2XIfCxYsSOrN5Te2bNkSWd5Zq2AKFAjv1u6//35s27YttPSF4xhbiKtWrQosrWzFTDmkQ2E0adIEv/zyi+k5r8/eq7zTp0/HvHnzksJfeumlpN5cfmPv3r2R5S0KxgeDBg1Cbm5uaOlnAvmlsk4Hmfws9YpgxowZ+PHHH23jr1+/Hps2bQpcjksvvRQtW7YMPN1MJg7lKmsVTBwebjrIlqEMt2TKe7WSc9u2bZEM2cSlnNSqVQsXXnihY7y4yOtEppTHqIi1giGiNUT0MxEtIKLZatjJRDSBiJYR0XgiOikK2TLlAwiaXbt2Aci/9+8Wq+fz2WefoVevXli2bFks5OnQoUOgizytKlwtfP/+/dixY0dg+TnlnZ8VQBy+0VgrGADHADRn5nrMfJEa1gvA98xcG8AkAKaeHvUFa968eXjooYfCljVf8Oyzz6Y1v2ytIPRWVVHe4+eff46RI0eGno++snNT8cWhchRSJ+4KhpAsYzsA76vH7wNob3qh7qMdNmwYBg4cGKxgMan44iJHukjlfuP+rIIyyQ5ikt9vBR/lM477+82PxF3BMIDxRDSHiLqpYeWZeQsAMPNmAOXMLtQK2/jx45GTkxO8YDFpYcVFjrDI1krD7L4WL14cWn7pKidOQ2TplCUdRL0A+NChQzjjjDMilcGOuCuYJszcAEBbAD2J6DIoSkePbWkdM2ZMWLIJMWb06NGYO3duZPlbVbRa5WpXyTZo0CAx15UtbN++3VN8P0ooijmYUaNGhZq+E3v27MHKlStNz7kpa2FTKLKcXcDMm9T/24hoDICLAGwhovLMvIWITgOw1exabZ+GsCqZoArun3/+iTJlygSSVn7A7XPv2LEj6tatiwULFoQskTmpfNTz5s3DmjVrsHXrVpx44om49NJL0y6P2/he0v3hhx88yRA3pk6dimLFirmygssUcnNzQ11yEdseDBEVI6IS6nFxAK0BLAIwFkAXNVpnAF+aXV+6dGkAQMOGDUORT/9h7du3z/dEadmyZT1bFDVt2hRdunTxlZ+QGkSEhQsXAgCeeOIJjBs3Ls95pwrXrpVtPNe6dWtcffXVqYjrCq+NpW3btvnayXHYsGGerwmTL7/8En/++afr+E2bNkXr1q1DlMg7bt6dXZzmzZujb9++ib+gia2CAVAewDQiWgBgFoCvmHkCgBcAtCKiZQBaADDdzlGz0nnrrbdCF/Szzz7DTTfd5Pt6ryttp06dim+//dZ3fpmE3cdx9913W64cN7s2qF7n77//DgB48cUX8eqrr7rKW8OogFasWJF4/7fddlsgwxlhr+QvV66cpTVhJs2ZtW/fHi+//HLUYliSk5OD999/3zlijImtgmHm1cxcVzVRPo+ZB6jh25m5JTPXZubWzOzaD3iQrf6oPySr/P/xj3/g008/TbM0yaxcuTLR0vfK+vXr0adPH8d4Q4YMsR0DD/IdMTO++uqrpHRTrbhr1aqVOP71119x991358kjqHuwk9PPPWzYsMH3tWGZKafjm9TyCGteo3bt2rjzzjsBAD179sR9993nO62o6yggxgomDIytgQcffND3FrJe7frtCLIgfP31147DdUWKFMEXX3wRWJ5mtGjRAvXq1fN17ahRo/DMM88ASH62ds+KmfHll6Yjpimzbt06XHvttQDyuiHyIp+eo0ePmpolv/vuuylIGT1RVmpmw11By/PXX38BAG6++WYAqX37/fr1SxpqXL58OaZOnepfQB1xsNazVTBEVN/uL11ChsXrr7+esW7E7VpSTgXr0KFDmD17dihyuZXBjkcffdR1XH0FsnHjRrRv3970XJA4KTk33HbbbahcuXJKeXmRwU4uszz8zCW5xUvZSMW5aZgEsTi1b9+++OOPP5LC46AYgsLJisxugJIBXBGgLGnl888/B+DtZd5zzz3YulUxWgtymCQOXdmgCeKe9u3b55iOn3yWL1+eZ2jKK/o8c3NzsXv3bpQqVcpTGtnqjdvufXz00Ueu09G+qZo1a6b0ffkpH0SEuXPnerIWIyIcOXIEBQsWdBXfblGt1f0eO3YMx44dQ6FCsTb+zYOtpMx8eboESTc//fST52s+/vjjhPlznFsZUch29OhRFCpUKNC8S5QogX379uUJ81JhWE3yp+q+XJ/usWPH8OKLL6bdhY4dfnokfvNxW6Ear0sXfntiS5Ys8eyU89ixY66fR9GiRV2lqbF//360atUKGzZsiP2eSHpcq0IiqgPgHACJJ8PMI8IQKh3EdRLRLanKkur1a9euzfP72LFjgabvhzDz1Kdt3AoiTKutdAyRuU3D7PzRo0c9petWlqgbcMby7IagZF61alVSg+Xxxx/H9OnTfaUX5bN0NclPRH0AvKH+XQ7gRQDXhihX2vDy8OO8uvqUU04JPQ/9mo90rZIPcy7AL+m0VrKiZcuWvir3VNm5095oM06NMMC/PHbl7pVXXvErjmuMa4aMnhAaN26M8ePHBzIfGCZurcg6QllzspmZ7wBwAYBI3OTHkTjMwYTlAl2PZj5phvEZrFmzJmRpFOyeXViVnd1mdumqYH/44QdLC0iz8mhnLelFZs2Kyorvv//eMY1M6MHY5T906FDP16Sav/EdzZo1K2mRrxGtTohS6btVMPuZ+RiAI0RUCop7FmcTmBgTdQFOlShaLm7SDbswBzEHk2q+dukGvU4iHe/ZrFfiN+3Ro0enKk4Se/bsSQpzu0AyjB6MlyUKbhSum/z192FmIBDX+sytgplLRKUBDAEwD8B8ADNDkyoG/P77767HYaPswejzfvJJ061xAiMOhbhv375o1apV4ndQcxReOHLkCJYsWWKbnjFs//79gSnfyZMnB5JOKli14t3gtQdj1hseMcLd9K/fd+1nDsbI9u3b85RVL2zevDnPb32v+V//+lfiWCtTdm54Yj8Hw8z3MvNOZh4MoBWAzupQWdZy5plnRuKKe9q0ab78PAHAgAGK1xy/BerIkSOYOnVqYpW2kTgomEmTJlm2CsNyDWNkyJAhOOecc2zjGJ+V3wW9ZoS9SBZI/dkFOV8ZhpfkdMztufle9u7da6rMjL0UvTxmC0q7deuWFBYH3E7yN9X+AFQBUFo9zlj0rqyvvvpqvP3220lxnMabndi7dy9+++03T9dcdtll+PDDDx3j2X0AY8eOxeDBgz3lCygr6Js2bYpKlSp5vjYswqgIUrUgtDNzthoi8zt0RkTYuXNnnjk2Ozfsbif+586d63oY6ODBg4l03TzvOXPmJJzNmuH1nQUxtLlp0yb06tUrEa75k/OTtrH3OmPGjMSxHqMJ8wMPPJCUVt26dfHSSy95yl9Ll4hiMWxth9shssd0f08D+ApA35BkSjvjxo3Ls6Wy9jE9/PDDiYLoh969e6NOnTqerzty5IhjHCdfVffff7/nfO1a2QsWLMjaxYFu8DoHE1R+RITGjRs79pgAZRjOuBbpyJEjpgqxYcOGeSpKO8qWLZvwkeYGJw/FfpSsmzA7xo4dixdeeMHTNW75+OOPHePs378fb7zxhum5dKxrOXz4cEp1mV/cDpH9Q/fXCkAdAOGbLaWRAwcOJI6nTZuWOHayzjJz9aDhtKCvd+/eAJI/Fi8foFXcIIezdu3ahfr14+kZyIsVWRwW+LmVwTgHoR+Tt0rbrIHQu3dvlCxZ0jS+XW9Hn8eePXtsvVanAhE5Dqf973//M73ObfqA9TObMmVKHvdCfnDzDXp11WOXluZT0cs3PnDgQJx55pmu4weFX2eX6wGcHaQg6cbu5eg/PDuT1AULFqBq1aq+K/P+/fu7lm327Nl5FN/69et9uZvQGDt2rEsp4WrLabfPgJk9727olyjXyAQ1RKa/1gmztO32GnIji9Zr9SK/k7zHjh3L48vLrCGmz0ebW9RjVfaNk+N2w4mA0phMxUFquuYl3SixL774Ak8//bRpnKjW8Lmdg3mDiF5X/wYBmArFkizjMK5Ad7IAsvtYzMwn3XD48GFXwyrTpk1LnGvWrBkuu+yyPOfffPNNyzycCr7XuaGg+PDDD3HqqacGmuajjz6KI0eOWD7Hzz77zJNpqRNudmY09nz9zjsQEfbv3w/g+DtzK/+aNWs8NSTGjBmTNLxVrly5PPIEhX5LCT9pW40cPPLII6bhTmtGrLjoootsHeIys+P7GDNmjO8y50Wx9+/f39Rl0YABA9K2Ls2IazNlKObJ86CYJz/BzLeGJlWIzJ07F/PmzUu8MLNxaL1Vh1Xh1z56wHuF5WRRpKXnpAQOHjzoaIhwySWXYPHixXjnnXcsd6xjZixYsADz5s2zlScINm3a5Cm+m0n+l19+2XT/mJkzZ2LOnDno1KkTtmzZ4ilfr3IYadOmTUr5mWH0g+ckU/Xq1W3PG8v2ddddl1ilbtUDc6MMnFb7m8mxZs2aPA02p3tzGlkwXv/11197kkljzpw5mD/fXVvaSuaJEyfmGX0w4nWITH+dm8bw0KFDI7GIBdzPwbyv+/uQmf05xYkBRIQGDRpg0aJFAJSPyi/ay924caPteSfctC6sVmc7KauZM2di8uTJ6NOnD/r165d0/o8//sCsWbNQv3592x6RE1FZqmj5Pv/886YyjBkzBkDe5/fyyy+bDrvY4fZdBv0cwnRc6Wb9jp0cVrjZ3dXYg69evXpijxU32MlTv359LFiwwPSc5v7feL1fAxY37yInJwdt27a1PD9y5Mg8c8Bm6bvJx+/yhjBx2g9mERH9YvWXLiHDwM5SS2/tYTW5rS+gVj0DuzkSPdpGVhqaUzunQtW3b19L5eSm6w4AVatWTShbPe+9957jtcb8okD/HtyaNI8aNcr1otQ///wT69ev9yWPm3Av6WiNCe0+zcxe9ecBmFqfuTVPdpInVYxKwK6XYRxpcJJHq2yN8Xr27GkaX9uGwyv678zvN7B9+3ZMmDDBMR+n8MWLF/vKP0ycejDXAPgHgO/Uv1vUv28B+BvUjBi7gqlN7mtb4/pJh4gwatQobN++3ZX5ohlmVjMaRmVgN2E+c+ZxZwt2hV8/3KdhNZZtRVQ7MepXNdvh9+Nv1aoVKleu7NoKyKni82IMYeT222/P83v48OGO15QtWzYpzK5c5ubmpmVnyDVr1iSG49z00tyYaQPJe87YrU3R8+WXX6JKlSqJ3z179nTdQNSnuWXLloTRgpcyt2LFCtd5WGE3bKjhZEIeNLYSMfNaZl4LoBUzP87Mi9S/JwC0To+I4WBWqCtXrpzyXiEAcOONN3paN+AGrbCff/75ecLtCp5Vt9uIWRpmSscObbFZlJi90+effx6A8wdqZbKrfZB+Wv1BoL8nO5N4QFnb4oalS5danvvpp5/yLEjUyzFz5kzUrFnTVR5O/Prrr57mdYw4zfM8/PDDnia2c3Nz80zmHzhwIPG8nfzWjt8QAAAgAElEQVTBae+/R48eOO2003Drrd6np512cbWag7H7bUa6PXC7neQnImqi+3GJh2tjidkL27Rpk6lPH6vuq13F4nYRW6rYFRi3H+63336bFOa10rSKP2nSJE/ppCqHHw4fPmy5WZr2DK2Go6ziW+F1LsfMOMEqjZUrV3pKW89bb71lex5QdgINCqden5f3bjYPOWPGDHz22Weu03DT+jdDL6fmH00zlQ6y7FrNl+nD//77b8d00j1X6vap3gkgh4jWENFaADkAuoYnVnQ8+OCDSWFXXnllUtjixYttW+12CsaL6agTdk753BYmMwXqZW7DLs6///1vVzIAivuOVD4AzULNbBiAmS3T1pS09izvu+8+XH311QCOPwe7tRJElLD4s8oj1XF6PUE4YgTyyqLf48TMBYnTe/Fq9q7Pw5j28uXLXW8/MWzYMNvhWeN9rFmzxnToa/z48bb5WClXZk7LZntxnMB3g6sdLZl5HoALiOgk9Xd8d95yIKiX72SFpC/YkydPRrNmzRK/vVjLpILWKrv33ntRvnz5UPOyqjjtFvoZSdXVe7169QCYzzvYYaz8P/30U8+TvnXq1HFlieTleVihf9bLly8HEeVZpe2kxMw2iytSpEji2I+CadmyZVIYESUtfNRjtl/JxRdfjNmzZ5umZcadd96JG264wfScWeW/bNkyPP/885g1a5alXPrrtbyt1j4xc1Ljw6k3ZFVOzJYcaL0SKwenZu/6wIEDllsyL1++HOXLl0+bUY6TFdmt6v+HiehhKD2ZO3W/Mw7Nekw/AR4UVq65mzdv7is9N4XAau0KkJ7ucIcOHXDs2DFLWfW9iaNHj1pW3E8//bQnX0m33XabJznN5NMs5dxO4FuhlSk3PuTceEVwylcvb+3atVGrVi3L82Zoiw7r1q1reo2ZgvE7hGT1vvfv349vvvkmkZ+GmXIxxjFiJ5tZb2/Hjh2Ww9764TY3zm4PHjyYNF/pVGZuueUW03AzKzA/IxRNmjQxDQfgySIyCJxKTXH1f0mLv4yjU6dOANybEHvBbnMhJ7PhKVOm5PkdxHi30zCXXUEEgLvuusuxVf7555/j8OHDtve2b98+AMDrr79u6Tn22WefNR3qeP31103jf/LJJ7ZyGTGTT1ucqZ0zmy8yen4wQ+vNvvvuu0mTz4cOHcLu3bsDm9NyOuc3rpOJq9/Gipvr/KatGbHYDUuaVdB2Ckm/YPeSSy5xlGHJkiWWHj3Mnv3ff/9tOVfiNLzq9tzPP/9sGb9QIVeDVoHhZEX2tvq/n9lfekTMDgYPHmxbsPVDaECyO5lUMWtJ2s0hERHeffddTJ061TFt430ZHSN27twZgPmEtRNmvpWA4OYhgOMfqZf5IjO+/fbbhNWaRs+ePXHSSeHvLq6f0NZXOn6GQszKXZgKxm86mjKwm8w3u3+7VrxZGe3UqRMmTpzoRkwAyjq2PXv25JnX0mjUqFFinZsRPwrGK7FSMBpE9CIRlSKiwkT0AxFt04bP8gupbhjl1Rvthg0bXJsZeyEIM2wNreD37t07z0dwwQUX5Im3ceNGTJo0KVCF6fWjMxv2Mab1008/pSSTmems0/oGK+yGC83uQ+uZA0jZLb2fORi7tILAzFpSm8Ox+zbNnpXdXJ9Z/IMHD3re5G3AgAH473//mxTuxweg17JuZ1kaZMPMDW4HVlsz824oCy/XADgDyt4w+YannnoqkHTcfnBXXnml60WEXnDjxkOPXeHWCuvkyZNt482cORMtWrSwPO+ml2SVtxesekN62ffu3eu7UjTKpLnpMebhBisX8qtWrXK891RbvESUVEkRka/nEuY84JAhQ2zPWw2RpQOr+Ti7UQyvPZjXXnvN87s2U3ph4ra/pMW7GsCnzLwrKt9TUeHkysEJNztMGs2XtbkLvxi9LwOw9Qxrhp2nAP3H6+ZD1peZbt26gZkxdOhQNG3qfXPUIIYNjNZjtWvXttw7xY9MbsbwvdKxY8ekxbZBs2jRIlOlb/fMrazF3NQTYVo0eVUwQclidd9+FEzhwoUDkQnw7xLHL257MF8T0VIAFwL4gYjKAgh+/CbGuLEQSpV27dqFnodX7rnnHstz2sc4Z84cV2npP6ChQ4eaLmpNJ1u2bMHq1asDW6MShAmyEytXrkxsOBUWZh6v/TYozfzcGfHq984tfnowQa10D1LB2OG1zMZyiIyZewG4BEADZj4M4G8A8asNQyRbemxeV5LbmSDrx6XdjFHbWdmFjdmHdejQIdSoUSNxf2FuCT1o0KBA0knH+gWzSpCI0KVLF89p9ejRwzGOW793bvbhMeL1eQXl1t7KHU8Y86qxRjOftfsDUAzAvwC8o/4+E8A1bq6N4g8AZ9PfW2+9FVha5513nqt4xYsXd4xToUKFlOVh5YWF/teuXTvLc3/99VfostWrVy8t9xnEX8WKFZPCihYtGrlcXv/69+/P1157beRyuP2bNWuW52sGDx5sGm5VdqtWrZo4b4Z6LrC62O0Q2XAAh6D0YgBgA4DkrdOEUPAzCW5FkEN9XjcPM8PO8WKQWO0PAiAtvYJMwswtSSa2vKOc5PeDHytGqyHsoUOHpipOILhVMDWZ+UUAhwGAmf8GkB1jRhlAkBWgWyecqRoYuOXss89OSz52w19WpsRBuja3U3Bxo2DBglGLEAiHDh1Ky9xpUHhdPGxHt27dAksrFdxakR0iohOhdLNARDUBBL8UXhBCwm77gcaNG5uGe/Vpli0E0TONA1YbAeZn3HimCBJHBUPKbO9gKBuOVSaiDwE0AdAlXNEEQRDyD+kcqmW29i4eJI5DZOrEz2MA/glFqXwMxZosN1TJbCCiq4hoKREtJ6InopJDEAQhKKx8moVBupSZ2yGy+QBqMPM3YQrjBiIqAGAQgBYANgKYQ0RfMnN6ZosFQRBCwM2aoaA4duyYbw/ZXiA3mkxdZHkGgLUA9kGZ4GdmDndJsbksjQD0YeY26u9eqiwv6OKIWZAgCIIFBw8exAknnJAUrvqhC2zszG0PJnlLx+ioCEDv72Q9gIsikkUQBCHjSJf5ttsdLdNreiAIgiCERqwUTMzYAKCK7nclNSwPl19+OX788UcAQOXKlT07eRQEQchWNAWTm5uL3Nzc0PJxNQcTJ4ioIIBlUCb5NwGYDeAmZl6ii8N6M7ymTZti8uTJWeNPTBAEIRV27dqFUqVKJYVHNQcTG5j5KBHdB2ACFDProXrlYoYfd/CCIAjZigyR2cDM3wGo7TZ+OrasFQRByBTSpWDCN4SOAemw9xYEQcgURMEESN26daMWQRBCp23btlGLIGQIomACoGXLltixYweuuOKKqEURAqR8+fJRixBL0mWwM3DgwLTkI4RHiRIl0pJPViuYiRMnonTp0rZxMs2KTgDOPfdcAMA///nPiCWJF+kqyzLknPkUK1YsLflISRFCwWzTKjtuuukm13F37twJABg9erSnPM466yxP8TONdCkYMxcjgmCGKBghJYoXL24aXqiQPwPF1157zfY8M2P+/Pm+0r7gggs8X/Pyyy/7yisK0qVgrPbPCYPChQunLS8heETBhMz9999veW7IkCFplMQ/L7zwguW5evXqWZ6rXLmy57yqVq3qKf51113nOu6aNWs8SgM8/PDDeX63b98ev/zyi+d09FSrVi2l663IxiGyJ598Mm15CcEjCiZk7Laf7dy5cxol8U+NGjUslYVZpeZnfFdbVWznbcHMpcXnn3/uOg8/e54bKV68OM4777yU0gjLo4SdgildunRGzjdmwpbHQVrvZct21RqiYACcf775rgNdu3ZNOe06depYnsuUwmSsmPTKxqzS+uijjzzn0b9/fwBAixYtPF/rluHDh6echpvWu1OcdFT0J598cp7fZ5xxRuh5hsHBg/HfmV0zOgmCkiVLBpaWFXPnzg09Dw1RMADuvPNO0/BKlSrl+f3OO+94TvvUU081DX/ppZdCHWoI2gxRXynqj83s6du1awfAW0td6/WUKFEC06dPR58+fWxl0DNnzpzEsd3H3qVLF09DamY0b97cMc5TTz1lez6sHow+3d69eyeOK1Wq5Nkgwo509oQOHDiQ0vVlypSxPFejRo2U0tYoVqxYUl3hlyeeMN+gt1+/fr7nNY1ceOGFgaTjBlEwOP7BPP7443nCjRXBXXfd5Tltq8rEanI8KJ5++ulQ0h00aBC++uqrxG8v1l926D+eSy65BH379k2KY1WxNWjQwHU+qbgNWrBggWOv1kr2dKM9z9zcXKxatQpVqlRxuMI96d47PhXatm2Lyy+/3PTcwoULXc+DnnPOObbnv/76a9vz11xzjat87r77btNwIsJll13mKo04IQoGykTxAw88kJjMfuihhwCYK4f27dunVba40bNnzzyeEewWsWrPb9y4cY7pFixY0LYyKViwIOrXr295fuLEiXnyNFK7tmvXdUm0atXKNm3geAU0ffp010Nku3fvxu7duz3L06lTJ1PvFHrrOk2G+vXrB26JFfVczpdffuk67uuvv255zstwVMeOHZPCPvjgAwBKuahVq5bpdVpDr1y5cnnCjb81sm2NUXbdjU+qVKliujrZrEL54osvbNPSV4J2QyX6tN12fa+99lpX8cLAWKloY/xuKpsKFSqknP/gwYNT6n04DUt1797dsjKaMGGC5XV///03AG/Wb9ozK1mypK8x9xIlSpgq9r/++itxXLBgQZQvXx4nnnii5/SdSNXNiBfLMGP5YmZP30FQjm7NlPQtt9zieJ1myGO8DysFk21biuRLBePUoiYi1KhRA82aNUs5LyszXn1BGjVqlOVcjR6r8dkOHTr4E84DZh+6Hqt5LCC1se5LL73UND8jRYsWBWDfowKAe++91/Kc0aR8x44djvJpFXjx4sXRo0cPx/hAelzdEBE2b97s2HjRWswVK1Z0nbaXHoyZ0YbT/esNY1JRZkH2BoxbfugVxNGjRy2vM1MYp59+umWjxUpmP4rHbL+XdJMvFUybNm3w22+/mZ6rWLEimjVrhpUrV9oqmIsvvtg0XP/xERGqVKliauKrr3SrVauWmJOxsjo799xzcckll+Af//hH0rlGjRpZyhkEZmO/xkrGbDK1ffv2aNasWaKgt2nTxneL2slooUmTJgCcraUaNmxoGm5WaRrdDNl95AUKFEBOTo5t3hoPP/wwNm3a5CquX9xWrlo8faXvZMigf1aHDh1KHN9zzz1Jcbt3754Udsopp+T5bexl6Of13CgzK+VoV/FruDWDLlCgQJ4KW28scvToUUs5q1evnmfOEgCWLFli2asPsgej1TFROkHNlwoGAE477TTT8PXr1yesoOxw0/W2KywtW7YEAOzfvx/169d3/JCmTZsGABg2bFjSObPKxC5v4zySNudkhdWzAuwrgNdeey2xduXVV1/FiBEj8N1339nmZYSIsGzZMtxwww2JsPvuu880nv6/1XknzBS4mzSc3t/pp5+eOC5UqJDtMzViVK5GOfTPRsNOkdv1Ntu0aYNhw4ahU6dOlnHKli2bONYPHb311luW13hBf39uFMzAgQPx2WefuUrPiNUSBbM0zHpTb775pqli1V9nnOC361mEMQdzyy23oEyZMrj++usDT9uJfKtgTjnlFN+Tlb///js+/vhj0+uNPRhjmBFtaMdJFi2tMmXKJBVQs/U0dukZ55FuueUWxzU5xklOpyEzIw899BDKlCnjyXrujjvuQPfu3VGrVq08H57XYZOhQ4fijTfeSPzWt6CNvVR95ekFL2XJSyXSqVMnPPfcc0nhen9gI0eOTDpvVamOGzcuMYy4evXqJLlbtWqF6tWr225xUblyZQwePNiV/A0bNnRcmGqUVd+7d/Ncich3xXzJJZfk+a33smB0iaMvd5rM9957LypXrhyY4YNdA8ksjxdeeCFP2TajZMmS2LZtWyDyeSXfKRhjgXKDcUKuZs2aSd18M4yFRZtnefHFFx3juuX666/3ZKZrRr169XD48OE8wx1GjLbzWmH3+mF5iT9s2DDTiVQ7BWOsaDp27IiuXbvmmQvQ9x40hZJqBeFW6ZUuXdq2l2Tkk08+wQMPPJAnrHXr1o7yWp1v06ZNoqxVq1YtEa9hw4Y47bTTEpWql/TtJu2rVavm2bWONu/mRo6CBQvi3HPPDWR/k2nTpmHp0qUJww29orOq4DWc5rqcGqP6fKwwU6IVK1a0XNeipR+l4UC+UzBO9ux+WLt2renCQA3tRRcpUiTPb7M4APD2228nnbcqJCNGjEhate0G/RAfEYGIbM1ZjTKfdNJJKFiwoO1CNjfp+MEujTAn0O0+VCszVSPNmjVzbTU4derUpLBFixbhhhtuCKzFXKFCBZxyyinIycnBpk2bEvN5XtL3K4tWBo3P1WyI7MEHHzRN48iRI6hdu7arNSJOw9pNmjRBkSJFcOKJJ2LXrl147LHH8pw3G53QKFKkSCCT6lZl7KyzzsLQoUOTwrWhdiPr1q2L3JwcyIcKJgzKlSvnajhMW+xmFkdfULTFVk7jw+vWrUORIkV8tVDsPhan+AAwe/ZsrFq1ChUrVgQzJ8778XbglSuvvNL02Wzbts3VHjHvvvuur3ytntO+ffvw7LPP+krTirfffjtPSz4s5s+fjyVLliS1jrX36UZh+63IvJQ7J+WgjTLY9ebXrVtn6lH71ltvTQorVapUnmFj4xyMmezjx49PrMcyot3H0qVLLeWzo0OHDknzdtdff73l+9F7FnCanwyTfKdg/HwMdi1OZk7MoxgxvtDx48eje/fuppWgmeWOvvDo02rSpAnq1q2bKERmllPGvFu0aJGnErF6DgMGDDANN1KhQgXT1eFO3g6CaFVdd911+Pnnn5PCy5Qp4+ojaty4McaNG4fLL78cgwYNcpXn9ddfj+rVq5ueK1asWNr8yrlpyDhx4YUXJiwby5UrZ7kmwy6fIN6j9l25MZ4w5me14t2O4sWLp+STzWkYrlGjRmjZsqVpz1NDb+xhpEGDBrZ77Xg1YJEeTAawY8cOz4vhevbsCSC5MihVqhQGDx6MM888M+mas88+OylMX3D0x1999VUeh3WFCxc2VVAvvfRS4vjOO+/E/v37E367rHowTzzxRGLluh6/Y/5xpU2bNpg0aVJCiTvJP2rUqJR2Afz4448BxGMhXa1atfDHH3/Yxmnbtq1rx6N+3v20adOSrK+uvvpqAO6syKzWhAHJ3rydWvBWeRi/P6cejIZdz9Puuv/7v/8zPW+1DMGpLMXhmxQF40Dp0qU9deWZ2XWrWM9FF12U5/c777yDAQMGmLqoL1iwYFKL+ccff0xq1T3yyCN5fp9wwgmJIYQwJ+fTkY4bLr30Ulcm54B/6zEnNNNQbaFeqvfvtgeTaj4XXnghvv/+e1dx/eTVpEmTRM9fuyezdSt+JqrtFIYXWY0KRn+tG99uxvUvxjTdMnPmTNdxzeZjtTyvvPJK13OFQRGMe07BFO3Fvvjii9i4caNj/NGjRyeGvbShplQtQZgZRJTUQ7L70MzOOW2T26FDB6xcudKXjGFhN1ShZ/369Z6NFdwSlkKNQ+tUIyilaaZggrAOs/IEbszfiHFe6osvvkD79u1Rq1atpMabkUOHDkWyG+epp56a8EBhNC664447cMcdd6RVHlEwAWFXcN34LAJgOjfjR7H069fPcTW010rhxhtvtN3TvlGjRq42/4pTxaihXwnuVr7GjRu7mrAN+n7t3MKfcMIJtqbmfnC7DkVP0aJFsW/fPixZssRV2kOGDMH27dsT+wjp02vevHnCqaRXucziar2mX3/91fE6431pvWHNgtIOo3Jx01D0+q2bxb/mmmswe/bsPHk6uU8Kk3w3RObng/fbPQ9yvN2L3P/+978d96fw+hwKFSqU8nqbdJGOCfcJEyZg1apVoaVv9n7ee++9hClsEI0GP2heA/R59e7dGz/88AMAxdlro0aNUKBAAdcbcbVv3x5du3Y17cHYeR3ww5tvvolffvkl0E3CgsJvfaG/7tVXX8X06dMDSTcIpAcTEh9++CFat24dWHpBVx5BDD34IR2VYKFChVz5oTLD7cfodkO3l19+Oc+wRFAfu3HODgjHzYjxfZl5DShdunSilTxr1izX92hM+4wzzjCdczSjSpUqeTZVs0tXz8knn+x63Viqbofccv7553tejKrH6DNP4/nnn8fy5ct9pxsE+U7BpGtM/Oabbw41/bilFydSubegn0uVKlXyTAj7TT8nJyePD7kbbrgBnTt3TmwpPG7cOBQuXNjUAjCdeJl3MG5xkJOTg1deeQX79u1zlY+ZCx0znJ75TTfdZJpnukYgtPk/fX7Vq1fH5s2bsX//fse0a9eujV27dmHPnj15wt0auIRJvhsi80McTEvd4EVOux6MmUPNoEiHYss05blgwQLHOD169LBdbNimTRvLVd2poJ//6NKlS6Bp33bbbXkqxcKFC6NkyZJJ5TgnJ8dyrZkZX3/9tatN7jTatm1rOn9oZS7tZe8fuzQBpcc3YsSIpHhly5ZNGuq1K9elSpXytOVCuhAFExB+fXO5JWglZyen2fYCQjCYvUe9Y8mgPAIEVQ6vvPLKQNKxMmBxM9TYo0cPT+W/RYsWaNOmTeK332dhZlm4fft2vP/++77SM+Piiy+2VAxODmUzoeErQ2QucPMi7VZDp8rChQt97Xxoh/Yc7HbdDIN09C4KFy6Mw4cPh55PGLid2wGi7al5zXv06NEpVYhuNuSzwu9z0s9taLL78funJxUrskzrmQP5sAcTVuvcbqfEVDHzn5QqWmENU+6o+Pnnn7Fs2bKU0rjqqqsCksYbqTqZNBtuyQa6du2KtWvXOsYzWzDrZYsIK4Kq3INYMOonrajIVz2Y7du3B94T0EiXL6qgCWPPdjvMrJ+CVvqp+JvS+OSTTwKQJBm3e8TXr1/fV/ph9qT1WPllc8Pll1/u+ZqCBQs6rp5fv3696fcdh4rYy3qinJwcnH322QnXOcY09u7d66mnGyX5SsGk2r11Q6Z1Y9PxTPSYWRkNGTIkrTJExcqVK127pLHa2llPusuavqJu27atr0Wd99xzj+MaLb/EcZLbiLaq3gzt+fbo0QNA8vvVGrHFixdHt27dcOONN4YkZXDkKwXjlzi0gILGzS6W6SKsXmUqhPFs9BtYBYFXdz9BE5YrlGz83oJ4H0WLFsWuXbsAZE6jTBSMC7KxwJu534iC+fPn227PGxWpeE4Wso+4jEwEsalZOonlJD8R9SGi9UQ0X/27SnfuSSJaQURLiCi4pfKCKW73h/FLvXr1slKBe0WrwOrUqROxJOESxbtOxQItKPJrGY+lglF5hZnrq3/fAQARnQ3gegBnA2gDIIdi9uaibOmE8Sjatm0beJqCNXXq1HFdhhYtWoTffvstZIkyn5ycHFcWaGHix2loNhDnITKzp90OwEhmPgJgDRGtAHARgJ/SKllM6dmzZ0ZMdArWeKlkzDap0wijoROEr7MoKtESJUpkhNWV8dnEZVguFeKsYHoS0W0A5gJ4hJl3AagIQL/7zgY1LFQypWVRp06drB9iEaJDvzreL9lQafohv953ZAqGiCYCKK8PAsAAngKQA+AZZmYiehbAywC6eUm/b9++iePmzZubbinsQVZX8WbOnBnaxlVC9vPf//435fKjWb8F3SiqXLkyOnTokFIakyZNwvnnn+8YL1MadEETRQ8mNzfXtQdrP0SmYJjZrdvXIQC0vUc3ANCvyqukhiWhVzDpwmrvbEFww6OPPppyGprb/Di2mP0ssMwvPPnkkyk1gv1ibHz369cv0PRjOclPRKfpfv4TgLb93FgANxLRCURUHcAZAGanW778RH5tTWYqYewJk2405Zju7X3NGDhwIIDUFfZzzz2HTz/91PL8888/n7ToOY6NBK/EdQ7mRSKqC+AYgDUAugMAMy8mok8ALAZwGMC9nA1vQRBijjQ0UqNatWqoVq1a1GKknVgqGGa+3eZcfwD90yhOVrQKBSFTEGWmUL9+fezduzdqMVIilgombowZMwbbtm2LWgxBcE3QHXup9NNPbm5uxg+TiYJxgXHbWyEYPvzww6hFSCLTP2ghWKJUrHaOMTMFGfsRbAnLIWbdunVx8803h5K2EDzprGjj1FuSBkdqiIIRbDnrrLOiFiGtxKlyy69olbq8i8xHhsgEW/LbRy4t1mT+9a9/5UsLKOB4+Zdy4Q9RMIKQhQRZIf7nP/8JLC035LdGTTYjCkZIO927d0eDBg2iFkMQHBFllxoyByOkncGDB6NbN0+u5TKOZ555JtL8pWIU4oAoGEEIgaeffjrScftsmDMQJZn5iIIRBEEQQkEUjCBkGSeffDIuvPDCqMXICrKhJxglMskvCFnG9u3boxYh6xBF4w/pwQiCECviNPcSJ1kyEVEwgiVRbIAkCNJbyB5EwQiW3HrrrVGLkHZKliwZtQiCkDWIghEsyY/DA/3798eSJUuiFiNfo5W7OJQ/6U2lhigYwZL8+HEVL1483zn4FISwEAUjCEIsuffee6MWIUF+bGwFgSgYQRBihTY0Vq9evYglEVJFFIwgCLFCegvZgygYQRAEIRREwQiWFC1aNGoRhHxIHKzHjIj5uj/EVYxgysKFC1GnTp2oxRCEWFC3bt2oRchIRMEIplxwwQVRiyAIQoYjQ2SCIAhCKIiCEQRBEEJBFIwgCIIQCqJgBEEQhFAQBSMIQqyIo5my4A9RMIIgCBYUKiSGtqlA2eiWgYg4G+9LEPIDBw4cwIknnhgLlzEHDx7E7Nmzcdlll0UtSlogIjBzYF1IUTCCIAgCgOAVjAyRCYIgCKEgCkYQBEEIBVEwgiAIQihEpmCIqCMR/UpER4movuHck0S0goiWEFFrXfhVRLSUiJYT0RPpl1oQBEFwS5Q9mEUArgMwWR9IRGcDuB7A2QDaAMghhQIABhacqYMAAAqsSURBVAG4EsC5AG4ioozaPD03NzdqEZIQmdwhMrknjnKJTNEQmYJh5mXMvAKA0WKhHYCRzHyEmdcAWAHgIvVvBTOvZebDAEaqcTOGOBYokckdIpN74iiXyBQNcZyDqQhgne73BjXMGL5eDRMEQRBiSKjLVIloIoDy+iAADOApZv4qzLwFQRCEaIl8oSUR/QjgEWaer/7uBYCZ+QX193cA+kBRTn2Z+SqzeIY0ZZWlIAiCD4JcaBkXRzv6GxoL4EMiehXKENgZAGZDGc47g4iqAtgE4EYAN5klFuQDEgRBEPwRmYIhovYA3gBQBsDXRLSQmdsw82Ii+gTAYgCHAdyr+n05SkT3AZgARdkMZeYlUckvCIIg2BP5EJkgCIKQncTRiiwl0rkYk4iGEtEWIvpFF3YyEU0gomVENJ6ITtKde11dQLqQiOrqwjur8i4jottTlKkSEU0iot+IaBERPRC1XERUhIh+IqIFqkx91PBqRDRLzeNjIiqkhp9ARCNVmWYSURVdWqaLcFOQrQARzSeisTGSaQ0R/aw+r9lqWNTl6iQi+lS9x9+I6OKIy1Qt9fnMV//vIqIHYvCc/o+UBeS/ENGHarmJtEwR0YPqd5f++oCZs+YPisL8HUBVAIUBLARwVoj5XQqgLoBfdGEvAHhcPX4CwAD1uA2Ab9TjiwHMUo9PBrASwEkASmvHKch0GoC66nEJAMsAnBUDuYqp/wsCmKXmNQpAJzX8LQDd1eMeAHLU4xugrIsCgHMALIAytFtNfdeU4jv8PwAfABir/o6DTKsAnGwIi/r9vQfgDvW4kJpupDLpZCsAYCOAylHKBOB09d2doCtLnaMsU1AWpf8CoAiUb28CgJrpek4pvdi4/QFoBOBb3e9eAJ4IOc+qyKtglgIorx6fBmCJejwYwA26eEugmHDfCOAtXfhb+ngByDcGQMu4yAWgGIC5UBbObgVQwPjuAHwH4GL1uCCArWbvE8C3WjyfslQCMBFAcxxXMNuilElNYzWAUw1hkb0/AKUArDQJj0uZag1gatQyQVEwa6FUxoWgGCy1irKcA+gIYIju978APKbdf9jPKduGyOKwGLMcM28BAGbejOPrgKxks1pYmjJEVA1KD2sWlMIUmVzqUNQCAJuhVOorAexk5mOGfPPIxMxHAewiolOClgnAq1A+NlZlPBXAjohlgirPeCKaQ0Td1LAo3191AH8S0XB1SOodIioWsUx6bgDwkXocmUzMvBHAywD+UNPZBWA+oi3nvwK4TB0SKwagLZSeXlqeU7YpmDhiZUURqik1EZUA8BmAB5l5r4kcaZWLmY8xcz0ovYaLoAzbuSVwmYjoagBbmHmhIX23eYX5/powcwMolUFPIroM0b6/QgDqA3iTmesD2AellR1pmQIAIioM4FoAn1rIkDaZiKg0FPdVVaH0ZooDuMpLEkHLxMxLoQyHTQQwDsrQ21GzqGHIlG0KZgOAKrrfldSwdLKFiMoDABGdBqV7rMlW2US2wGVWJxE/A/A/Zv4yLnIBADPvBpALoDGA0qQ4MTWmn5CJiAoCKMXM221k9UMTANcS0SoAHwO4AsBAACdFKBMAgJk3qf+3QRnivAjRvr/1ANYx81z192goCicOZaoNgHnM/Kf6O0qZWgJYxczb1R7JF1DKWZTlHMw8nJkbMHNzADuhzMum5Tllm4KZA3UxJhGdAGXccGzIeRKSF4p2UY+7APhSF347ABBRIyjd5i0AxgNopVrpnAxlzHZ8ijINA7CYmQfGQS4iKqNZqRDRiWpaiwH8CKCTGq2zQabO6nEnAJN04Teq1jfVcXwRrmeYuTczV2HmGlDKySRmvjVKmQCAiIqpvU8QUXEo8wuLEOH7U9NbR0S11KAWAH6LUiYdN0FpIGhEKdMfABoRUVEiIhx/TlGXqbLq/ypQPNh/hHQ9p1Qn2OL2B6VLugyKF+ZeIef1ERTrlYNQCtcdUCb4vldlmACgtC7+ICgWIT8DqK8L76LKuxzA7SnK1ARKF3ghlO7wfPWZnBKVXADOU+VYCMWi5Sk1vDqAn9T0RwEorIYXAfCJmvcsANV0aT2pyroEQOuA3mMzHJ/kj1QmNX/t3S3SynCU709N6wIoDbiFAD6HYk0UtUzFoBhllNSFRS1TH7Uc/ALgfSjWrFGXqSlQ5mIWAGiezuckCy0FQRCEUMi2ITJBEAQhJoiCEQRBEEJBFIwgCIIQCqJgBEEQhFAQBSMIgiCEgigYQRAEIRREwQgZibqYdlGc01bdm78RhEwpytGOiM7S/e5HRFdEKZOQPxAFI2QygS/iUl12BJl2oDLq5PNCeyhu2wEAzNyHmSfZxBeEQBAFI2QyhVTPvr8S0XdEVAQAiOgCdQOnhUQ0Wuem5kciqq8en0pEq9XjzkT0JRH9AGV1cwIimkxE5+t+TyWi8wxxihDRMFI2mZpHRM11p6uo+S4jon+r8YsR0dekbJT1CxF1UsPrE1Gu6kX5W52vqB+J6FVSNh97iojW6PIuRkR/EFFBIupGRLPVdD9VXZY0huIM8kVSPCFXJ8Ur8j/V61uo4T8T0bukOI8EEa0mor7q/fyscxMjCK4RBSNkMmcCeIOZ60Bxjd5BDR8B4DFmrgvFRUYfi+v1vYt6AP7JzJcb4gyF4gIIRHQmgCLMbBw+6wngGDOfD+BmAO+rvvAAoCEU/08XAOikKrirAGxg5nrqNd+R4qD0DQAdmLkhgOEAntflUZiZL2LmZwAsIKJmavg1AL5jxbniaDVOPSj7otzJzDOh+Jd6jJnrM/NqLUFVIQ+HshnWBVDcmvTQ5bmVmS+EskfIYxbPUBAsEQUjZDKrdJX9PADViKgUlJ32pqnh7wNo6iKticy8yyT8UwBXq0NTXaHs7GjkUii7YoKZlwFYA0Br8U9k5p3MfACKD69LofgYa0VE/YnoUmbeA6A2gDoAJpKyb85TUFy+a4zSHX8CZQ8UQHHUqZ07n4imkLKF983QDYtZUBvKM1yp/jY+qy/U//OguKAXBE8UiloAQUiBg7rjowCKqsdWe1gcwfFGVVHDuX1mFzDzfiKaCGUeoxOAC13Ipc8/aX8SZl6h9mTaAviPOjQ3BsCvzNzEIk29fGMBPKd6tb0Qx73wDgdwLTP/SkSdoTjx9CKrEe35HoXUFYIPpAcjZDJJlSMr+81sJyKtor4NwGT1eA2ABupxJ9ijT3sogNcBzLbo5UwFcAsAqHMVlaF4qQWUnkppUrYpaA9gOhFVALCfmT8C8BKUvVWWASirukgHERUionPMBGPmfVC2nR4I4Cs+7rG2BIDN6jzKLbpL9kDZ9tjIMgBViaiG+vs2KHv1CEIgSKtEyGSsLLS6ABisVuqroM6hQKnMPyGiuwB84zZtZp5PRLuh9BDMyAHwljo0dRhAZ2Y+rGwJgtlQhsYqQtkAbj4RtQbwXyI6BuAQgB5q/I4A3lCNEgoCeA3Kvjlm9zkKylCZvpfytJrfViju4Uuq4SMBDCGi+6Hs0c7qfR0kojsAfKYOAc4B8Lbx/gXBL+KuXxAcIKLToWxI5mWbZ0HI98gQmSDYQES3AZgJoHfUsghCpiE9GEEQBCEUpAcjCIIghIIoGEEQBCEURMEIgiAIoSAKRhAEQQgFUTCCIAhCKIiCEQRBEELh/wH1SDJ5Eg000gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f61213f9dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xsort, res, color='k')\n",
    "plt.ylabel(\"residual\")\n",
    "plt.xlabel(\"hourly observation\")\n",
    "plt.title(\"Residuals (autoregressive model)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(using Goodreads data from examples above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "542338"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(interactions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactionsTrain = interactionsWithPrevious[:500000]\n",
    "interactionsTest = interactionsWithPrevious[500000:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FPMC implementation we built above allowed us to control which terms (user/item or item/item) were included"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelFPMC = FPMC(5, 0.001, 1, 1)\n",
    "modelMF = FPMC(5, 0.001, 1, 0)\n",
    "modelFMC = FPMC(5, 0.001, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactionsTestPerUser = defaultdict(set)\n",
    "itemSet = set()\n",
    "for _,u,i,j,_ in interactionsTest:\n",
    "    interactionsTestPerUser[u].add((i,j))\n",
    "    itemSet.add(i)\n",
    "    itemSet.add(j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AUCu(model, u, N):\n",
    "    win = 0\n",
    "    if N > len(interactionsTestPerUser[u]):\n",
    "        N = len(interactionsTestPerUser[u])\n",
    "    positive = random.sample(interactionsTestPerUser[u],N)\n",
    "    negative = random.sample(itemSet,N)\n",
    "    for (i,j),k in zip(positive,negative):\n",
    "        sp = model.predict(userIDs[u], itemIDs[i], itemIDs[j]).numpy()\n",
    "        sn = model.predict(userIDs[u], itemIDs[k], itemIDs[j]).numpy()\n",
    "        if sp > sn:\n",
    "            win += 1\n",
    "    return win/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AUC(model):\n",
    "    av = []\n",
    "    for u in interactionsTestPerUser:\n",
    "        av.append(AUCu(model, u, 10))\n",
    "    return sum(av) / len(av)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration 10, objective = 1.1740164\n",
      "iteration 20, objective = 0.9275751\n",
      "iteration 30, objective = 0.74895376\n",
      "iteration 40, objective = 0.7000136\n",
      "iteration 50, objective = 0.6867862\n",
      "iteration 60, objective = 0.68066156\n",
      "iteration 70, objective = 0.67874753\n",
      "iteration 80, objective = 0.67833763\n",
      "iteration 90, objective = 0.67842907\n",
      "iteration 100, objective = 0.6778376\n",
      "FPMC AUC = 0.7603135912093875\n",
      "iteration 10, objective = 1.2148852\n",
      "iteration 20, objective = 0.926385\n",
      "iteration 30, objective = 0.87147164\n",
      "iteration 40, objective = 0.7567261\n",
      "iteration 50, objective = 0.6929133\n",
      "iteration 60, objective = 0.68015\n",
      "iteration 70, objective = 0.67653024\n",
      "iteration 80, objective = 0.67463166\n",
      "iteration 90, objective = 0.6742325\n",
      "iteration 100, objective = 0.6739234\n",
      "MF AUC = 0.749866757590709\n",
      "iteration 10, objective = 1.2447531\n",
      "iteration 20, objective = 1.1213548\n",
      "iteration 30, objective = 1.0458206\n",
      "iteration 40, objective = 0.76125264\n",
      "iteration 50, objective = 0.7084387\n",
      "iteration 60, objective = 0.69079274\n",
      "iteration 70, objective = 0.68107307\n",
      "iteration 80, objective = 0.67796034\n",
      "iteration 90, objective = 0.67755425\n",
      "iteration 100, objective = 0.6774621\n",
      "FMC AUC = 0.7502356931668346\n"
     ]
    }
   ],
   "source": [
    "for model,name in [(modelFPMC,\"FPMC\"), (modelMF,\"MF\"), (modelFMC,\"FMC\")]:\n",
    "    for i in range(100):\n",
    "        obj = trainingStep(model, interactionsTrain)\n",
    "        if (i % 10 == 9): print(\"iteration \" + str(i+1) + \", objective = \" + str(obj))\n",
    "    print(name + \" AUC = \" + str(AUC(model)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FISM implementaion, using a factorization machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "nItems = len(itemIDs)\n",
    "nUsers = len(userIDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "fismInter = random.sample(interactions,100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Factorization machine design matrix. Note that we have two sets of features (the user history, and the target item). Both are of dimension nItems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = scipy.sparse.lil_matrix((len(fismInter), 2*nItems))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(len(fismInter)):\n",
    "    _,u,i,r = fismInter[n]\n",
    "    item = itemIDs[i]\n",
    "    history = itemsPerUser[u]\n",
    "    for j in history:\n",
    "        if i == j: continue # Exclude the target item from the history\n",
    "        X[n,itemIDs[j]] = 1.0 / (len(history) - 1) # One-hot encoding, normalized by history length\n",
    "    X[n,nItems + item] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = numpy.array([r for _,_,_,r in fismInter])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fairly slow and memory-hungry (every row contains a copy of a user's history). Could possibly be implemented faster in Tensorflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "fm = als.FMRegression(n_iter=1000, init_stdev=0.1, rank=5, l2_reg_w=0.1, l2_reg_V=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train,y_train = X[:90000],y[:90000]\n",
    "X_test,y_test = X[90000:],y[90000:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FMRegression(init_stdev=0.1, l2_reg=0, l2_reg_V=0.5, l2_reg_w=0.1, n_iter=1000,\n",
       "             random_state=123, rank=5)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fm.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = fm.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MSE(predictions, labels):\n",
    "    differences = [(x-y)**2 for x,y in zip(predictions,labels)]\n",
    "    return sum(differences) / len(differences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6111794992027808"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSE(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(still using Goodreads data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactions.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactionsPerItem = defaultdict(list)\n",
    "for t,u,i,r in interactions:\n",
    "    interactionsPerItem[i].append((t,u,r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.shuffle(interactions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression-based approach. Just collect past K interactions (ratings) as features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feat(t,u,i):\n",
    "    older = [r for (q,u,r) in interactionsPerItem[i] if q < t] # Collect previous ratings\n",
    "    f = []\n",
    "    for k in range(1,6):\n",
    "        try:\n",
    "            f += [0,older[-k]] # Previous rating\n",
    "        except Exception as e:\n",
    "            f += [1,0] # Or missing value indicator if we don't have history of length k\n",
    "    if len(older):\n",
    "        f += [0,sum(older)/len(older)] # Add feature for the average (going beyond just last k)\n",
    "    else:\n",
    "        f += [1,0] # Missing value indicator if no interaction history\n",
    "    return f + [1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [feat(t,u,i) for t,u,i,_ in interactions]\n",
    "y = [r for _,_,_,r in interactions]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test = X[:400000],X[400000:]\n",
    "y_train, Y_test = y[:400000],y[400000:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.65119630e+10,  2.30226517e-02,  1.77267075e-01,  4.43439484e-02,\n",
       "        1.67497635e-01,  4.66625690e-02,  1.63910866e-01,  4.96263504e-02,\n",
       "        1.32824421e-01,  4.92523909e-02,  1.65119630e+10,  2.88013697e-01,\n",
       "        1.90127087e+00])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = sklearn.linear_model.LinearRegression(fit_intercept=False)\n",
    "model.fit(X_train, y_train)\n",
    "theta = model.coef_\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = model.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3810078092245268"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSE(y_pred, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Factorization machine-based approach. Copy the same features from the model above, but also include a user term. In theory, this should allow us to learn how sensitive a particular user is to herding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "XF = scipy.sparse.lil_matrix((len(interactions), nUsers + len(X[0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(len(interactions)):\n",
    "    _,u,i,r = interactions[n]\n",
    "    user = userIDs[u]\n",
    "    XF[n,user] = 1\n",
    "    for j in range(len(X[n])): # Copy features from previous model\n",
    "        XF[n,nUsers+j] = X[n][j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "fm = als.FMRegression(n_iter=1000, init_stdev=0.1, rank=5, l2_reg_w=0.1, l2_reg_V=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = numpy.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train,y_train = XF[:400000],y[:400000]\n",
    "X_test,y_test = XF[400000:],y[400000:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FMRegression(init_stdev=0.1, l2_reg=0, l2_reg_V=0.5, l2_reg_w=0.1, n_iter=1000,\n",
       "             random_state=123, rank=5)"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fm.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = fm.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5863953730808005"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MSE(y_pred, y_test)"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}