{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pymc\n",
      "import daft\n",
      "import csv\n",
      "from scipy.sparse import coo_matrix\n",
      "from sklearn.decomposition import ProjectedGradientNMF\n",
      "from itertools import product, izip"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Welcome"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "**Bayesian networks are a method for building your own machine learning models.**\n",
      "\n",
      "- With Bayesian networks, we model our data problem as a **causal network**, or a **story** involving **hidden** and **observed** variables.\n",
      "   - We come up with **a story** that we think **explains our data**.\n",
      "   - We **use Bayes's rule to find** _posterior_ probability distributions of the **hidden variables** given our observed variables (data).\n",
      "   - We use the full posterior distributions for our next action (decision, recommendation, prediction).\n",
      "- **Advantages**:\n",
      "   - **Confidence estimates.**\n",
      "   - **Flexibility.** We can change the story and re-train the model.\n",
      "     - Classical machine learning methods are **inflexible**: Code and theory only works for its specific problem.\n",
      "     - Your real world data probably has some important \n",
      "Difficult to adapt to your particular problem, which may not have been studied by researchers.\n",
      " - **Disadvantages**:\n",
      "     - **Slow** (unless you do a lot of math)\n",
      "     - **You have to think.**\n",
      "- **For many of your data problems, you may want or need to fit a custom machine learning model.**\n",
      "     - Some important aspect of your data probably hasn't been studied by researchers.\n",
      "     - I don't agree with this picture: http://scikit-learn.org/stable/_static/ml_map.png\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Bayesian Statistics"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have some data $X$. How did it get there?\n",
      "\n",
      "**We don't know.** There is no one right answer, only answers that are more or less supported by the evidence.\n",
      "\n",
      "A textbook example: Somebody tests HIV positive.\n",
      "\n",
      "Story:\n",
      "\n",
      "- Grab a person off the street in New York. (0.03% of NY is diagnosed with HIV)\n",
      "- Give them an HIV test.\n",
      "  - For HIV positive patients, test returns positive 99% of the time. (True positive rate.)\n",
      "  - For HIV negative patients, test returns negative 0.1% of the time. (False positive rate.)\n",
      "\n",
      "Does the person have HIV? **We don't know**. But we can compute probabilities, given this story, using Bayes' rule.\n",
      "\n",
      "Let $+$ denote the event the 'test positive', $-$ denote the event 'test negative', and $H$ denote the event 'HIV positive'. Then using Bayes's rule,\n",
      "$$\\begin{eqnarray*}\n",
      "p(H|+) & = & \\frac{p(+|H)p(H)}{p(+)} \\\\\\\\\n",
      " & = & \\frac{p(+|H)P(H)}{p(+|H)p(H)+p(+|\\overline{H})p(\\overline{H})} \\\\\\\\\n",
      " & = & \\frac{(0.99)(0.0003)}{(0.99)(0.0003)+(0.001)(0.9997)} \\\\\\\\\n",
      " & \\approx & 0.229\n",
      "\\end{eqnarray*}$$\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Bayesian Statistics, More Formally"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To do Bayesian statistics, we need\n",
      "\n",
      " - A **probability model**, a story of how the data came to be.\n",
      " - **Hidden variables** $\\theta$, about which which we want to infer probabilities.\n",
      " - **Observed variables** $X$.\n",
      " - A **prior** on our hidden variables $p(\\theta)$: What do we believe before seeing any data?\n",
      " - A **likelihood** $p(X | \\theta)$ relating hidden variables to observed variables\n",
      "\n",
      "Then Bayes' rule ties them together:\n",
      "\n",
      "$$\\begin{eqnarray*}\n",
      "p(\\theta|X) & = & \\frac{p(X|\\theta)p(\\theta)}{p(X)} \\\\\\\\\n",
      "& = & \\frac{p(X|\\theta)P(\\theta)}{\\int p(X|\\theta)p(\\theta) d\\theta} \\\\\\\\\n",
      "\\end{eqnarray*}$$\n",
      "\n",
      "Bayesian inference is\n",
      "\n",
      " - Useful: hidden variables\n",
      " - Hard: the denominator is high-dimensional integral. You can't just use Simpson's rule; it will take exponential time.\n",
      "  - We use Markov Chain Monte Carlo (MCMC) with the PyMC package."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Movie Recommendations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have $N$ users and $M$ items. Given a small samples of ratings $r_{ij}$, predict $\\hat{r}_{i'j'}$ for pairs $(i', j')$ we don't observe.\n",
      "\n",
      "Our toy dataset: http://mlcomp.org/datasets/736, 190 training ratings (81 test), $N = 6, M = 91$.\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def loadData(filename):\n",
      "    data = np.loadtxt(filename)\n",
      "    users = data[:,0].astype(int)\n",
      "    items = data[:,1].astype(int)\n",
      "    ratings = data[:,2]\n",
      "    \n",
      "    nExamples = len(users)\n",
      "    \n",
      "    return (nExamples, users, items, ratings)\n",
      "\n",
      "TRAIN_FILE = '736/train'\n",
      "TEST_FILE  = '736/test'\n",
      "\n",
      "(nTrain, users, items, ratings) = loadData(TRAIN_FILE)\n",
      "(nTest,  testUsers, testItems, testRatings) = loadData(TEST_FILE)\n",
      "\n",
      "N = np.max(users)\n",
      "M = np.max(items)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Always look at our data!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure()\n",
      "plt.hist(ratings)\n",
      "plt.figure()\n",
      "plt.hist(testRatings)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "(array([ 47.,  18.,   9.,   2.,   0.,   3.,   1.,   0.,   0.,   1.]),\n",
        " array([ 0.75432862,  1.2940267 ,  1.83372478,  2.37342286,  2.91312094,\n",
        "        3.45281902,  3.9925171 ,  4.53221518,  5.07191326,  5.61161134,\n",
        "        6.15130942]),\n",
        " <a list of 10 Patch objects>)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEZFJREFUeJzt3G9MlfX/x/HXVXCnKUXf4qDnuJ3Gn/AoAsuwG7nhFFtr\n8aXZKCpjSXdqtVyusrYK20ps64b9cWtGjdWWeifwRjCqdcq5Jv3BtV/kpMIGRzhLkX1Ra6h8fjeq\ns0hBPVxw4Lyfj+1snOucc30+Hy98enFxjp5zzgkAYMYVqZ4AAGBmEX4AMIbwA4AxhB8AjCH8AGAM\n4QcAYyYNf19fn1atWqUlS5Zo6dKlev311yVJQ0NDqqysVGFhodauXavh4eHEa7Zu3aqCggIVFRWp\no6NjemcPALhs3mTv4x8cHNTg4KBKS0t18uRJ3XTTTWppadF7772n6667Tk8//bS2bdumEydOqLGx\nUd3d3brvvvv09ddfKxaLac2aNTp8+LCuuIIfLABgtpi0yLm5uSotLZUkzZs3T4sXL1YsFtPevXtV\nV1cnSaqrq1NLS4skqbW1VbW1tcrMzFQ4HFZ+fr46OzuneQkAgMtxyafiR44cUVdXl1asWKF4PK5A\nICBJCgQCisfjkqSjR48qFAolXhMKhRSLxXyeMgBgKi4p/CdPntS6deu0fft2zZ8/f9xjnufJ87wJ\nXzvZYwCAmZdxsSecOXNG69at0/r161VdXS3pz7P8wcFB5ebmamBgQDk5OZKkYDCovr6+xGv7+/sV\nDAbP22d+fr5+/vlnv9YAACbk5eXpp59+mvJ+Jj3jd86pvr5ekUhEGzduTGyvqqpSc3OzJKm5uTnx\nD0JVVZV27dql0dFR9fb2qqenR+Xl5eft9+eff5ZzLm1vL774YsrnwNpYH+tLv5tfJ8yTnvHv379f\nH3zwgZYtW6aysjJJf75dc/PmzaqpqVFTU5PC4bD27NkjSYpEIqqpqVEkElFGRoZ27NjBpR4AmGUm\nDf+tt96qsbGxCz726aefXnD7c889p+eee27qMwMATAveYD8NKioqUj2FaZPOa5NY31yX7uvzy6Qf\n4Jq2QT1PKRgWAOY0v9rJGT8AGEP4AcAYwg8AxhB+ADCG8AOAMYQfAIwh/ABgzEX/k7bp8s4776Rq\naN1www1avXp1ysYHgFRK2Qe4rrqqfqaHlSSNjQ3rP//5P/X3H0rJ+ACQLL8+wJWyM/7Tp1N1xn9I\nzlWnaGwASD2u8QOAMYQfAIwh/ABgDOEHAGMIPwAYQ/gBwBjCDwDGEH4AMIbwA4AxhB8AjCH8AGAM\n4QcAYwg/ABhD+AHAGMIPAMYQfgAwhvADgDGEHwCMIfwAYAzhBwBjCD8AGEP4AcAYwg8AxhB+ADCG\n8AOAMYQfAIwh/ABgDOEHAGMIPwAYQ/gBwBjCDwDGEH4AMIbwA4AxhB8AjCH8AGDMRcO/YcMGBQIB\nFRcXJ7Y1NDQoFAqprKxMZWVlamtrSzy2detWFRQUqKioSB0dHdMzawBA0i4a/oceekjt7e3jtnme\npyeffFJdXV3q6urS7bffLknq7u7W7t271d3drfb2dj366KMaGxubnpkDAJJy0fCvXLlS2dnZ5213\nzp23rbW1VbW1tcrMzFQ4HFZ+fr46Ozv9mSkAwBdJX+N/4403VFJSovr6eg0PD0uSjh49qlAolHhO\nKBRSLBab+iwBAL7JSOZFjzzyiF544QVJ0vPPP69Nmzapqanpgs/1PG+CvTT84+uKv24AgL9Fo1FF\no1Hf95tU+HNychJfP/zww7rzzjslScFgUH19fYnH+vv7FQwGJ9hLQzJDA4AZFRUVqqioSNzfsmWL\nL/tN6lLPwMBA4uuPPvoo8Y6fqqoq7dq1S6Ojo+rt7VVPT4/Ky8t9mSgAwB8XPeOvra3VF198oWPH\njmnRokXasmWLotGoDh48KM/zdMMNN+jtt9+WJEUiEdXU1CgSiSgjI0M7duyY5FIPACAVPHeht+dM\n96CeJ2nGh/3LIS1cWK1Y7FCKxgeA5Hied8F3VF4uPrkLAMYQfgAwhvADgDGEHwCMIfwAYAzhBwBj\nCD8AGEP4AcAYwg8AxhB+ADCG8AOAMYQfAIwh/ABgDOEHAGMIPwAYQ/gBwBjCDwDGEH4AMIbwA4Ax\nhB8AjCH8AGAM4QcAYwg/ABhD+AHAGMIPAMYQfgAwhvADgDGEHwCMIfwAYAzhBwBjCD8AGEP4AcAY\nwg8AxhB+ADCG8AOAMYQfAIwh/ABgDOEHAGMIPwAYQ/gBwBjCDwDGEH4AMIbwA4AxhB8AjCH8AGAM\n4QcAYwg/ABhz0fBv2LBBgUBAxcXFiW1DQ0OqrKxUYWGh1q5dq+Hh4cRjW7duVUFBgYqKitTR0TE9\nswYAJO2i4X/ooYfU3t4+bltjY6MqKyt1+PBhrV69Wo2NjZKk7u5u7d69W93d3Wpvb9ejjz6qsbGx\n6Zk5ACApFw3/ypUrlZ2dPW7b3r17VVdXJ0mqq6tTS0uLJKm1tVW1tbXKzMxUOBxWfn6+Ojs7p2Ha\nAIBkJXWNPx6PKxAISJICgYDi8bgk6ejRowqFQonnhUIhxWIxH6YJAPBLxlR34HmePM+b9PELa/jH\n1xV/3QAAf4tGo4pGo77vN6nwBwIBDQ4OKjc3VwMDA8rJyZEkBYNB9fX1JZ7X39+vYDA4wV4akhka\nAMyoqKhQRUVF4v6WLVt82W9Sl3qqqqrU3NwsSWpublZ1dXVi+65duzQ6Oqre3l719PSovLzcl4kC\nAPxx0TP+2tpaffHFFzp27JgWLVqkl156SZs3b1ZNTY2ampoUDoe1Z88eSVIkElFNTY0ikYgyMjK0\nY8eOSS8DAQBmnuecczM+qOdJmvFh/3JICxdWKxY7lKLxASA5nufJj2TzyV0AMIbwA4AxhB8AjCH8\nAGAM4QcAYwg/ABhD+AHAGMIPAMYQfgAwhvADgDGEHwCMIfwAYAzhBwBjCD8AGEP4AcAYwg8AxhB+\nADCG8AOAMYQfAIwh/ABgDOEHAGMIPwAYQ/gBwBjCDwDGEH4AMIbwA4AxhB8AjCH8AGAM4QcAYwg/\nABhD+AHAGMIPAMYQfgAwhvADgDGEHwCMIfwAYAzhBwBjCD8AGEP4AcAYwg8AxhB+ADCG8AOAMYQf\nAIwh/ABgDOEHAGMIPwAYQ/gBwJiMqbw4HA4rKytLV155pTIzM9XZ2amhoSHdc889+vXXXxUOh7Vn\nzx5dc801fs0XADBFUzrj9zxP0WhUXV1d6uzslCQ1NjaqsrJShw8f1urVq9XY2OjLRAEA/pjypR7n\n3Lj7e/fuVV1dnSSprq5OLS0tUx0CAOCjKZ/xr1mzRsuXL9fOnTslSfF4XIFAQJIUCAQUj8enPksA\ngG+mdI1///79WrBggX777TdVVlaqqKho3OOe58nzvAle3fCPryv+ugEA/haNRhWNRn3fr+f+fa0m\nSVu2bNG8efO0c+dORaNR5ebmamBgQKtWrdKhQ4fGD+p5knwZNgmHtHBhtWKxQxd/KgDMIp7nnXd5\nPRlJX+o5ffq0RkZGJEmnTp1SR0eHiouLVVVVpebmZklSc3OzqqurpzxJAIB/kr7UE4/Hddddd0mS\nzp49q/vvv19r167V8uXLVVNTo6ampsTbOQEAs4dvl3oua1Au9QDAZUv5pR4AwNxE+AHAGMIPAMYQ\nfgAwhvADgDGEHwCMIfwAYAzhBwBjCD8AGEP4AcAYwg8AxhB+ADCG8AOAMYQfAIwh/ABgDOEHAGMI\nPwAYQ/gBwBjCDwDGEH4AMIbwA4AxhB8AjCH8AGAM4QcAYwg/ABhD+AHAGMIPAMYQfgAwhvADgDGE\nHwCMIfwAYAzhBwBjCD8AGEP4AcAYwg8AxhB+ADCG8AOAMYQfAIwh/ABgDOEHAGMIPwAYQ/gBwBjC\nDwDGEH4AMIbwA4AxhB8AjJmW8Le3t6uoqEgFBQXatm3bdAwBAEiS7+E/d+6cHnvsMbW3t6u7u1sf\nfvihfvzxR7+HmdWi0WiqpzBt0nltEuub69J9fX7xPfydnZ3Kz89XOBxWZmam7r33XrW2tvo9zKyW\nzt986bw2ifXNdem+Pr/4Hv5YLKZFixYl7odCIcViMb+HQRrIyrpWnuel7JaVdW2q/wiAlMjwe4ee\n513S87Ky7vR76EsyNjaijAx+pz0bjIyckORSOP6lfa8ivWVlXfvX92JqzJ+frf/9b2hmB3U+++qr\nr9xtt92WuP/KK6+4xsbGcc/Jy8tz+vNvPDdu3Lhxu8RbXl6eL532nHNOPjp79qxuvPFGffbZZ1q4\ncKHKy8v14YcfavHixX4OAwBIku+XejIyMvTmm2/qtttu07lz51RfX0/0AWAW8f2MHwAwu83obznT\n5YNd4XBYy5YtU1lZmcrLyyVJQ0NDqqysVGFhodauXavh4eHE87du3aqCggIVFRWpo6MjVdOe0IYN\nGxQIBFRcXJzYlsx6vv32WxUXF6ugoEBPPPHEjK5hIhdaW0NDg0KhkMrKylRWVqa2trbEY3NpbZLU\n19enVatWacmSJVq6dKlef/11Selz/CZaX7ocwz/++EMrVqxQaWmpIpGInn32WUkzcPx8+U3BJTh7\n9qzLy8tzvb29bnR01JWUlLju7u6ZGt5X4XDYHT9+fNy2p556ym3bts0551xjY6N75plnnHPO/fDD\nD66kpMSNjo663t5el5eX586dOzfjc57Ml19+6b777ju3dOnSxLbLWc/Y2Jhzzrmbb77ZHThwwDnn\n3O233+7a2tpmeCXnu9DaGhoa3GuvvXbec+fa2pxzbmBgwHV1dTnnnBsZGXGFhYWuu7s7bY7fROtL\np2N46tQp55xzZ86ccStWrHD79u2b9uM3Y2f86fbBLvevK2R79+5VXV2dJKmurk4tLS2SpNbWVtXW\n1iozM1PhcFj5+fnq7Oyc8flOZuXKlcrOzh637XLWc+DAAQ0MDGhkZCTxE9CDDz6YeE0qXWht0vnH\nT5p7a5Ok3NxclZaWSpLmzZunxYsXKxaLpc3xm2h9Uvocw6uuukqSNDo6qnPnzik7O3vaj9+MhT+d\nPtjleZ7WrFmj5cuXa+fOnZKkeDyuQCAgSQoEAorH45Kko0ePKhQKJV47V9Z9uev59/ZgMDir1/nG\nG2+opKRE9fX1iR+j5/rajhw5oq6uLq1YsSItj9/f67vlllskpc8xHBsbU2lpqQKBQOKy1nQfvxkL\n/6V+sGsu2L9/v7q6utTW1qa33npL+/btG/f4358Mnchc+7O42HrmmkceeUS9vb06ePCgFixYoE2b\nNqV6SlN28uRJrVu3Ttu3b9f8+fPHPZYOx+/kyZO6++67tX37ds2bNy+tjuEVV1yhgwcPqr+/X19+\n+aU+//zzcY9Px/GbsfAHg0H19fUl7vf19Y37F2ouWbBggSTp+uuv11133aXOzk4FAgENDg5KkgYG\nBpSTkyPp/HX39/crGAzO/KQv0+WsJxQKKRgMqr+/f9z22brOnJycxF+mhx9+OHHpba6u7cyZM1q3\nbp3Wr1+v6upqSel1/P5e3wMPPJBYX7odQ0m6+uqrdccdd+jbb7+d9uM3Y+Ffvny5enp6dOTIEY2O\njmr37t2qqqqaqeF9c/r0aY2MjEiSTp06pY6ODhUXF6uqqkrNzc2SpObm5sQ3aFVVlXbt2qXR0VH1\n9vaqp6cncR1uNrvc9eTm5iorK0sHDhyQc07vv/9+4jWzzcDAQOLrjz76KPGOn7m4Nuec6uvrFYlE\ntHHjxsT2dDl+E60vXY7hsWPHEpepfv/9d33yyScqKyub/uPn/++oJ/bxxx+7wsJCl5eX51555ZWZ\nHNo3v/zyiyspKXElJSVuyZIliXUcP37crV692hUUFLjKykp34sSJxGtefvlll5eX52688UbX3t6e\nqqlP6N5773ULFixwmZmZLhQKuXfffTep9XzzzTdu6dKlLi8vzz3++OOpWMp5/r22pqYmt379eldc\nXOyWLVvm/vvf/7rBwcHE8+fS2pxzbt++fc7zPFdSUuJKS0tdaWmpa2trS5vjd6H1ffzxx2lzDL//\n/ntXVlbmSkpKXHFxsXv11Vedc8n15HLWxwe4AMAY/ptKADCG8AOAMYQfAIwh/ABgDOEHAGMIPwAY\nQ/gBwBjCDwDG/D8q4jCfn0CZXQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x111c21e50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAAEACAYAAAB8nvebAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEXRJREFUeJzt3V1sVOXah/H/wtbNVj6khK5BIamptlKodhRtYjQO1ika\nAykbBBFxAuqJ8UCTd2s9EjzQMcYY/DgyxFR5A5ujWg1WqTCCEFO1JSZqqBrQgu1E7R6oBVLarveA\n12otdD6Y1cVdrl8ySZlZM8/dBi5W1swDjud5ngAApkwKegAAQPaINwAYRLwBwCDiDQAGEW8AMIh4\nA4BBBZkcVFJSomnTpumSSy5RYWGhWltb1dPTo1WrVunHH39USUmJtm/friuuuMLveQEAyvDM23Ec\nJRIJtbe3q7W1VZIUj8cVjUbV0dGhmpoaxeNxXwcFAPwp48smf9/L09TUpFgsJkmKxWJqbGzM72QA\ngHPK+Mz7rrvu0sKFC/Xmm29KkpLJpFzXlSS5rqtkMunflACAETK65r1v3z7Nnj1bv/zyi6LRqK67\n7roRjzuOI8dxfBkQADBaRvGePXu2JGnWrFlatmyZWltb5bquuru7FQqF1NXVpeLi4lHPu+aaa/TD\nDz/kd2IAmOBKS0v1/fffj3lM2ssmJ06cUG9vrySpr69PH330kSorK7V06VI1NDRIkhoaGlRXVzfq\nuT/88IM8zzN7e/bZZwOf4WKcnfmDvzF/sLdMTnrTnnknk0ktW7ZMkjQwMKA1a9aotrZWCxcu1MqV\nK7V58+bhjwoCAMZH2nhfffXVOnDgwKj7i4qK1NLS4stQAICxscNyDJFIJOgRcmZ5don5g8b8Fz7H\n8zzf/jMGx3Hk48sDwISUSTs58wYAg4g3ABhEvAHAIOINAAYRbwAwiHgDgEEZ/dsmE0FLS4s++WSP\n7+tceuml+ve//0eTJ0/2fS0AF6+L5nPe99yzSs3NQ5IqfV3n0ks36auv9qu8vNzXdQBMXJm086I5\n8z5jhaRVvq7wj3/8r6+vDwAS17wBwCTiDQAGEW8AMIh4A4BBxBsADCLeAGAQ8QYAg4g3ABhEvAHA\nIOINAAYRbwAwiHgDgEHEGwAMIt4AYBDxBgCDiDcAGES8AcAg4g0ABhFvADCIeAOAQcQbAAwi3gBg\nEPEGAIOINwAYRLwBwKCM4j04OKhwOKwlS5ZIknp6ehSNRlVWVqba2lqlUilfhwQAjJRRvDdt2qSK\nigo5jiNJisfjikaj6ujoUE1NjeLxuK9DAgBGShvvI0eOaMeOHXrkkUfkeZ4kqampSbFYTJIUi8XU\n2Njo75QAgBHSxvvJJ5/USy+9pEmT/jw0mUzKdV1Jkuu6SiaT/k0IABilYKwH33//fRUXFyscDiuR\nSJz1GMdxhi+nnM2GDRuGv45EIopEIrnMCQATViKROGdjz2XMeO/fv19NTU3asWOHTp06pePHj2vt\n2rVyXVfd3d0KhULq6upScXHxOV/jr/EGAIz29xPbjRs3pn3OmJdNnn/+eXV2durQoUPatm2b7rzz\nTr3zzjtaunSpGhoaJEkNDQ2qq6s7v8kBAFnJ6nPef1weqa+v186dO1VWVqZdu3apvr7el+EAAGc3\n5mWTv7rjjjt0xx13SJKKiorU0tLi21AAgLGxwxIADCLeAGAQ8QYAg4g3ABhEvAHAIOINAAYRbwAw\niHgDgEHEGwAMIt4AYBDxBgCDiDcAGES8AcAg4g0ABhFvADCIeAOAQcQbAAwi3gBgEPEGAIOINwAY\nRLwBwCDiDQAGEW8AMIh4A4BBxBsADCLeAGAQ8QYAg4g3ABhEvAHAIOINAAYRbwAwiHgDgEHEGwAM\nIt4AYBDxBgCDxoz3qVOnVF1draqqKlVUVOiZZ56RJPX09CgajaqsrEy1tbVKpVLjMiwA4Iwx4z15\n8mTt3r1bBw4c0FdffaXdu3fr008/VTweVzQaVUdHh2pqahSPx8drXgCAMrhsctlll0mS+vv7NTg4\nqBkzZqipqUmxWEySFIvF1NjY6O+UAIAR0sZ7aGhIVVVVcl1XixYt0vz585VMJuW6riTJdV0lk0nf\nBwUA/Kkg3QGTJk3SgQMHdOzYMS1evFi7d+8e8bjjOHIcx7cBAQCjpY33H6ZPn657771XX375pVzX\nVXd3t0KhkLq6ulRcXHzO523YsGH460gkokgkcj7zAsCEk0gklEgksnqO43med64Hf/31VxUUFOiK\nK67QyZMntXjxYj377LP68MMPNXPmTD399NOKx+NKpVJnfdPScRyN8fLj6p57Vqm5+V+SVvm6ztSp\n5fr88yaVl5f7ug6AiSuTdo555t3V1aVYLKahoSENDQ1p7dq1qqmpUTgc1sqVK7V582aVlJRo+/bt\neR0cADC2MeNdWVmptra2UfcXFRWppaXFt6EAAGNjhyUAGES8AcAg4g0ABhFvADCIeAOAQcQbAAwi\n3gBgEPEGAIOINwAYRLwBwCDiDQAGEW8AMIh4A4BBxBsADCLeAGAQ8QYAg4g3ABhEvAHAIOINAAYR\nbwAwiHgDgEHEGwAMIt4AYBDxBgCDiDcAGES8AcAg4g0ABhFvADCIeAOAQcQbAAwi3gBgEPEGAIOI\nNwAYRLwBwCDiDQAGEW8AMChtvDs7O7Vo0SLNnz9fCxYs0KuvvipJ6unpUTQaVVlZmWpra5VKpXwf\nFgBwRtp4FxYW6pVXXtHXX3+tzz77TG+88Ya+/fZbxeNxRaNRdXR0qKamRvF4fDzmBQAog3iHQiFV\nVVVJkqZMmaJ58+bp6NGjampqUiwWkyTFYjE1Njb6OykAYFhW17wPHz6s9vZ2VVdXK5lMynVdSZLr\nukomk74MCAAYLeN4//7771q+fLk2bdqkqVOnjnjMcRw5jpP34QAAZ1eQyUGnT5/W8uXLtXbtWtXV\n1Uk6c7bd3d2tUCikrq4uFRcXn/W5GzZsGP46EokoEomc99AAMJEkEgklEomsnuN4nueNdYDneYrF\nYpo5c6ZeeeWV4fufeuopzZw5U08//bTi8bhSqdSoNy0dx1Galx8399yzSs3N/5K0ytd1pk4t1+ef\nN6m8vNzXdQBMXJm0M+2Z9759+7RlyxZdf/31CofDkqQXXnhB9fX1WrlypTZv3qySkhJt3749P1MD\nANJKG+/bbrtNQ0NDZ32spaUl7wMBANJjhyUAGES8AcAg4g0ABhFvADCIePvgppuqhzcu+XWbNq0o\n6G8TQIAy2qSD7PT1HZPk7+fbe3vZ0QpczDjzBgCDiDcAGES8AcAg4g0ABhFvADCIeAOAQcQbAAwi\n3gBgEPEGAIOINwAYRLwBwCDiDQAGEW8AMIh4A4BBxBsADCLeAGAQ8QYAg4g3ABhEvAHAIOINAAYR\nbwAwiHgDgEHEGwAMIt4AYBDxBgCDiDcAGES8AcAg4g0ABhFvADCIeAOAQWnjvX79ermuq8rKyuH7\nenp6FI1GVVZWptraWqVSKV+HBACMlDbe69atU3Nz84j74vG4otGoOjo6VFNTo3g87tuAAIDR0sb7\n9ttv14wZM0bc19TUpFgsJkmKxWJqbGz0ZzoAwFnldM07mUzKdV1Jkuu6SiaTeR0KADC2gvN9Acdx\n5DjOOR/fsGHD8NeRSESRSOR8lwSACSWRSCiRSGT1nJzi7bquuru7FQqF1NXVpeLi4nMe+9d4AwBG\n+/uJ7caNG9M+J6fLJkuXLlVDQ4MkqaGhQXV1dbm8DAAgR2njvXr1at166606ePCg5s6dq7feekv1\n9fXauXOnysrKtGvXLtXX14/HrACA/5f2ssnWrVvPen9LS0vehwEAZOa837BEUArGfKM4X6ZOnaHj\nx3t8XwdAdoi3WQOSPN9X6e31/y8IANnj3zYBAIOINwAYRLwBwCDiDQAGEW8AMIh4A4BBxBsADCLe\nAGAQ8QYAg4g3ABhEvAHAIOINAAYRbwAwiHgDgEHEGwAMIt4AYBDxBgCDiDcAGES8AcAg4g0ABhFv\nADCIeAOAQcQbAAwi3gBgEPHGRWXatCI5juP7bdq0oqC/VUxwBUEPAIyn3t7/SvLGYR3H9zVwcePM\nGwAMIt4AYBDxBgCDuOaNNArkOP5fv506dYaOH+/xfR1goiDeSGNAvMEHXHi4bAIABhFvADDovOLd\n3Nys6667Ttdee61efPHFfM0ETAAFbAbKEhuospNzvAcHB/X444+rublZ33zzjbZu3apvv/02n7Nd\nABJBD3AeEkEPcF4SiUTQI5ynP94r8Pd2ZtNR/gXx8/9zA1U+brvP+ZhfP7PxlnO8W1tbdc0116ik\npESFhYW6//779e677+ZztgtAIugBzkMi6AHOi/1422b/558IegDf5Rzvo0ePau7cucO/njNnjo4e\nPZqXoQAAY8v5o4Lj8dnffCosnKR//vNFFRZuyfg5p04d1OTJX2a1zsmTR7IdDQCylnO8r7rqKnV2\ndg7/urOzU3PmzBlxTGlp6QUX+ZMn27M6vr//uxxXGo/vO90aG8dpnTyt8rffKxs35mv+USv59LrB\nrOPXnzH/fv5jyef3cu75L7Qu/V1paWnaYxzP83LagTEwMKDy8nJ9/PHHuvLKK3XLLbdo69atmjdv\nXi4vBwDIQs5n3gUFBXr99de1ePFiDQ4O6uGHHybcADBOcj7zBgAEx5cdlpY376xfv16u66qysjLo\nUXLS2dmpRYsWaf78+VqwYIFeffXVoEfKyqlTp1RdXa2qqipVVFTomWeeCXqkrA0ODiocDmvJkiVB\nj5KTkpISXX/99QqHw7rllluCHicrqVRKK1as0Lx581RRUaHPPvss6JEydvDgQYXD4eHb9OnTx/7z\n6+XZwMCAV1pa6h06dMjr7+/3brjhBu+bb77J9zK+2bNnj9fW1uYtWLAg6FFy0tXV5bW3t3ue53m9\nvb1eWVmZqZ+/53leX1+f53med/r0aa+6utrbu3dvwBNl5+WXX/YeeOABb8mSJUGPkpOSkhLvt99+\nC3qMnDz00EPe5s2bPc878/snlUoFPFFuBgcHvVAo5P3000/nPCbvZ97WN+/cfvvtmjFjRtBj5CwU\nCqmqqkqSNGXKFM2bN08///xzwFNl57LLLpMk9ff3a3BwUEVFdrYzHzlyRDt27NAjjzwiz/AVSYuz\nHzt2THv37tX69eslnXlfbvr06QFPlZuWlhaVlpaO2Evzd3mPN5t3LhyHDx9We3u7qqurgx4lK0ND\nQ6qqqpLrulq0aJEqKiqCHiljTz75pF566SVNmmT333xzHEd33XWXFi5cqDfffDPocTJ26NAhzZo1\nS+vWrdONN96oRx99VCdOnAh6rJxs27ZNDzzwwJjH5P132IX++cmLxe+//64VK1Zo06ZNmjJlStDj\nZGXSpEk6cOCAjhw5oj179pjZqv3++++ruLhY4XDY5JnrH/bt26f29nZ98MEHeuONN7R3796gR8rI\nwMCA2tra9Nhjj6mtrU2XX3654vF40GNlrb+/X++9957uu+++MY/Le7wz2bwDf50+fVrLly/Xgw8+\nqLq6uqDHydn06dN177336osvvgh6lIzs379fTU1Nuvrqq7V69Wrt2rVLDz30UNBjZW327NmSpFmz\nZmnZsmVqbW0NeKLMzJkzR3PmzNHNN98sSVqxYoXa2toCnip7H3zwgW666SbNmjVrzOPyHu+FCxfq\nu+++0+HDh9Xf36///Oc/Wrp0ab6XwTl4nqeHH35YFRUVeuKJJ4IeJ2u//vqrUqmUJOnkyZPauXOn\nwuFwwFNl5vnnn1dnZ6cOHTqkbdu26c4779Tbb78d9FhZOXHihHp7eyVJfX19+uijj8x88ioUCmnu\n3Lnq6OiQdOa68fz58wOeKntbt27V6tWr0x6X9/8GzfrmndWrV+uTTz7Rb7/9prlz5+q5557TunXr\ngh4rY/v27dOWLVuGP+olSS+88ILuvvvugCfLTFdXl2KxmIaGhjQ0NKS1a9eqpqYm6LFyYvESYjKZ\n1LJlyySduQyxZs0a1dbWBjxV5l577TWtWbNG/f39Ki0t1VtvvRX0SFnp6+tTS0tLRu81sEkHAAyy\n+5Y4AFzEiDcAGES8AcAg4g0ABhFvADCIeAOAQcQbAAwi3gBg0P8BYKDK10vtenMAAAAASUVORK5C\nYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x111c21890>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From the test set, it appears like we have 1-6 stars, but something weird is happening in the training set! The histogram says there are outliers, so let's take a look:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "outlierIdxs = ratings > 6\n",
      "print np.flatnonzero(outlierIdxs)\n",
      "print ratings[outlierIdxs]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 83 163 165 167 168]\n",
        "[    9.54361436  1379.98131862     8.76865101  2907.19215935  1380.47446215]\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I delete the outliers from my data file; much better. **Always** make a new copy; do not overwrite your original data!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "CLEAN_TRAIN_FILE = '736/trainClean'\n",
      "(nTrain, users, items, ratings) = loadData(CLEAN_TRAIN_FILE)\n",
      "\n",
      "plt.hist(ratings)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "(array([ 95.,  43.,  17.,   5.,   5.,   3.,  13.,   1.,   1.,   2.]),\n",
        " array([ 0.75722523,  1.23237711,  1.70752898,  2.18268086,  2.65783273,\n",
        "        3.13298461,  3.60813648,  4.08328836,  4.55844023,  5.03359211,\n",
        "        5.50874398]),\n",
        " <a list of 10 Patch objects>)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEhRJREFUeJzt3W9sFPW+x/HPcNrkRP5D6K5QTE1pbRcQSmoxRs1i3R7/\nQUBrI6BsKOoDYqLEXEVjYvWBXaNGwT8xUWJqTEAe1WqAaDULyF+9xWgEbylpY8GyN9qzUqlY2s59\nwLVHpC3tdnaXfnm/kialndnfd4p5d5zdYR3XdV0BAEa9MekeAADgDYIOAEYQdAAwgqADgBEEHQCM\nIOgAYMSgQa+srJTP59PcuXP7vtbe3q5QKKT8/HyVlZUpHo/3fa+6ulp5eXkqKCjQp59+mrypAQAX\nGDToq1ev1o4dO877WiQSUSgUUmNjo0pLSxWJRCRJhw8f1ocffqjDhw9rx44dWrt2rXp7e5M3OQDg\nPIMG/aabbtLkyZPP+1pdXZ3C4bAkKRwOq7a2VpL00Ucfafny5crMzFROTo5mzZqlgwcPJmlsAMDf\nDfsaeiwWk8/nkyT5fD7FYjFJ0k8//aTs7Oy+7bKzs3XixAmPxgQAXMyInhR1HEeO4wz6fQBAamQM\ndwefz6eTJ0/K7/erra1NWVlZkqQZM2aotbW1b7vjx49rxowZF+w/a9YsHTt2bAQjA8DlJzc3V01N\nTYNuM+wz9CVLlqimpkaSVFNTo6VLl/Z9fcuWLerq6lJzc7OOHj2qkpKSC/Y/duyYXNc1+/Hss8+m\nfQaOj+O7HI/P8rG5rjukE+FBz9CXL1+unTt36ueff9bMmTP1/PPPa/369aqoqNCmTZuUk5OjrVu3\nSpICgYAqKioUCASUkZGht956i0suAJBCgwZ98+bN/X69vr6+368//fTTevrpp0c+FQBg2LhT1GPB\nYDDdIyQVxze6WT4+y8c2VI7ruil9gwvHcZTiJQFg1BtKOzlDBwAjCDoAGEHQAcAIgg4ARgz7TtHL\nxeOPP6NodF9K1lq27DY988x/pWQtAHbxKpcB5OUVq6mpUtI1SV7poK6/fpf27due5HUAjGZDaSdn\n6IMqkVSc5DW6JO1K8hoALgdcQwcAIwg6ABhB0AHACIIOAEYQdAAwgqADgBEEHQCMIOgAYARBBwAj\nCDoAGEHQAcAIgg4ARhB0ADCCoAOAEQQdAIwg6ABgBEEHACMIOgAYQdABwAiCDgBGEHQAMIKgA4AR\nBB0AjCDoAGAEQQcAIwg6ABhB0AHAiISDXl1drdmzZ2vu3LlasWKF/vjjD7W3tysUCik/P19lZWWK\nx+NezgoAGERCQW9padE777yjhoYGfffdd+rp6dGWLVsUiUQUCoXU2Nio0tJSRSIRr+cFAAwgoaBP\nmDBBmZmZ6uzsVHd3tzo7OzV9+nTV1dUpHA5LksLhsGpraz0dFgAwsISCPmXKFD3++OO66qqrNH36\ndE2aNEmhUEixWEw+n0+S5PP5FIvFPB0WADCwjER2OnbsmF577TW1tLRo4sSJuvfee/XBBx+ct43j\nOHIcp9/9q6qq+j4PBoMKBoOJjAEAZkWjUUWj0WHtk1DQv/76a91www2aOnWqJOnuu+/Wvn375Pf7\ndfLkSfn9frW1tSkrK6vf/f8adADAhf5+svvcc89ddJ+ELrkUFBRo//79+v333+W6rurr6xUIBLR4\n8WLV1NRIkmpqarR06dJEHh4AkICEztDnzZunVatWqbi4WGPGjNGCBQv08MMPq6OjQxUVFdq0aZNy\ncnK0detWr+cFAAzAcV3XTemCjqMUL5mQvLxiNTW9Lak4yStt1/XXb9S+fduTvA6A0Wwo7eROUQAw\ngqADgBEEHQCMIOgAYARBBwAjCDoAGEHQAcAIgg4ARhB0ADCCoAOAEQQdAIwg6ABgBEEHACMIOgAY\nQdABwAiCDgBGEHQAMIKgA4ARBB0AjCDoAGAEQQcAIwg6ABhB0AHACIIOAEYQdAAwgqADgBEEHQCM\nIOgAYARBBwAjCDoAGEHQAcAIgg4ARhB0ADCCoAOAEQQdAIxIOOjxeFzl5eUqLCxUIBDQgQMH1N7e\nrlAopPz8fJWVlSkej3s5KwBgEAkH/dFHH9Udd9yhI0eO6Ntvv1VBQYEikYhCoZAaGxtVWlqqSCTi\n5awAgEEkFPRff/1Vu3fvVmVlpSQpIyNDEydOVF1dncLhsCQpHA6rtrbWu0kBAINKKOjNzc2aNm2a\nVq9erQULFuihhx7S6dOnFYvF5PP5JEk+n0+xWMzTYQEAA0so6N3d3WpoaNDatWvV0NCgsWPHXnB5\nxXEcOY7jyZAAgIvLSGSn7OxsZWdn67rrrpMklZeXq7q6Wn6/XydPnpTf71dbW5uysrL63b+qqqrv\n82AwqGAwmMgYAGBWNBpVNBod1j6O67puIovdfPPNevfdd5Wfn6+qqip1dnZKkqZOnaonn3xSkUhE\n8Xi83zP3BJdMqby8YjU1vS2pOMkrbdf112/Uvn3bk7wOgNFsKO1M6Axdkl5//XWtXLlSXV1dys3N\n1Xvvvaeenh5VVFRo06ZNysnJ0datWxN9eADAMCUc9Hnz5umrr7664Ov19fUjGggAkBjuFAUAIwg6\nABhB0AHACIIOAEYQdAAwgqADgBEEHQCMIOgAYARBBwAjCDoAGEHQAcAIgg4ARhB0ADCCoAOAEQQd\nAIwg6ABgBEEHACMIOgAYQdABwAiCDgBGEHQAMIKgA4ARBB0AjCDoAGAEQQcAIwg6ABhB0AHACIIO\nAEYQdAAwgqADgBEEHQCMIOgAYARBBwAjCDoAGJGR7gEgHTwYleM4KVlr/PjJOnWqPSVrAUgtgn4J\n6O09I8lNyVodHan5xQEg9UZ0yaWnp0dFRUVavHixJKm9vV2hUEj5+fkqKytTPB73ZEgAwMWNKOgb\nNmxQIBDou1wQiUQUCoXU2Nio0tJSRSIRT4YEAFxcwkE/fvy4tm3bpgcffFCue+5yQV1dncLhsCQp\nHA6rtrbWmykBABeVcNDXrVunl156SWPG/OchYrGYfD6fJMnn8ykWi418QgDAkCQU9E8++URZWVkq\nKirqOzv/O8dxUvbKDQBAgq9y2bt3r+rq6rRt2zadOXNGp06d0gMPPCCfz6eTJ0/K7/erra1NWVlZ\n/e5fVVXV93kwGFQwGExkDAAwKxqNKhqNDmsfxx3oFHuIdu7cqZdfflkff/yxnnjiCU2dOlVPPvmk\nIpGI4vH4BU+MOo4z4Fn9pSQvr1hNTW9LKk7yStsl3aFUvWxRGh0/fwDnG0o7PblT9M9LK+vXr9dn\nn32m/Px8ffHFF1q/fr0XDw8AGIIRn6EPe0HO0P+GM3QAF5eyM3QAQPoRdAAwgqADgBEEHQCMIOgA\nYARBBwAjCDoAGEHQAcAIgg4ARhB0ADCCoAOAEQQdAIwg6ABgBEEHACMIOgAYQdABwAiCDgBGEHQA\nMIKgA4ARBB0AjCDoAGAEQQcAIwg6ABhB0AHACIIOAEYQdAAwgqADgBEEHQCMIOgAYARBBwAjCDoA\nGEHQAcAIgg4ARhB0ADCCoAOAEQQdAIxIKOitra1atGiRZs+erTlz5mjjxo2SpPb2doVCIeXn56us\nrEzxeNzTYQEAA0so6JmZmXr11Vf1/fffa//+/XrzzTd15MgRRSIRhUIhNTY2qrS0VJFIxOt5AQAD\nSCjofr9f8+fPlySNGzdOhYWFOnHihOrq6hQOhyVJ4XBYtbW13k0KABjUiK+ht7S06NChQ1q4cKFi\nsZh8Pp8kyefzKRaLjXhAAMDQZIxk599++0333HOPNmzYoPHjx5/3Pcdx5DhOv/tVVVX1fR4MBhUM\nBkcyBgCYE41GFY1Gh7WP47qum8hiZ8+e1V133aXbb79djz32mCSpoKBA0WhUfr9fbW1tWrRokX74\n4YfzF3QcJbhkSuXlFaup6W1JxUleabukOySl6mcyOn7+AM43lHYmdMnFdV2tWbNGgUCgL+aStGTJ\nEtXU1EiSampqtHTp0kQeHgCQgITO0L/88kvdfPPNuvbaa/suq1RXV6ukpEQVFRX68ccflZOTo61b\nt2rSpEnnL8gZ+t9whg7g4obSzoSuod94443q7e3t93v19fWJPCQAYIS4UxQAjCDoAGDEiF62iNEo\nY8CXk3pp/PjJOnWqPenrAPgPgn7Z6VYqnoDt6Ej+Lw0A5+OSCwAYQdABwAiCDgBGEHQAMIKgA4AR\nBB0AjCDoAGAEQQcuQRMmTOl7T4FkfkyYMCXdhwoPcWMRcAnq6Pi3uAEMw8UZOgAYQdABwAiCDgBG\nEHQAMIKgA4ARBB0AjCDoAGAEQQcAIwg6ABhB0AHACIIOAEYQdAAwgqADgBEEHQCMIOgAYARBBwAj\neIMLJEmGHCdVb56QKels0lcZP36yTp1qT/o6QKIIOpKkW6l4x51znJSsxbv74FLHJRcAMIKgA4AR\nXHIBhiyVzwsAw0fQgSFL9fMCwPB4fsllx44dKigoUF5enl588UWvHx4ABjVhwhQ5jpOSjwkTpqT7\ncM/jadB7enr0yCOPaMeOHTp8+LA2b96sI0eOeLnEKBBN9wBJFk33AEkWTfcASRZN9wBJE41GJUkd\nHf/Wuf+TSv7HubUuHZ4G/eDBg5o1a5ZycnKUmZmp++67Tx999JGXS4wC0XQPkGTRdA+QZNF0D5Bk\n0XQPkDR/Bv1y5uk19BMnTmjmzJl9f87OztaBAwe8XAKAp2zdAPbcc88l9fEvdZ4G3dIrADIzx2js\n2HX6xz8mDWu/M2f+R//8538Pefvu7v9VZ+dwpwO8YukGsKr//7DToeHyNOgzZsxQa2tr359bW1uV\nnZ193ja5ubmmwt+frq6jCeyVyp/JSNca6lnQaDqmvxrs+EbrMf3V34/PwjH96c9jS90xpapnubm5\nF93GcV3Xs1+Z3d3duuaaa/T5559r+vTpKikp0ebNm1VYWOjVEgCAAXh6hp6RkaE33nhD//rXv9TT\n06M1a9YQcwBIEU/P0AEA6ZPSf8vF8k1HlZWV8vl8mjt3brpHSYrW1lYtWrRIs2fP1pw5c7Rx48Z0\nj+SZM2fOaOHChZo/f74CgYCeeuqpdI+UFD09PSoqKtLixYvTPYrncnJydO2116qoqEglJSXpHsdz\n8Xhc5eXlKiwsVCAQ0P79+/vf0E2R7u5uNzc3121ubna7urrcefPmuYcPH07V8km3a9cut6GhwZ0z\nZ066R0mKtrY299ChQ67rum5HR4ebn59v6u/v9OnTruu67tmzZ92FCxe6u3fvTvNE3nvllVfcFStW\nuIsXL073KJ7Lyclxf/nll3SPkTSrVq1yN23a5Lruuf9G4/F4v9ul7Azd+k1HN910kyZPnpzuMZLG\n7/dr/vz5kqRx48apsLBQP/30U5qn8s4VV1whSerq6lJPT4+mTLm0bukeqePHj2vbtm168MEH5Rq9\nymr1uH799Vft3r1blZWVks49Vzlx4sR+t01Z0Pu76ejEiROpWh4eamlp0aFDh7Rw4cJ0j+KZ3t5e\nzZ8/Xz6fT4sWLVIgEEj3SJ5at26dXnrpJY0ZY/NfzHYcR7feequKi4v1zjvvpHscTzU3N2vatGla\nvXq1FixYoIceekidA9y8krK/XeuvPb9c/PbbbyovL9eGDRs0bty4dI/jmTFjxuibb77R8ePHtWvX\nLlO3kX/yySfKyspSUVGR2bPYPXv26NChQ9q+fbvefPNN7d69O90jeaa7u1sNDQ1au3atGhoaNHbs\nWEUikX63TVnQh3LTES5tZ8+e1T333KP7779fS5cuTfc4STFx4kTdeeed+vrrr9M9imf27t2ruro6\nXX311Vq+fLm++OILrVq1Kt1jeerKK6+UJE2bNk3Lli3TwYMH0zyRd7Kzs5Wdna3rrrtOklReXq6G\nhoZ+t01Z0IuLi3X06FG1tLSoq6tLH374oZYsWZKq5TFCrutqzZo1CgQCeuyxx9I9jqd+/vlnxeNx\nSdLvv/+uzz77TEVFRWmeyjsvvPCCWltb1dzcrC1btuiWW27R+++/n+6xPNPZ2amOjg5J0unTp/Xp\np5+aerWZ3+/XzJkz1djYKEmqr6/X7Nmz+902ZW9wYf2mo+XLl2vnzp365ZdfNHPmTD3//PNavXp1\nusfyzJ49e/TBBx/0vTRMkqqrq3XbbbelebKRa2trUzgcVm9vr3p7e/XAAw+otLQ03WMljbXLn7FY\nTMuWLZN07vLEypUrVVZWluapvPX6669r5cqV6urqUm5urt57771+t+PGIgAwwuZT3gBwGSLoAGAE\nQQcAIwg6ABhB0AHACIIOAEYQdAAwgqADgBH/B8SiM+BzQwRDAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1128ed090>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that we've taken a look at our data, let's tell a story.\n",
      "\n",
      "The classic _matrix factorization_ method for recommendations represents users and items as $d$-dimensional vectors.\n",
      "\n",
      "- Each dimension represents some abstract \"taste\": horror? chick flick? comedy?\n",
      "- A rating $r_{ij}$ is just the dot product of the respective user and item vectors: $\\mathbf{u}_i ^\\top \\mathbf{v}_j$.\n",
      "\n",
      "In Bayesian land,\n",
      "\n",
      "- For user $i = 1, \\ldots, N$, draw $\\mathbf{u}_i \\sim \\mathcal{N}(0, \\tau_0 I_d)$\n",
      "- For item $j = 1, \\ldots, M$ draw $\\mathbf{v}_j \\sim \\mathcal{N}(0, \\tau_0 I_d)$\n",
      "- For $i = 1, \\ldots, N$ and $j = 1, \\ldots, M$,\n",
      "\n",
      "    - Compute $z_{ij} = \\mathbf{u}_i ^\\top \\mathbf{v}_j$\n",
      "    - Draw $r_{ij} \\sim \\mathcal{LN}(z_{ij}, \\tau_1)$.\n",
      "\n",
      "The Bayesian network looks like\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import rc\n",
      "rc(\"font\", family=\"serif\", size=36)\n",
      "rc(\"text\", usetex=True)\n",
      "\n",
      "pgm = daft.PGM([5, 4], grid_unit=5, node_unit=3)\n",
      "pgm.add_node(daft.Node(\"tau0\", r\"$\\tau_0$\",    2, 3.5, fixed=True))\n",
      "pgm.add_node(daft.Node(\"u\", r\"$\\mathbf{u}_i$\", 1, 2.5))\n",
      "pgm.add_node(daft.Node(\"v\", r\"$\\mathbf{v}_j$\", 3, 2.5))\n",
      "\n",
      "pgm.add_node(daft.Node(\"z\", r\"$z_{ij}$\", 2, 2.5))\n",
      "pgm.add_node(daft.Node(\"r\", r\"$r_{ij}$\", 2, 1.5, observed=True))\n",
      "\n",
      "pgm.add_node(daft.Node(\"tau1\", r\"$\\tau_1\", 0.25, 1.5, fixed=True))\n",
      "\n",
      "pgm.add_edge(\"tau0\", \"u\")\n",
      "pgm.add_edge(\"tau0\", \"v\")\n",
      "pgm.add_edge(\"u\", \"z\")\n",
      "pgm.add_edge(\"v\", \"z\")\n",
      "pgm.add_edge(\"z\", \"r\")\n",
      "pgm.add_edge(\"tau1\", \"r\")\n",
      "\n",
      "# [start_x, start_y, xlen, ylen]\n",
      "pgm.add_plate(daft.Plate([0.5, 0.5, 2, 2.5], label=r\"$i = 1, \\ldots, N$\"))\n",
      "pgm.add_plate(daft.Plate([1.5, 0.25, 2, 2.7], shift=0.5, label=r\"$j = 1, \\ldots, M$\"))\n",
      "\n",
      "pgm.render()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<matplotlib.axes.Axes at 0x111c94f10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAJFCAYAAAAfy3tvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XtclHX+///nIBCeAY9gWYllZSQC5iG1k2imYuah1Szb\nTERtNz+1hW3bbu1WatZuu5XndMvc1LASMyuBPFS7KuBpy1OiuSniioJngeH6/bFf+GUeghGu9zUz\nj/vtxh/AMNcTfPme51xzzXW5LMuyBAAAAKDKAkwHAAAAALwVZRoAAADwEGUaAAAA8BBlGgAAAPAQ\nZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBl\nGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUa\nAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoA\nAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGgAA\nAPAQZRoAAADwEGUaAAAA8BBlGgAAAPAQZRoAAADwEGUaAAAA8BBlGkCNmjlzpqKiohQQEFDpj549\ne5qODQBApQSaDgDAdw0ePFiLFy+Wy+WSy+Wq+LplWWd9/mOhoaGKi4uzKyIAAJeEPdMAakRKSooy\nMzM1c+ZMHTlyRG63u+IjNjZW2dnZZ32t/KOgoEATJ040HR8AgEpxWZZlmQ4BwLekp6crOTlZOTk5\natCgwVnfKywsVHh4uMrKygylAwCg+rBnGkC1mzBhgtLT088p0tL/ijaHcQAAfAXHTAOodpmZmect\n0pK0cOFCxcfH/+x9pKenKz09XVFRUSosLFRBQYEmTZpU3VEBALgklGkA1e5CRVqSFi9erJkzZ170\n51NTUzVz5kx9/vnnFV/bsGGD4uPjlZWVVW05AQC4VBzmAcA26enpkvSze6aTkpLOKdzt27dXq1at\nNGXKlBrLBwBAVVGmAdjm/fffl8vlUkxMzAVvM3PmTDVq1EhXXXXVOd/r0aMHZ/oAADgKZRqAbRYt\nWqRWrVpd9DYrVqy44G1atWqlwsJC7dmzpwbSAQBQdZRpALbIzc1VUVGRYmNjL3q7nJyci5bp8tsA\nAOAElGkAtkhNTZUkdejQ4aK327179wW/Fx4eLul/xRwAACegTAOwxcKFCyX977hnAAB8BWUaQI0r\nLCzUhg0bfvbNh5VVUFBQDakAALh0lGkANa78lHg/9+bDymrUqFG13A8AAJeKMg2gxrlcLoWFhSkl\nJcV0FAAAqhVXQARQ4wYOHKiBAwdW6ratWrXS4cOHz/u98q9X1x5uAAAuFXumAThKjx49tGvXrvN+\nr/wsHryJEQDgFJRpAI4SFxd3wVPf5ebmKioqSg0aNLA5FQAA50eZBuAoQ4YMkSRt2LDhnO+9//77\nGj16tN2RAAC4IMo0AEdp2LChZs2adc6bFdPT01VYWKjf/OY3hpIBAHAul2VZlukQAPBTGRkZWrFi\nhaKiolRYWKjDhw9r4sSJpmMBAHAWyjQAAADgIQ7zAAAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxE\nmQZgxIEDB/T888/rxhtv1OWXX65rrrlGw4YN0/r1601HAwCg0jibBwBbnTp1SiNHjtQHH3wgl8ul\n06dPV3wvICBAISEhuvLKK7V48WJdf/31BpMCAPDzKNMAbHPy5El17dpVW7duPatE/5TL5VK9evW0\natUqtW/f3saEAABUDWUagG0GDhyoTz755KJF+sfCw8P13XffKSwsrIaTAQDgGY6ZBmCLPXv2VKlI\nS9Lp06c1d+7cGkwFAMCloUwDsMWbb76psrKyKv3MyZMn9eqrr1b55wAAsAtlGoAtFixYoOLi4ir/\nXFFRkbZv314DiQAAuHSUaQC2OHr0qEc/FxgYqMOHD1dzGgAAqgdlGoAtgoKCPPo5y7IUEhJSzWkA\nAKgelGkAtrjhhhs8+rni4mJdffXV1ZwGAIDqQZkGYIvf/OY3ql+/fpV+JiAgQH369FF4eHgNpQIA\n4NJwnmkAtnC73WrevLkOHTpU6Z+pU6eOMjIy1KlTpxpMBgCA59gzDcAWtWrV0rRp06r0M/Hx8RRp\nAICjUaYB2CIvL0+DBw+WJNWtW/eib0isW7eurrvuOq1evVozZsywKyIAAFVGmQZQ4/Ly8hQZGSlJ\nKi0t1ZYtWzR27FjVq1dPDRo0UIMGDdSwYUOFhISoZ8+eWrJkibZu3aqHH35YycnJFGoAgGNxzDSA\nGvXTIl2rVq2K7506dUrr1q1TYWGh6tSpo7Zt21bcttzIkSM1Z84cTZ8+XaNHj7Y1OwAAPyfQdAAA\nvutiRVqSateurVtvvfWi9/HWW29JkpKTkyWJQg0AcBTKNIAa8XNFuioo1AAAp6JMA6h21Vmky1Go\nAQBORJkGUK1qokiXo1ADAJyGMg2g2tRkkS5HoQYAOAllGkC1sKNIl6NQAwCcgjIN4JLZWaTLUagB\nAE5AmQZwSUwU6XIUagCAaZRpAB4zWaTLUagBACZRpgF4xAlFuhyFGgBgCmUaQJU5qUiXo1ADAEyg\nTAOoEicW6XIUagCA3SjTACrNyUW6HIUaAGAnyjSASvGGIl2OQg0AsAtlGsDP8qYiXY5CDQCwA2Ua\nwEV5Y5EuR6EGANQ0yjSAC/LmIl2OQg0AqEmUaQDn5QtFuhyFGgBQUyjTAM7hS0W6HIUaAFATKNMA\nzuKLRbochRoAUN0o0wAq+HKRLkehBgBUJ8o0AEn+UaTLUagBANWFMg3Ar4p0OQo1AKA6UKYBP+eP\nRbochRoAcKko04Af8+ciXY5CDQC4FJRpwE9RpP9/FGoAgKco04Afokifi0INAPAEZRrwMxTpC6NQ\nAwCqijIN+BGK9M+jUAMAqoIyDfgJinTlUagBAJVFmQb8AEW66ijUAIDKoEwDPo4i7TkKNQDg51Cm\nAR9Gkb50FGoAwMVQpgEfRZGuPhRqAMCFUKYBH0SRrn4UagDA+VCmAR9Dka45FGoAwE9RpgEfQpGu\neRRqAMCPUaYBH0GRtg+FGgBQjjIN+ACKtP0o1AAAiTINeD2KtDkUagAAZRrwYhRp8yjUAODfKNOA\nl6JIOweFGgD8F2Ua8EIUaeehUAOAf6JMA16GIu1cFGoA8D+UacCLUKSdj0INAP6FMg14CYq096BQ\nA4D/oEwDXoAi7X0o1ADgHyjTgMNRpL0XhRoAfB9lGnAwirT3o1ADgG+jTAMORZH2HRRqAPBdlGnA\ngSjSvodCDQC+iTINOAxF2ndRqAHA91CmAQehSPs+CjUA+BbKNOAQFGn/QaEGAN9BmQYcgCLtfyjU\nAOAbKNOAYRRp/0WhBgDvR5kGDKJIg0INAN6NMg0YQpFGOQo1AHgvyjRgAEUaP0WhBgDvRJkGbEaR\nxoVQqAHA+1CmARtRpPFzKNQA4F0o04BNKNKoLAo1AHgPyjRgA4o0qopCDQDewWVZlmU6hDcKDw/X\nkSNHTMcA4CemT59OoQYAB6JMe8jlcok/HX4Oe6RRHVwulyQKNQA4UYDpAICvokijOj388MNKTk7W\njBkzTEcBAPwIx0wDNYAijerGMdQA4EyUaaCaUaRRUyjUAOA8lGmgGlGkUdMo1ADgLJRpoJpQpGEX\nCjUAOAdn8/AQZ/PAj1W2SHNKRQB2CAsL0+HDh03HAPwCZdpDlGmUq8oeaeYGnrjY3IwcOVJz5szh\ntHk4C2sNYB8O8wAuAYd2wDQO+QAAsyjTgIco0nAKCjUAmEOZBjxAkYbTUKgBwAzKNFBFFGk4FYUa\nAOxHmQaqgCINp6NQA4C9KNNAJVGk4S0o1ABgH8o0UAkUaXgbCjUA2IMyDfwMijS8FYUaAGoeZRq4\nCIo0vB2FGgBqFmUauACKNHwFhRoAag5lGjgPijR8DYUaAGoGZRr4CYo0fBWFGgCqH2Ua+BGKNHwd\nhRoAqhdlGvh/KNLwFxRqAKg+lGlAFGn4Hwo1AFQPyjT8HkUa/opCDQCXjjINv0aRhr+jUAPApaFM\nw29RpIH/oVADgOco0/BIQUGB9u7dq7y8PO3fv195eXk6ceKESktL5XK5FBgYqPr16ysyMlIRERGK\njIzUlVdeqQYNGpiOLokibYJlWfrhhx+0b9++irnJz8/XmTNnKv4NAgMD1ahRo3Pm5rLLLjMd3+c5\ntVAXFxfr+++/P2utOXTokEpKSuR2uxUYGKjg4GA1a9asYmYiIyN1xRVXKCAgwHR8AH7AZVmWZTqE\nN3K5XPKXP93p06e1evVqrVu3TtnZ2crOzlZRUZGuuuqqigeviIgI1atXT4GB/3t+VlpaqsLCQuXl\n5VV87NmzRxEREYqLi1NcXJw6deqkzp07V/yMXUwWaX+am0OHDikzM1PZ2dnKyspSTk6OQkJCdMUV\nV1TMTbNmzRQSEqLAwECVlZWpuLhYBQUFFcWp/OP666+vmJuuXbvqhhtukMvlMv0r2sbOuRk5cqTm\nzJmj6dOnGynUO3bs0OrVqyvWmn//+99q3ry5WrRooYiICEVERKhx48YKDg5WrVq15Ha7dfr0aeXn\n51fMzQ8//KDjx4+rffv2iouLU3x8vG6//XY1b97c9t/HFH9aawDTKNMe8vWF6uDBg1q2bJmWLl2q\njIwMRUdHq0uXLhWFJioqqsplxu12a9u2bRUPkl9++aX27Nmj3r17KzExUXfddVeN77k2vUfa1+dm\n27ZtSktL09KlS7V582bdeuutuvnmmyvmpmnTplW+z5MnT2rTpk0Vc5OZmanAwEAlJiYqMTFRXbt2\nVVBQUA38Ns5h99zYWajdbrf++c9/Ki0tTWlpaTp27JjuuOMOxcfHKz4+XjExMapbt26V7/fQoUPK\nyclRdna21q1bp5UrV6pNmzYVc9O2bVuffkLm62sN4CSUaQ/54kJlWZZWrlypN998U+np6UpISFBi\nYqJ69+6txo0b18g2f/jhB3388cdKS0vT119/rYEDB2rcuHGKjY2t9m2ZLtKSb87N6dOntWjRIk2d\nOlX/+c9/1L9/f/Xr10+33367QkJCqn17lmVp8+bNFeVr3759SkpK0qhRo9SiRYtq354TmJibmi7U\n+fn5euuttzRjxgyFhYVVzE1sbGyNHJ5RXFys1atXKy0tTUuWLFGjRo00duxYDR061KOy7nS+uNYA\njmXBI770pztz5ow1bdo067rrrrPatm1rvfnmm1ZRUZHtOfLz862XXnrJatmypdWxY0drwYIFltvt\nrpb73r9/vyXJkmSVlpZWy316wpfm5uDBg1ZKSorVpEkTq1evXtaSJUuM/G23bNlijR071goLC7MG\nDhxorVu3zvYMNc3U3Dz88MOWJGv69OnVdp+bNm2yhg4daoWGhlqPPPKIlZ2dXW33XVlut9v69NNP\nrX79+lnh4eHW+PHjrX379tmeoyb50loDOB3/2zzkCwuV2+225s+fb7Vq1crq1auXtXLlSqusrMx0\nLKu0tNRasmSJdfPNN1sxMTHW8uXLLymXU4q0ZfnG3Bw9etR67rnnrEaNGlnjxo2zdu7caTqSZVn/\ny/X6669bLVq0sAYNGmRt27bNdKRqY3JuqqtQ5+bmWsOHD7eaNWtmTZkyxTp8+HA1Jbw0e/bssZ54\n4gkrPDzcmjBhgmNyXSpfWGsAb8H/Ng95+0K1cuVKq127dlbHjh2tL774wnSc8yorK7MWL15stWnT\nxrrtttuszZs3V/pnS0pKrJKSEkcVacvy7rlxu93WtGnTrGbNmlnDhw+3cnNzTUc6rxMnTliTJk2y\nGjdubCUlJVkFBQWmI10y03Pz00J9+vTpSv9sUVGR9dhjj1mNGjWynnvuOevo0aM1FfOS/Oc//7Ee\neeQRq3Hjxtarr77qiPXiUpieGcCf8L/NQ966UB07dsx69NFHrRYtWlipqamO2BP9c0pKSqzp06db\nTZo0sV544QWruLj4orcvKyuz+vTpY7Vu3dpRRdqyvHdudu3aZd16661W586drY0bN5qOUymHDx+2\nfvWrX1mRkZHWkiVLTMe5JE6Ym/JCHRMTY4WGhloHDx782Z/5/PPPrZYtW1ojR46s1O2dYOvWrdYd\nd9xhdezY0dq6davpOB5zwswA/oL/bR7yxoVq5cqVVqtWrawHH3zQK1/K3Lt3r9WrVy8rNjbW2rJl\nywVvt2DBAqtOnTqWy+WyJDnqd/W2uSkrK7PeeOMNq3HjxtYrr7zimCclVbFq1SorKirKGj58uKNm\noSqcMjd16tSxJFlBQUFWnz59Lni7o0ePWklJSVbLli2tzz77zMaE1cPtdltTp061GjdubE2ePNkr\n594pMwP4A85o7wcsy9Krr76qX/ziF/rrX/+qt99+W2FhYaZjVdkVV1yh5cuXa+zYsbr99tu1cOHC\nc25z8OBBJSUl6eTJk7IsS7Vq1VJaWpqBtN7vxIkTuu+++/T222/ryy+/1BNPPOGVF7fp3r27Nm3a\npNDQUMXHx+vf//636UheqaCgQCUlJXK5XCopKdEXX3yhRYsWnXO77777Th07dlRJSYk2b96snj17\nGkh7aQICAjRmzBitX79ey5YtU2JiooqKikzHAuBUptu8t/KWP92pU6esBx980IqJibG+//5703Gq\nzcaNG60rr7zSeuaZZyrO+FFWVmb17t3bCgoKsoKCgqwGDRpY7733nqMOZfGWudmzZ48VExNjjRgx\nwjp16pTpONVm3rx5VuPGjb3usA+nzM3WrVut6Ojoij3UDRo0OOvwjRUrVlhNmza1pk2bZjBl9Sou\nLrYeffRR67rrrrN27NhhOk6lOWVmAH/A/zYPecNCdfDgQatjx47W4MGDrePHj5uOU+3y8/Otbt26\nWf3797dOnjxpLViwwHK5XFZISIjVu3dvKz8/33TEc3jD3Pzzn/+0mjdvbr366quOeiJSXdauXWu1\naNHCmjx5sukoleakuSktLbUmTpxo1a5d23K5XBWHe0ydOtVq3ry5tXLlSsMJa8aMGTOspk2bWpmZ\nmaajVIqTZgbwdVy0xUNOPyH+gQMHdOedd6p///568cUXffZKX8XFxfrlL3+pvXv3auPGjbIsS7Nn\nz9Z9993nyN/Z6XOzevVqDRo0SHPmzFHfvn1Nx6kx+/btU8+ePTVgwAD96U9/cuSs/JgT52bbtm0a\nMmSItmzZomHDhmnt2rVasWKFrr76atPRaswXX3yh++67T++8847uuusu03EuyokzA/gqyrSHnLxQ\n5efn67bbbtOwYcP07LPPmo5T49xut+6880598803Wr9+va666irTkS7IyXOzZs0a3XvvvVqwYIHu\nvPNO03Fq3H//+1/16NFDffv21Ysvvmg6zkU5dW7cbrf69eunr776St9++63PXoHyx/75z3+qf//+\nji/UTp0ZwBfxBkQfU1hYqISEBN13331+UaQlqVatWsrMzFTv3r2VlJSkkpIS05G8TlZWlgYOHKj3\n3nvPL4q0JDVp0kTp6elasmSJXnrpJdNxvNK0adO0Y8cOvynSktS5c2ctWbJEDz74oFatWmU6DgAH\nYM+0h5z4rN/tdqtPnz5q3bq1Xn/9dce/dF3dSktLNWDAAF1xxRWaOnWq6Tjn5cS52b9/vzp27KjX\nX39d99xzj+k4tsvLy1PHjh31l7/8RQMHDjQd57ycODefffaZHnroIX399dc+fWjHhaSnp2v48OH6\n+uuv1apVK9NxzuHEmQF8FXumfUhKSopKS0v12muv+V2RlqTAwEDNnz9fK1eu1PTp003H8QqnT5/W\ngAEDlJyc7JdFWpIiIiL00UcfKTk5WZs2bTIdxyvs2LFDDzzwgBYtWuSXRVqSevTooWeeeUb9+/fX\nsWPHTMcBYBB7pj3ktGf98+bN0/PPP6+1a9eqUaNGpuMY9d133+mWW27RokWLdOutt5qOcxYnzY1l\nWXrooYd05swZvffee375BOzHFi5cqJSUFK1fv15NmjQxHecsTpqboqIidezYUU888YRGjRplOo5R\nlmUpKSlJhw4d0uLFixUQ4Jz9U06aGcDXUaY95KSFKjc3VzfffLNWrlypG2+80XQcR/jss8/0yCOP\naMuWLQoNDTUdp4KT5mb+/PmaOHGi1q1bpzp16piO4whPPfWUvvvuOy1evNhRTy6cNDcPPfSQgoOD\nNXPmTNNRHKG4uFjdunXTiBEjNHbsWNNxKjhpZgBfR5n2kFMWqrKyMt1xxx3q16+fnnjiCdNxHGXM\nmDE6c+aM5syZYzpKBafMTV5enmJiYvTJJ58oLi7OdBzHOHPmjGJjY/W73/1OQ4cONR2nglPmZtmy\nZfrVr36lzZs3q169eqbjOMbWrVvVrVs3rV+/3jGHvThlZgB/QJn2kFMWqjfeeEP/+Mc/tGbNGq+8\n1HNNOnbsmKKjozVt2jT17t3bdBxJzpgby7J0zz33KDo6Wi+88ILRLE60fv169e3bV5s3b1azZs1M\nx5HkjLk5cuSIoqOjNW/ePN1+++1GszjRyy+/rOXLlysjI8MRh3s4YWYAf0GZ9pATFqp9+/apXbt2\n+uqrr9SmTRujWZwqIyNDv/zlL7Vt2zZHHMrghLlJTU3V888/r6ysLF122WVGszjVb3/7W+Xm5mrB\nggWmo0hyxtyMHj1atWrVcuyZckxzu9265ZZbNGrUKI0cOdJ0HEfMDOAvKNMecsJCNWrUKDVq1EiT\nJk0ymsPpBg0apA4dOiglJcV0FONzU1JSohtuuEHTp0/3m/NJe+LkyZO65pprtGTJEsXHx5uOY3xu\ntm3bpm7dumnHjh0KCwszlsPp1q1bp3vvvVc7duww/uTd9MwA/oQy7SHTCxUPbpW3fft2de3a1RF/\nK9NzM336dC1evFgrVqwwlsFbzJgxQ6mpqY74W5mem4EDB6pjx4566qmnjGXwFk558m56ZgB/Qpn2\nkOmFyikLtrdwyl58k3PjtL2tTldSUqK2bdtq6tSp6tGjh9EsJufGSXtbvYFTdnSYfowC/All2kMm\nF6qtW7fq9ttvV25uLg9ulfTDDz8oOjpau3fvNnqqPJNz88YbbygzM1MffPCBke17o/fee08zZszQ\nypUrjeYwOTf9+/fXXXfdpTFjxhjZvjd6+OGH1apVK/3ud78zloEyDdiHMu0hkwvVr3/9azVo0IAz\nMVTR0KFD1blzZ/361782lsHU3FiWpbZt22ratGmOu5CNk5WUlKhly5ZKT09X27ZtjeUwNTd79+5V\n+/bttXfvXtWtW9f27XurjRs3KjExUbm5uQoMDDSSgTIN2Mf8+XtQJSdOnND8+fOVlJRkOorXGTt2\nrKZOneqXDzCrVq2Sy+VS9+7dTUfxKkFBQRo1apTfXp5+5syZGj58OEW6imJiYnT55Zdr2bJlpqMA\nsAFl2sv84x//ULdu3dSyZUvTUbxO165dFRQUpMzMTNNRbDd16lSNHTvWUVf18xZJSUmaP3++jh07\nZjqKrYqLizV79mwO7/DQ2LFj9eabb5qOAcAGlGkvM3/+fD388MOmY3gll8ulkSNH6h//+IfpKLY6\nduyYli9fruHDh5uO4pUuv/xyde7c2e/2MqanpysqKkrXXXed6SheadCgQVq3bp3y8/NNRwFQwyjT\nXqSgoEA5OTlKSEgwHcVrJSYm6uOPP1ZZWZnpKLb5/PPP1aVLFzVs2NB0FK/Vv39/paWlmY5hq7S0\nNN1zzz2mY3itkJAQ9ezZ0++ehAH+iDLtRZYvX6477rhDtWvXNh3Fa7Vq1UpNmzbVunXrTEexzdKl\nS5WYmGg6hlfr27evPv30U5WUlJiOYgvLspibapCYmKilS5eajgGghlGmvUhaWhoPbtWgX79+WrJk\niekYtnC73Vq2bJn69u1r63ZzcnIUFxengIAAxcfHKyMj45zbpKSknPfrThQZGanWrVtr9erVpqPY\nIicnR/Xq1VObNm1MR/FqvXv3VkZGhk6dOmU6CoAaRJn2EpZlaeXKlbYd4pGSkqKAgIBzPsLDw8+5\nbWpq6nlvGxDgzPHq2bOnVq1aZTqGLb755huFh4fryiuvtG2bqampio+P14YNG+RyuSoOTZoyZYok\nqbCwUCkpKdq4caNXXdLcn+bGzrUmNzf3gutH+cfGjRvP+bkLrVGzZs2yJXdlNGrUSNdff72ysrJM\nRwFQg8ycABNV9sMPPyggIECXX365Ldu7+eab1aNHD2VnZ+vIkSMVXz/f2SCioqLUo0cP7d69W7t2\n7brobZ0gNjZWmzdvVmlpqbFzwNolOztbHTp0sG17hYWFSkpK0uDBg5WQkKDw8HAdPnxYu3bt0ooV\nK7Rw4ULt2rVLjRs3VnZ2tm25qkNcXJxmz55tOoYtsrOz1atXL1u21ahRI7388suSpJdeekmFhYUV\n3xs8eLA6dOigq6+++pyf69mzp4qKijRz5kxJUkJCglq1amXrvFdGfHy8srOz1a1bN9NRANQQLtri\nIbtPiP/hhx9q9uzZtr+Zpaio6KxL4oaFhamgoOCCt//x3miXyyW3212j+TzVpk0bpaamKjo62tbt\n2j0348aNU+vWrfV///d/tmwvIyNDRUVFuvfee8/7/fT0dI0ZM0bZ2dlq0KCBLZmqy969e9WhQwcd\nOHDA9ieKds/Ntddeqw8++EA33nijbduUpFmzZmn06NEVn/fo0UOff/75BW9fWFio8PBwDR48WAsX\nLrQjYpXNmTNHmZmZevfdd23dLhdtAezjzNfhcY7s7GzFxcXZvt2qngHCW84YERcX53V7Rj1h99zc\neeedFyzSOTk5GjJkiFasWOF1RVqSrrjiCpWVlWnfvn2mo9SooqIi7d+/38gp8UaNGqXQ0NCKz9PT\n07V79+4L3r58r/TTTz9d49k85S9rDeDPKNNeYuvWrbbvJbJDQkKChgwZYvt2b7zxRn377be2b9du\nTpmbnJwcJSUlKScnR1ddddVFb2tqJn6Oy+XSjTfeqK1bt5qOUqO2b9+uNm3aGDsE6qfFOCUl5YK3\nnThxohINq0uEAAAgAElEQVQSEhQTE1PTsTx2ww036LvvvlNpaanpKABqCGXaS+Tl5alFixamY1Sr\nwsJCZWdn6+abb7Z92y1atFBeXp7t27XTiRMnVFxcfNZhOiakp6drwoQJyszM/NkiXdWZiIuLU3Jy\ncjWkrBx/mBvTa01SUtJZn6empqqoqOic282cOVNFRUUXLdsXYufcBAUFqVGjRjp48KAt2wNgP8q0\nl9i/f78iIiJMx6hWoaGhOnz4sH7zm9/Yvu2IiAi/KEURERFG3wianp6umTNn6vPPP6/UoR1VmYnc\n3Fzt3r1b8fHx1RG1UiIiIrR//37btmeC6bWmYcOGeuqpp8762sSJE8+53eTJkxUVFaU77rijSvdv\nam58fb0B/Bll2gtYlqUDBw74XJk2KTIy0i9KUWRkpLHtp6amKj09XYsWLaqR+2/VqpUOHz6sRx55\npEbu/3wiIyN9vhTl5eUZnRvp3EM9Xn755bP2TpcfSz158uQq37epufH19QbwZ5RpL3D8+HEFBQVx\n5cNq1KRJE/33v/81HeOS7dmzR++8847+/e9/n3NM5qFDh9SkSRMjuWbOnKnCwkJNmjTJyPZryvnm\nxrIs7du3T2lpaV5zEZqLMTk35Ro2bKhBgwad9bXyNxtK/9srHRYWdsE3uzqNr6w3AM7Pt0+y6yOK\ni4sVHBxsOoZPCQ4OVklJyVnntLVLdW7znXfe0QsvvKDLLrtMxcXFat26tbp27aouXbro4MGDCgoK\nqrZtVdasWbPUuHHjixadwYMH6+WXXz7v+YOdLDg4WEePHlVaWprWrl2rVatWacuWLSouLlZAQIBa\ntGhRY5eqt2tWjx8/7oj1ZvLkyUpNTa34fOLEiXryySeVk5OjjIyMinNTe4Py9QaAb6JMewG32+3Y\nqwn+1PneKPRTs2bN0q5du5STk6PJkyerffv2NiQ7W61atc45h7ZdamKb5Q/U3377rb799tuKvXh2\nn0c7NTVV4eHhFy3Subm52rBhw1lFurIz8eO93Tk5OZoxY4athXzfvn1atmzZBc/3vnPnzhqbKTtn\n1e1223oYxPlcffXV6tGjh9LT0yX9799+1qxZ+vzzz+Vyuar0XgvTc1OrVi3HnnMfwKXzjobm5wID\nA71iIa7MnrPU1FRFRUVp0qRJuvrqqzVq1Cgbkp3L7XYrNDRUlmXZ+iGpWu/vz3/+s4KDgxUSEqIG\nDRrolltu0dNPP62PPvpIb775pq6//nrb/qblp7+bOHGievbsed5DHnJzc5WQkHDWXsXKzkR5ISr/\nCA0NPesCH3Zo0aKFBgwYoFWrVunVV19V//79FRkZqaCgINWpU0c33nijV8zNxT4eeeQR3Xbbbbb+\nXS/kp8dEp6SkaPHixee8QfFinDA3brfb56+2Cvgz/nd7gaCgIBUXFxvbfmhoaKWK8vnecf9TWVlZ\nFXuIsrKyjO1xP3PmjJFDIKrbL37xC7Vu3VqxsbGKjIw868wdH3zwgc6cOWNLjvLLiGdmZiomJkaL\nFy9WQkKCYmNjNXr0aIWFhWnFihWaNWuW4uLiztpzXdmZmDRpkn77299WfH7kyBHt2bOnRn+vnzpz\n5oxq166t7t27q3v37nr88cclSceOHdOGDRt02WWX2ZqnJpheb36sffv2io2NVU5OjqT/zZnL5arS\nRVqcMje+sN4AOD/2THuBevXqye126+TJk0a2/+MLaBQWFp73UI709HQtXrxYsbGxZ339x1cv2717\ntxISEiT9/y/133fffTWU+uIOHjyopk2bGtl2dYqIiFC/fv3UokWLc06B17RpU9vObTtkyBClpqZW\nXDxj4MCBmjFjhnJycjR69GgNGTJEs2bNUlRU1Fl7rCs7E0VFRWrUqNFZp9fLyso6Z95q2oXmpn79\n+urevbs6duxoa56aYOfcVMZPi3NSUlKlr6Dp9LkB4Bso017A5XIZPU/p5MmT1apVK0n/e6l58ODB\nysjIUE5OjlJTUzV69Gj16tVLqampZx3XaVmWYmNjFRUVpVmzZunqq6/WnXfeKUmaMWOGpHMv0GAX\nJ5z+q6bZdRq3jIwMDRky5JwLsowaNari5fiwsDClpKRo586dZxWbys5Ew4YN9eSTT1Z8npOTo6Ki\noooibpfyc3f7MqedE3ngwIEV64/L5arSRVqcNDe+vt4A/ozDPLxE+cUioqKibN92w4YN9d1332nx\n4sWaMWNGxXGv0v8OAbnvvvt05MgRNWjQQC6X66w9pEePHtXRo0eVm5t71n3OnDlTCQkJld7DVN1M\nX5jCDuWlyLKsGr1wS3kZPp/y41QroyozUf6mtB49elQuZDXZv3+/2rVrZ+s27ebEC9NMnjxZgwcP\nVlJS0s9eRfNiTM6Nr683gD+jTHuJyMhI7du3z2iGgQMHauDAgRe9zeeff/6z95Oenu7xZYCry759\n+3z+wa127dqqXbu2CgoK1LhxY9NxLqqqM7FixQqFhYVdUrHyhD/MjRPWmp8aOHCgysrKLvl+TMzN\nmTNnVFhYaPzc3QBqDod5eIkbbrhBW7ZsMR2jWsyYMUNhYWEVlwGeMmWK7Rm2bNmitm3b2r5du3nL\n3FR1JjIyMmzfu2hZll/MzXXXXacdO3b45HmRTczNN998o2uvvVa1atWydbsA7EOZ9hJxcXHKzs42\nHaNaLF68+Kw3NRYUFNieITs7W3FxcbZv127eMjdVmYnyMzvYfdzrnj17FBIS4vN7puvVq6eWLVvq\n22+/NR2lWpmaG39ZawB/Rpn2EvHx8crKyqo456y3K39AmzJlipKTk23d9pEjR5Sfn682bdrYul0T\nyufGG1R2Jkwd95qdna34+Hhbt2mKN81NZZmam6ysLL+ZG8BfUaa9RPmFIb7//nvTUS7ZjBkzNH36\ndE2YMEFxcXG2H/eanZ2tmJgYv3jZNS4uzitKUVVmwtTx0llZWX6zh9Fb5qYqmBsANcVl+cquTpu5\nXC7b9xLff//96t69u+1X7/I1Tz75pGrXrq0//vGPtm/b7rkpKytTZGSkvvrqKyNngqkJAQEBGjx4\nsBYuXGjrdtu3b6+//e1v6tatm63bleyfm02bNmnAgAHatWtXjZ4Jxk4m5ubgwYO69tprlZ+fb/sF\nfUw8RgH+ij3TXiQxMVFpaWmmY3i9pUuXKjEx0XQMWwQEBKhv375aunSp6SjVovylersv9rN37179\n5z//UefOnW3drik33XST3G63vvnmG9NRqoWpuVm2bJkSEhJ84sqYAC6MMu1F7rrrLq1Zs0YnTpww\nHcVr7dixQ8eOHbP9CmgmefOTsJSUFM2aNavi8xkzZigqKuqsy5Hb4eOPP1afPn0UGOgfZxN1uVxK\nTEz02idhTpmbtLQ0v3niDvgzyrQXadiwoTp27KhPP/3UdBSvtWTJEvXt21cBAf4z+j169FBWVpaR\ns6ZcipycHE2ZMqXiLAypqanKzMzUihUrbM/y0UcfqV+/frZv16TExER9+OGHpmNUmVPm5uTJk8rM\nzNTdd99t63YB2M9/GoWPGD58uN566y3TMbySZVmaPXu2hg8fbjqKrerUqaPExES98847pqNUSWxs\nrAYNGqTQ0FAlJycrIyNDu3fvtv0NZLt371ZOTo7flaLbb79d+/fv1+bNm01HqRKnzM2CBQvUrVs3\nNWrUyNbtArAfb0D0kKk3d5w6dUotW7bUv/71L595Q5ld0tPT9fjjj2vTpk3G3lRlam6+/vprjRgx\nQtu3b/ervfLVYcKECSopKdGrr75qLIOpufnjH/+ovLw8TZs2zfZtezPLshQfH68XXnhBvXv3NpKB\nNyAC9uFR1cvUrl1bDz30kGbMmGE6iteZOnWqxo0b5zNnJ6iKzp07q169ehVvxELlnD59WnPmzLH9\nXOhO8cgjj2jBggU6evSo6SheZd26dSosLFSvXr1MRwFgA8q0F0pOTtbcuXN17Ngx01G8xu7du7Vy\n5Urdf//9pqMY4XK5NHbsWL322mumo3iV9957T+3bt9c111xjOooRkZGRSkhI4NCyKvrrX/+q5ORk\nXgUC/ASHeXjI9Etow4cP1zXXXKM//OEPxjJ4kwceeEBRUVF67rnnjOYwOTenT59WmzZtNH/+fHXt\n2tVIBm9y5swZtWnTRu+++67xv5fJudm0aZN69eqlnTt3qn79+kYyeJPNmzerZ8+exv9eph+jAH9C\nmfaQ6YVq9+7dio+P19atW9W0aVNjObyBUx7cJPNz8/bbb2vWrFlas2aNXx7uUhWvvfaaMjIyHHF6\nONNzw5P3yuvTp4969eqlX//610ZzmJ4ZwJ9Qpj3khIXq17/+tVwul/76178azeF0Tnlwk8zPjdvt\nVrt27TRp0iT17dvXWA6nO3r0qK655hqlp6crOjradBzjc5Obm6sOHTrw5P1nrF69WiNGjNC2bduM\nX6jF9MwA/oQy7SEnLFT5+fmKjo7WihUr1K5dO6NZnGrp0qUaP368vv32W+MPbpIz5mbZsmUaP368\nNm3apDp16hjN4lSPPfaYjh49qrlz55qOIskZczN+/HgVFhbq73//u9EcTlVcXKwOHTooJSVFw4YN\nMx3HETMD+AvKtIecslDNnTtXf/vb37Ru3ToFBQWZjuMohw8fVnR0tObPn6/bbrvNdBxJzpmbYcOG\nqXnz5vrzn/9sOorjrF69WkOHDtWWLVsUHh5uOo4kZ8zN8ePHFR0drTfeeEN9+vQxmsWJnnvuOWVl\nZWnp0qWOOITKCTMD+AvKtIecslBZlqU+ffqoU6dO+v3vf286jqM8+OCDatiwoV5//XXTUSo4ZW4K\nCgoUHR2tRYsWGX9znZOcOHFC7dq105///GdHXQbaKXPzxRdf6IEHHtCWLVsUFhZmOo5jbNy4UQkJ\nCdq4caNatGhhOo4k58wM4A8o0x5y0kL1ww8/qH379lq+fLni4+NNx3GEDz74QE8++aQ2bdqkevXq\nmY5TwUlz89FHH+nJJ59UVlaWGjZsaDqOI4wdO1ZHjx7Vu+++azrKWZw0N+PGjVNhYaHeffddR+yB\nNe3kyZPq1KmTHn/8cT300EOm41Rw0swAvo4y7SGnLVQffPCBxo8fr3Xr1ql58+am4xi1ZcsW3XHH\nHVq2bJluvvlm03HO4rS5GTdunL7//nstWbJEtWrVMh3HqNmzZ2vKlClau3atQkNDTcc5i5Pm5sSJ\nE+rataseeOABPf7446bjGGVZloYOHarAwEDNmzfPUU8unDQzgK+jTHvIiQvV888/r08//VQrV650\nxJvtTDh06JBuvvlm/fGPf9Tw4cNNxzmH0+ampKREPXv2VMeOHTVp0iTTcYz58ssvde+992rNmjVq\n06aN6TjncNrcfP/99+rUqZPmzp2ru+66y3QcY1566SV99NFHWrVqlWrXrm06zlmcNjOAL6NMe8iJ\nC1VZWZmGDBmiunXr6u9//7uj9pLY4cyZM+rdu7c6dOigyZMnm45zXk6cm/InIM8//7weeOAB03Fs\nt2fPHnXp0kVz58517OWfnTg35U9AVq1apeuvv950HNstWbJE48aN09q1ax1znPSPOXFmAF/FtU59\nSEBAgN5++21t27ZN//d//+dXC2lJSYl+8YtfKDQ0VC+99JLpOF6lcePGWrp0qZ588kl99NFHpuPY\nat++fbrzzjv1zDPPOLZIO1XXrl31yiuvqGfPntq1a5fpOLZKT0/XqFGj9OGHHzqySAOwF2Xax9St\nW1effvqpVq9erSeffNIvCnVxcbGGDRum4uJiLViwwO+P/fVE27ZttWzZMo0ePVpLliwxHccWP/zw\ng+644w6NHj1a48aNMx3HKz344IN65plndOedd/pNoU5PT9fQoUO1ePFidejQwXQcAA5AmfZBYWFh\nSk9P1+rVqzVu3DiVlpaajlRjTp48qUGDBqm4uFgffPCBgoODTUfyWnFxcRWF+r333jMdp0bt2rVL\n3bt316hRo/TUU0+ZjuPVkpOT9fTTT+u2227TN998YzpOjUpLS9OwYcP0wQcfqFu3bqbjAHAIyrSP\nCg8P14oVK7Rz50717dtXhYWFpiNVu71796pr164KCwtTamqq377psjrFx8drxYoVevrpp/Xss8+q\nrKzMdKRql5GRoS5duiglJUW/+c1vTMfxCaNHj9bkyZN1++23a+nSpabjVDvLsjR58mSNGTNGH3/8\nMUUawFko0z6sYcOGWr58ua677jp17NhR27dvNx2p2nz11Vfq1KmThg8frr///e9c/bEaRUdHa926\ndfriiy9077336tixY6YjVQvLsvT666/r/vvv14IFCzR69GjTkXzKsGHDtHTpUo0ZM0aTJk3ymUPM\nTp06peHDh+v999/X2rVrHXe6TQDmUaZ9XGBgoF577TU99dRT6tatm+bNm+fVD3Jut1uvvPKKBgwY\noDlz5ujxxx/3u7OW2KFp06bKzMxUkyZN1KlTJ2VnZ5uOdEmOHDmiBx54QDNnztTXX3+t22+/3XQk\nn9SxY0etXbtWixcv1sCBA3Xw4EHTkS7JN998o65du8qyLK1Zs0aXX3656UgAHIgy7SdGjhypTz/9\nVFOmTNE999yjvLw805GqbPv27erWrZuWLVumf/3rX359fls7BAcHa+bMmXrmmWd0991369lnn9WZ\nM2dMx6qyjz/+WNHR0QoPD9e//vUvtWrVynQkn9aiRQutWbNG1157rW666SYtWrTIdKQqKy0t1cSJ\nE3XbbbcpOTlZ8+fPd9x5pAE4B2Xaj8TGxiorK0s33XSTYmJiNHv2bK94c+KpU6c0adIkde3aVfff\nf78yMjIoRDZxuVwaNmyYNm3apC1btig+Pl5r1qwxHatS9u/frwcffFCPPfaY5s+fr7/97W+qW7eu\n6Vh+ISQkRJMmTdKSJUv0hz/8QYMGDdL3339vOlalrF+/Xp07d9YXX3yh7OxsjRo1ile/AFwUZdrP\nBAcH609/+pOWL1+ud955RzfddJM++ugjRx76UVpaqrfeekvXXnut1q5dq7Vr12rcuHEKCGBs7da8\neXN9+OGHeuaZZ/TAAw+ob9++2rJli+lY53XkyBFNmDBB0dHRat68uTZv3qxbb73VdCy/1LFjR23Y\nsEFt27ZVbGysxo8fr//+97+mY53X9u3bNWjQIA0YMEBjxozRZ599ppYtW5qOBcAL0Er8VGxsrFat\nWqVXXnlFf/jDH9SlSxctWbJEbrfbdDSdPn1a8+bN00033aR58+bp/fff14cffsjeaMNcLpd+8Ytf\naPv27UpISFBCQoKGDx+unJwc09EkSQcPHtQLL7yga6+9VgUFBdq0aZNefvll9kYbFhISoueff17f\nfvut3G63rr/+ej377LPat2+f6WiSpC1btuiRRx5R165d1aFDB+3YsUMPP/wwe6MBVBpl2o+5XC7d\nfffd2rBhgx577DFNmjRJrVq10osvvqj8/Hzb8+Tm5iolJUUtW7bU/Pnz9Ze//EVffPGFOnXqZHsW\nXNhll12mxx57TDt37tSNN96oAQMGqFOnTnrnnXd0+vRpW7NYlqUvv/xSw4YNU5s2bbRnzx6tWbNG\ns2bN4s1iDtOsWTO9/vrrWrdunQoKChQdHa1BgwYpMzPT9lfGyi/w1L17d91111264oortGPHDqWk\npKhOnTq2ZgHg/VyWE1/f9wIul8uRh0ZcqpycHE2bNk3vv/++4uPjlZiYqH79+unqq6+u9m1ZlqVv\nv/1WS5cuVVpamnbs2KGHHnpIycnJat26dbVvzwl8cW7cbrc++eQTTZ06VWvXrlWvXr2UmJio3r17\nKzQ0tNq3V1paqi+//LJibgICAjRmzBiNGDFCYWFh1b49J/DFuTl27JjeffddTZ06VceOHVNiYqIS\nExPVvXv3Grn40rFjx/TZZ58pLS1Nn3zyidq1a6exY8cqMTHRJ0+t6YszAzgVZdpDvr5QnThxQitW\nrFBaWpo+/vhjNW3aVF26dFFcXJzi4uIUHR1d5YuknDhxQhs3blR2drays7P15ZdfqrS0tOJB9NZb\nb/X5Kxj6+tzk5eXp448/VlpamlatWqWYmBh16NChYm6uueaaKh/zfvDgwYqZyc7O1urVq3X11VdX\nzE27du18/iV5X54by7L0zTffKC0tTWlpaRVn7Smfmbi4OEVERFT5Pnft2lUxM1lZWcrKytItt9yi\nxMRE9e3bV1dccUUN/UbO4MszAzgNZdpD/rRQud1uZWVlaf369RUPTjt37lSTJk0UERGhyMhIRURE\nqF69egoMDJT0v72HRUVF2r9/v/Ly8pSXl6cjR46obdu2FQ+QHTt2VHR0tM8XoR/zp7k5ceKEvvrq\nK2VlZVXMzX//+9+KeYmMjFSzZs0UEhKi5cuXKyEhQS6XS4cOHVJeXl7F7JSWlp5VrG655Ra/O4TD\nn+bmwIED+vLLLytKcPk5zsvnJiIiQk2aNFFQUJBq1aolt9ut06dPKz8/v2Ju9u/fr7CwsLPmplu3\nbqpfv77h384+/jQzgGmUaQ/5+0J1+vTpipJc/nHixAmVlJQoICBAgYGBql+//lllu3nz5j75cmpV\n+PvcHD169KyinJ+frzNnzujpp5/WsGHDFBMTo0aNGlWUpoiICDVt2tSvnnCdjz/PjWVZFU+wymfn\n0KFDKi0tldvtVmBgoIKDg9WsWbOKmYmMjFTDhg1NRzfKn2cGsBtl2kMsVPAEc3N+LpdL77//vgYN\nGmQ6iiMxN6gqZgawD2fzAAAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZ\nBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkG\nAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYA\nAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAA\nADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAA\nPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8RJkGAAAAPESZBgAAADxEmQYAAAA8\nRJkGYMz27dvVqFEjSdLgwYP15z//2XAiAACqhjINwJiwsDAdO3ZMkhQSEqLatWsbTgQAQNVQpgEY\n07RpU9WvX1+SFBwcrLi4OMOJAACoGso0AKNiYmIkSSdPntRNN91kOA0AAFVDmQZg1K233qqAgABd\nccUVCgkJMR0HAIAqoUwDMKpDhw4qKytTp06dTEcBAKDKKNMAjCo/Trpbt26GkwAAUHUuy7Is0yG8\nkcvlEn86VJW/zs3Ro0f14YcfauXKldq7d68OHDiggoICHT9+XGVlZTpz5oyCgoJUq1Yt1a5dW6Gh\noWrWrJlatGih2NhYDRw4UNdcc43pX8MYf50beI6ZAexDmfYQCxU84S9zs3fvXv3lL3/RJ598or17\n9+r06dMKDg5WeHi4wsLCFB4ermbNmqlZs2aqXbu2atWqJbfbrZKSEhUUFOjAgQM6dOiQCgsLdfjw\nYR0/fly1atVS06ZNdfPNN2v8+PG67bbbTP+atvGXuUH1YWYA+1CmPcRCBU/48txs375dzzzzjDIz\nM3XkyBE1bNhQbdu2VVxcnDp16qQGDRp4fN+lpaXauHGj1q9fr02bNmn//v0KCgpSbGysnnjiCQ0a\nNKgafxPn8eW5Qc1gZgD7UKY9xEIFT/ja3JSVlWnOnDl68cUXtWfPHjVv3lxdunTRPffco7CwsBrb\nbmlpqTIyMpSenq4dO3aoTp06euCBBzRp0qRLKu1O5Wtzg5rHzAD2oUx7iIUKnvCVuSkrK9Mf/vAH\nvfrqqyouLlZsbKxGjhypFi1a2J6luLhY//jHP7RixQodO3ZMd955p+bPn6+mTZvanqWm+MrcwD7M\nDGAfyrSHWKjgCV+Ym6lTpyolJUVnzpxR//79df/99yswMNB0LEnSunXrNG3aNB0+fFhDhgzRW2+9\npTp16piOdcl8YW5gL2YGsA9l2kMsVPCEN8/Nxo0bdffddys/P189evTQ6NGjFRwcbDrWeWVmZmr2\n7Nk6ffq0pkyZoscee8x0pEvizXMDM5gZwD6UaQ+xUMET3jg3ZWVlGjdunGbMmKHrr79ezz77rOrV\nq2c6VqXMmzdPqampuuGGG5SRkeG1h35449zALGYGsA9l2kMsVPCEt83Nzp071a1bNx0+fFiPPvqo\n7rjjDtORqiwvL0+///3vdejQIb355ptKSkoyHanKvG1uYB4zA9jHEVdAnDlzpqKiohQQEFDpj549\ne5qODfi0Dz/8UDfccIPq16+vd9991yuLtCRFRERo1qxZ6t+/v5KTk/XLX/7SdCQAgA8x/q6hwYMH\na/HixXK5XHK5XBVftyzrrM9/LDQ0tOISxACq329/+1tNmjRJPXv21KOPPmo6TrV46KGH1LZtW730\n0kvKycnR2rVrFRISYjoWAMDLGT3MIyUlRbNnz9bkyZM1ZMiQs84PGx8fr9mzZysmJsZUvIviJTR4\nwhvmZsCAAUpLS9OvfvUr9ejRw3ScanfgwAE98cQTCgoK0rZt29S4cWPTkX6WN8wNnIWZAexjrEyn\np6crOTlZOTk551xkobCwUOHh4SorKzMRrVJYqOAJp89Nr169lJGRoYkTJ+r66683HafGFBcXa9y4\ncTp58qS2b9+u5s2bm450UU6fGzgPMwPYx9gx0xMmTFB6evp5r1aWnp7OYRyAzfr06aPMzEy98sor\nPl2kJSk4OFhvvvmm6tWrp+uuu06HDx82HQkA4KWMlenMzExdddVV5/3ewoULFR8fX+n7ys3NVXJy\ncjUlA/zPiBEj9Nlnn2ny5Mlq3bq16Ti2CA4O1uuvv67atWvruuuuU3FxselIAAAvZKxMn2+PdLnF\ni9cjpn4AACAASURBVBdXes90amqq4uLidOTIkeqKBviVN954Q/PmzdPvf/97XXvttabj2Kq8UJ86\ndUpdu3Y1HQcA4IUccWq8H0tPT5ekn90zPWHCBA0ZMkRHjhxRVFSUHdEAn/PVV1/pscce0/DhwxUb\nG2s6jhEhISF65ZVXlJOTwytcAIAqc9xFW0aPHq3Zs2fL7XZX+md69uypsLAwLVy4sAaTnY03d8AT\nTpqbwsJCRUREqF27dvrd735nOo5x//znPzVx4kS98847Gj58uOk4Z3HS3MA7MDOAfRy3Z3rRokVq\n1apVlX6GBQOoul69eqlu3br67W9/azqKI3Tu3Fl33323Ro4cqcLCQtNxAABewlFlOjc3V0VFRX77\ncjNgl7feekvr16/X888/r4AARy0DRiUlJal+/frq3bu36SgAAC/hqEfR1NRUSVKHDh0MJwF8V2Fh\nocaMGaN+/frpyiuvNB3HUQICAvTcc89p7dq1evvtt03HAQB4AUeV6fJjnn3xqmuAUwwaNEgNGjTQ\nqFGjTEdxpKuuukq9e/fWmDFjVFpaajoOAMDhHFOmCwsLtWHDBrlcLsdeQhzwdjt37lRmZqbGjx9v\nOoqjjR49WpZl6YknnjAdBQDgcI4p0+WnxKvqmw8BVN7QoUPVsmVLnrD+jICAAA0dOlRTp07VyZMn\nTccBADiYY8q0y+VSWFiYUlJSTEcBfNLatWuVk5Ojxx9/3HQUr3DvvfeqTp06GjlypOkoAAAHc9x5\npj2RkJCg8PBwzjMNxzM5Nx07dlRBQYFee+01I9v3RsuXL9fMmTN16tQpBQYGGsvBeoOqYmYA+zhm\nzzSAmlNYWKj169drxIgRpqN4lV69eqlWrVqaMmWK6SgAAIeiTAN+YMKECapfv77at29vOopXCQgI\nUJcuXfTXv/7VdBQAgENRpgE/8O6776pnz56mY3ilhx9+WAcPHtTatWtNRwEAOJBPlOnDhw/ryJEj\npmMAjpSVlaWTJ0/qvvvuMx3FK4WGhioiIoJDPQAA52XuHTWXaNasWVqxYoVyc3Mrzk8dHx+vVq1a\nKSEhgQtSAP/PX/7yFzVr1kwhISGmo3itzp07KzMz03QMAIAD+cTZPEzgndLwhIm5adKkiTp37qyk\npKQa39akSZN04sQJ5efna8SIEbrlllskScePH9ff//535efnV9x20KBBateuXY1nqg5HjhzRgw8+\nqO+//14tW7a0ffusN6gqZgawj08c5gHg/A4dOqRDhw5pwIABNb6tN954Q71799af/vQndenSRZMn\nT9amTZt04MABPfvss7r77rv1pz/9SRMmTJAkPfvsszpx4kSN56oOYWFhCg0N1auvvmo6CgDAYbz2\nMA8AP2/RokUKCQlRkyZNanQ7Bw4cUH5+fsWe5ubNm0uSUlNTdfz4cb344ouqU6eOJGnu3LnatGmT\nXC6X8vPzveaqp1FRUVq1apXpGAAAh2HPNODD0tPT1bRp0xrfzldffaWuXbtWfH7gwAFJ0qZNm/TQ\nQw9VFGlJioyMlCR16dLFa4q0JN10003Kzc01HQMA4DCO3zN9+PBhrVmzRoWFhapTp47at2+v1q1b\nm44FeIUNGzbY8v/lyy+/1Isvvljx+a5duyRJMTEx5xwXfe+99+ree++t8UzVrUuXLpo7d66Ki4sV\nHBxsOg4AwCEcW6Y3bdqkSZMm6aOPPlJwcLDcbrcCAgJUWlqqdu3aacKECUpMTJTL5TIdFXCs/fv3\na/DgwTW+nZiYmLP2Pu/cuVOSdNddd3l0f88++6zq1q1bcXy1EzRv3ly1atXSp59+qsTERNNxAAAO\n4cjDPObOnasuXbpo0aJFOn36tI4ePaoTJ07o2LFjOnXqlP71r3/p/vvv14gRI+R2u03HBRyprKxM\nxcXFtpwx48eXKT9w4IBOnjwpl8ulmJiYKt/X8ePH9d133+naa6+t1O3Hjx+vN998s8rb8US9evWU\nlZVly7YAAN7BcXumU1NT9eijj+rkyZMXvd3/1969B0dV3n8c/+wCSYoZSIJW6qU/clFbrQpJbB2l\nlxEI2qljJSTSzjhDNSZSnGHKUBLttJ3RtoQExT/UmWy8MNY2BbKpHR0LyYYWHNTaZEMcL0XIBmoR\nnSLdKLdCzP7+yOw2Idd9snvO2d33aybjZrNnz/dkHg+fPPs9zzl58qS8Xq/S09PV2NhoUXVA4jh8\n+LAkafbs2Zbud+/evZKkiy++eNhs9WRlZmaqqalpUq8NX/h42223Rb0fE5mZmTp06JAl+wIAJAZH\nzUyfOXNG99xzz4RBOuzUqVP6/e9/rzfeeCPOlQGJ56233tKMGTMs3+++ffskyWhWOlpz585VU1OT\nli5dGvd9SYN/mHzwwQeW7AsAkBgcFaa3b98e9SLzZ86cYe1XYBT79+9Xenq65fvt7u6WZE2YtlpO\nTk5kpRIAACSHhelNmzbpxIkTUW0zMDCgl19+WX19fXGqCnCuPXv2aNmyZXrssce0Z88effbZZ5Gf\nffzxx5bfQvzgwYOSZNwv7XSzZ8+O+hwFAEhujuqZ7u3tNdrO5XLpvvvuU05OTowrGt/9999v6f6Q\nHGI5bl588UX9+9//1iuvvKL09HSdPn1aX/ziF1VcXKxjx45ZvtpNeFbapF96586dOnr0qHp6erRy\n5Url5+eP+roTJ06oublZ0uASfKtXr47cJCbepk+fzkXPAIBhXKFo+yriaObMmTp9+nTU27lcrqjb\nQ4BkMXT8p6WlaWBgQP39/crMzFRGRoaee+45y2r5+c9/ru7ubpWWlg5b4WMie/fuVWZmpq6//no9\n8cQT6unp0ebNm0e8LhykV65cKUmqra3VyZMn9cgjj8TqEMb1zDPP6I033tDRo0ct2V8Y5zhEizED\nWMdRbR6mqw5kZGTo0KFDCoVCln1JsnR/fCXHV6zHzVNPPSVJuvTSS3XHHXdow4YN2rVrlz799FNV\nVFTYNjMdbYvHgQMHIkv4hVtFRtPc3Kzy8vLI9ydPntTHH39sUKmZzz//XNOmTbNsfwAA53NUm8cP\nf/hDPfHEEzp79mxU211++eX68pe/HKeqAOeqrKzU3XffrczMzBE/mzlzps6dO2dZLeEL81wuV1Rr\nW3/00UdasGBB5HEgEIjMPA918uRJzZo1a8TNYcLbWuH06dPc/RAAMIyjZqYfeOABud3RlZSZmama\nmhruhIiUNG3atFGDtCRdeeWVRm1TU3HBBRdEfavwuXPnRsL3jh07JI1+58Tz3/vgwYM6deqUpRc6\nHjt2TBdeeKFl+wMAOJ+jwnRubq5uv/12feELX5jU691ut2bNmqW77rorzpUBiefqq6+2dGY6vOZz\nNL3S59uxY8eIW5OPxY4l+Pr6+nTZZZdZtj8AgPM5KkxL0vPPP6/rrrtuwkA9bdo0ZWdna/fu3UZ3\nWQOS3TXXXBO5pXgi2Ldvn06dOqXS0tJJvz4zM1MXX3xxnCv7n88++0z/93//Z9n+AADO57gwnZGR\nod27d+vuu+9WRkbGiKA8Y8YMZWRk6Oabb9a+fftUUFBgU6WAs82cOVNut1s9PT12lzIpO3bsiKzo\nIUktLS3jvr67uzuq3uxYOHHihL7yla9Yuk8AgLM5LkxLUnp6uhoaGnT06FH95je/0U033aSrr75a\nxcXF+vGPf6y3335bu3fv5uNWYALZ2dl688037S5jUl577TUtXLgw8v2nn3465mvDK35Y2eJx5swZ\nnTlzRrfffrtl+wQAOJ+jVvM4X1ZWltasWaM1a9bYXQqQkK666iq99957dpcxaeFw3NLSottuu23M\n19nRL93R0aEZM2bokksusWyfAADnc+TMNIDYWLhwoY4cOWJ3GZOyevVq/fnPf9aWLVuUn58/bi+0\nHf3SHR0d+tKXvmTZ/gAAicHRM9MApmb58uWqr69Xf3+/pk939v/uS5cu1dKlSyf12u7ubt18881x\nrmi4999/X4WFhZbuEwDgfMxMA0nshhtu0IwZM7Rr1y67S4mZffv2SZK++c1vWrbPgYEBHTlyRJWV\nlZbtEwCQGAjTQJJbsGCBfD6f3WUY27Jli3bu3Bn5fseOHZo7d65uuukmy2p49dVX5Xa7R72ZDAAg\ntRGmgSR3zz33JMzyeOc7ePCgWlpaIvXv3btXb731lh555BFL69i5c6euv/76qO/QCgBIfvzLACS5\nlStX6ty5c+ro6LC7lKgVFBTopptu0gUXXKAnn3xS3d3devrppy298HBgYED79++f0p0dAQDJyxUK\nhUJ2F5GIXC6X+NUhWnaNmxtvvFHHjh3T448/bvm+E90rr7yip59+WqdOnbLtIk7ON4gWYwawDjPT\nQArYtGmTAoGA+vr67C4l4fzxj3/Ud7/7XcevhgIAsAdhGkgBCxcu1IUXXqjnnnvO7lISyuHDh/XR\nRx9p8+bNdpcCAHAowjSQItasWaM9e/bo7NmzdpeSMJ566ildeeWVys3NtbsUAIBDEaaBFPHggw8q\nPT1dDQ0NdpeSEA4dOqT33ntPzz77rN2lAAAcjDANpAi32626ujq1t7fr1KlTdpfjeI8++qjmz59v\n+Z0WAQCJhTANpJBVq1Zpzpw5euyxx+wuxdG6u7t1+PBhNTU12V0KAMDhCNNAinnmmWf05ptv6t13\n37W7FEfq7+/Xxo0bdeutt+qqq66yuxwAgMOxzrQh1vCECaeMm6VLl+r111/XCy+8wF39zlNfX6/O\nzk4dP35caWlpdpcjyTnjBomDMQNYh39FgRT0pz/9SZ9//jlLvp3nnXfe0auvvqo//OEPjgnSAABn\nI0wDKSgjI0NNTU3avXu33njjDbvLcYRTp07p4Ycf1q233qrvfe97dpcDAEgQtHkY4iM0mHDauKms\nrNSzzz6rJ598Updeeqnd5dhmYGBAq1evVn9/vz744APH3e3QaeMGzseYAaxDmDbEiQomnDhuiouL\ntX//fj333HPKyMiwuxxb1NbWqqOjQ4cOHdLcuXPtLmcEJ44bOBtjBrAObR5AinvttdeUnp6utWvX\namBgwO5yLNfU1KTXXntNra2tjgzSAABnI0wDKS4tLU379u3TJ598knKBuqWlRU1NTfJ4PPrWt75l\ndzkAgAREmAagyy67TG+//baOHDmSMoF627Zt2rJlizZv3qyKigq7ywEAJCh6pg3RjwYTTh83vb29\n+trXvqacnBw9+uijmjlzpt0lxUVjY6NeeuklPfnkk1q1apXd5UzI6eMGzsOYAaxDmDbEiQomEmHc\nfPjhhyosLNSJEydUX1+vyy+/3O6SYqa/v18/+9nPtH//fv32t7/VD37wA7tLmpREGDdwFsYMYB3C\ntCFOVDCRKOOmv79f3/72t/W3v/1Na9euTYp+4nBP+Llz5/T666/rmmuusbukSUuUcQPnYMwA1qFn\nGsAI06dP1969e7V69Wpt2rRJjzzyiM6ePWt3WcZeeukl3Xvvvbrooov04YcfJlSQBgA4GzPThvir\nHyYScdzs2rVLd955p86dO6f169eruLjY7pImLRgM6pe//KUOHTqkmpoa/frXv7a7JCOJOG5gL8YM\nYB3CtCFOVDCRqOOmv79fK1asUEtLi66++mqtW7dOF154od1ljWlgYEDPPvusXn75ZV122WVqb29X\nfn6+3WUZS9RxA/swZgDrEKYNcaKCiUQfN3v37tXdd9+tQ4cO6Rvf+IZ+8pOfOG7FD6/Xq6amJrnd\nbv3iF79QTU2N3SVNWaKPG1iPMQNYhzBtiBMVTCTLuPF6vVq1apU++eQTLViwQPfee6+tq36cOXNG\nv/vd79TW1qb//ve/euCBB1RfX6/p06fbVlMsJcu4gXUYM4B1CNOGOFHBRLKNm+eff14PP/ywenp6\nNHfuXN1xxx0qKSlRWlqaJft/55139MILL+jdd99VZmamVq5cqQ0bNjhutnyqkm3cIP4YM4B1CNOG\nOFHBRLKOmwMHDmjt2rVqa2vT2bNndckll+jGG2/U97//fWVlZcVsP/39/frrX/+qtrY2HTx4UOfO\nndNXv/pV/epXv9Kdd94Zs/04TbKOG8QPYwawDmHaECcqmEiFcbNnzx49/vjj+stf/qJgMKi0tDRd\ndNFFys3N1bXXXqt58+Zp3rx5484e9/f361//+pf++c9/6p133tH777+vo0eP6uTJk0pLS9P8+fN1\nzz336Ec/+pFls+B2SoVxg9hizADWIUwb4kQFE6k2bj799FO9+OKL2rFjhzo7O3XkyBGdPn1aAwMD\ncrlcmj59utxut1wul6TBVTj6+/sjP09LS1NOTo6uu+46fec731FpaamuuOIKm4/Keqk2bjB1jBnA\nOoRpQ5yoYIJxM6i/v1/79+/XP/7xD508eVJnz57VjBkzlJaWpnnz5unaa69VZmam3WU6BuMG0WLM\nANYhTBviRAUTjBuYYNwgWowZwDrcThwAAAAwRJgGAAAADBGmAQAAAEOEaQAAAMAQYRoAAAAwRJgG\nAAAADBGmAQAAAEOEaQAAAMAQYRoAAAAwRJgGAAAADBGmAQAAAEOEaQAAAMAQYRoAAAAwRJgGAAAA\nDBGmAQAAAEOEaQAAAMAQYRoAAAAwRJgGAAAADBGmAQAAAEOEaQAAAMAQYRoAAAAwRJgGAAAADBGm\nAQAAAEOEaQAAAMAQYRoAAAAwRJgGAAAADBGmAQAAAEOEaQAAAMAQYRoAAAAwRJgGAAAADBGmAQAA\nAEOEaQAAAMAQYRoAMIzf75fb7R71q6+vz/h9g8Gg8vPzR33flpaWGB4BAFiHMA0kGJ/Pp+zsbJWX\nl9tdStxVVVWpvr7e7jJSTmFhoQYGBtTT06P169crLy9PkuRyudTR0WH8vtXV1XK5XJH3am5uVjAY\n1MDAgJYtWxaT2gHAaoRpIMGUlZWpr69P7e3tdpcSF4FAQB6PR/n5+WpsbNTx48ftLill5ebmas6c\nOaqqqpIkhUIhBQIBo/fy+/0qKCiIbF9ZWally5Zp1qxZMasXAOxAmAYSTFVVVSRoJpOioiLl5OSo\nvLxcbrdbZWVldpcESX//+9+1fPnyyPc9PT1G71NbW6vFixdHvl+yZMmUawMAJ5hudwEAolNbW6va\n2lq7y4i5zs7OYd/X1dXZVAmG8vv9ys3NVW5urnp7e9Xb2xv1e9TV1emhhx5Sa2urpMEWj6HBGgAS\nGTPTAIBRBQIBFRYWSlLkv9G2eQQCAfX29mr+/Plqa2uTJOXl5dHeASBpMDMNABiVz+fT17/+dUmK\nXIQYbZiuqanR008/LUmRixeZlQaQTJiZBgCMyufzRYJvQUGBpMHl7SarublZJSUlmjVrlgKBQGRZ\nPfqlASQTwjQAYFR+v1/z58+X9L+ZaUmT7pvetm2bKioqJA0Gc2mwXzrcMgIAyYAwDThceHavpKRE\nxcXFSbskHpwlEAgMC9BDH0+m1aO6unrYRaThfumsrCzNmzcvdoUCgM0I04CDVVVVqbOzU62trWpt\nbVVxcbGWLFmirq4uu0tDkvP5fCopKYl8n5ubK2lya02H76A4NDSHZ6bplwaQbAjTgEPV1dXpiiuu\n0IYNGyLPhWcHt27daldZSBFD+6XDZs+eLWnitaZra2uHjVv6pQEkM8I04ECBQEA+n0/r1q0b9nw4\nxJis9QtEY2i/dFj4j7nxxp/H49H9998/7Lmh/dLMTANINoRpwIGqq6tVU1Mz4vnw0mJD+1ejFf4I\nPtZf06ZNM64JznJ+v3RYfn5+5OejCQaD8vv9uuWWW4Y9T780gGTGOtOAAxUUFIwIJJLU1dUll8s1\npY/KCwsL5ff7p1LeqLKysmL+nrCH3+8f1i8dNtFa0zU1NaPeuZJ+aQDJjDANONDQftOwcCCRNGrQ\njsb5H98DQ/l8vhGtGtL/ZqZHW2va5/OpoKBgxJ0N6ZcGkOxo8wASRPijctboRbz5fL5R/+Aab61p\nj8czosc//F4S/dIAkhdhGkgQfFQOKwSDwTFbdoqLiyUNLo83tFWourpaDz300Kjb0C8NINkRpoEE\nEYt+aWAiPp9vzDEWXhpP+t/MdLiNY6zWIf4IBJDsCNNAAohlv7TP52M1D4ypra1t3D/YwjdvCS/T\nONZFh9LgLDf90gCSHRcgAgkglv3SixcvViAQGPUisqlgNY/k0N7eroaGhjF/npeXp97eXvX09Mjr\n9WrFihUjLjoMo18aQCogTAMJINYflSdr72pdXZ2Kioq0aNEio+09Ho/mzJmj0tLShNx+qsc/Xr90\nWF5entrb2xUIBLRt27Zx78ZJvzSAVECYBhIA/dITq6qqUmNjoySpoaFB9913X1Tbl5WVyev1SpK2\nb98edaC1e/upHr80uCRjTk7OuK8pKCiQNNgzPbT9aDThn4cvXASAZETPNOBwseyXThTBYDAy49nW\n1hbpux1PZ2dn5PF4bQpj6erqijw2uV273dubHn/4roXV1dWqr6+Xz+eT1+sd88Ys4Z7p9evXjzrb\nHAwGFQgEVFdXFzmOUCikrq6umLcWAYAjhGCEXx1MmIyb9evXh1wuV6i4uDgOFTnH8uXLQ9nZ2SGX\nyxVyuVwht9sd+Qo/l52dHSopKRl1e5/PF8rOzg4VFBSECgoKot7/0O37+voctf1kxo3J8ff09ER+\nt0N/5+HHXq93xDaBQGDMsVhZWTnm+4W/6uvrJ1UbpoZ/owDruEKhUMjeOJ+YXC6X+NUhWibjpqio\nSF1dXaqrqxv1phgYqaSkRK2trXaXETPRjptkO35Ej3+jAOvQ5gE4WDAYjPRLsxoCAADOQ5gGHGzb\ntm2SBldDGOumGBhuMitSJLNUP34AsBphGnCAsrIyFRQUDLsITZI2btw47L+YmMfj0YoVK+wuwzap\nfvwAYDXCNGAzj8cjr9er3t7eYSsyhFdDqKqqUkVFhY0VJhafz6dly5bZXYZtUv34AcBqXIBoiIs7\nYCInJ0f/+c9/7C4DQJLLzs7W8ePH7S4DSAmEaUOEacRSfX29tm7dqpycHB0/flz5+fnauHEjd42L\nQldXl9rb25NyxZPJnG+S+fgBwMkI04YI0wCswvkGAJyLnmkAAADAEGEaAAAAMESYBpJEc3Oz3G73\niK+CgoJJbZ+fnz/q9m63W3fddVecqwcAIDERpi1SV1en7Oxs1dfX211K3FVVVaXEcTrN8uXLFQwG\nFQgEtH79+sjzgUBAXq93wu27urrU2dmp5cuXSxoM183NzQoEAtq6dWvc6gYAIJERpi1SU1Ojvr6+\nyB3tkk0gEJDH41F+fr4aGxtZkskms2bN0rx58xQKhYbdfryhoWFS2y5YsEC1tbWSpM7OTi1btowV\nRQAAGAdh2iKLFy+OBM1kUlRUpJycHJWXl8vtdqusrMzukiCpvb1dDQ0NysvLkzR4I4/e3t5JbRsI\nBFRWVqZZs2bFs0QAAJICYdoira2tOnDggObPn293KTHV2dmp48ePq6OjQxUVFcrJybG7JGgwEOfm\n5qqqqiry3GRmpyWpra2NHmkAACaJMA0kGZ/PpyVLlkiSKisrI897PJ5Jbd/e3j6sRQQAAIyNMA0k\nmba2tkiYnj17duSCwmAwqPb29gm3DwaDtHgAADBJhGkgyZw/szy01WPjxo3jbuvz+ZiVBgAgCoRp\nIMkEAoFhK3AsWrRIWVlZkgbDcl9f35jb+v1+FRcXx7tEAACSBmE6Tvx+v0pKSlRcXKzi4uJJfbwO\nTNXQfumhHnzwwcjj8XqnmZkGACA6hOk4CAQCqqmpUXNzszo6OpSbm6slS5aMOyMIxMLQfumhhl6I\nON6qHufPagMAgPERpuOgvLxczc3NkYu4gsGgpMFZPyCeurq6Rp1Znj17duT5QCAw6iclfr9fhYWF\nca8RAIBkQpiOsebmZq1YsWLYaggdHR2SpOzsbLvKQoro6OgYc2a5uro68ni02Wmfz6eSkpJ4lQYA\nQFIiTMdYbW3tsI/U/X6/+vr65HK5dMsttxi/r9/vl9vtjvnXtGnTYnHYcAC/3z9uv/PQCxGbm5tH\ntB3RLw0AQPSm211AslmyZMmwWenwDGB4rV9ThYWF8vv9U3qP0YTDFRLfZGaWKysrVVdXJ2nwQsSf\n/vSnkZ/RLw0AQPQI0zG2YcOGYd83NjbK5XINW+vXVLLdihyx5fP5JrzLYVVVVSRMNzQ0RMI0/dIA\nAJihzSOOwsEmKytrSi0ewGRMZmY5Nzd32IWIXV1dkuiXBgDAFGE6jsItHkN7qIF4iGZmeeinJOFP\nUuiXBgDADGE6TsKzfrFq8QDGE83McmlpaaRX3uv1qq+vj35pAAAMEabjJDwrXVhYOCykTNTTOhaf\nz8dqHhhTtDPL4U9LQqGQqqurVVRUFK/SAABIalyAGCfNzc2SNGJW2uPxGLV9LF68WIFAIHIDmFhh\nNY/kEO3M8tALET0eT+QxAACIDmE6DgKBgHp7e+VyuVReXh55PnxDF1PJ+jF8XV2dioqKtGjRIqPt\nPR6P5syZo9LSUsu3n+q+p3rs0mC/dF5eXlTbhC9E9Pl8crlc9EsDAGCIMB0H4dnjvLy8YWtO19bW\nateuXXaV5UhVVVVqbGyUNNgac99990W1fVlZmbxeryRp+/btUYfaqWw/1X1P9djDNmzYEHWYDu8/\nfIt7ll0EAMAMPdNxUFhYOKJ9ory8XHV1dcPCdbIJBoPaunWrJKmtrW3EHfZG09nZGXk82i2uJxJe\n2k2Sent7Ld1+qvueyrEHg0H5/X5VVVXJ6/Vq27Zt8nq9UbUBhS9EZFYaAABzrlAoFLK7iETkcrk0\n3q+ut7dXZWVlkqScnBzV1NQk5VrTZWVlam9vj4Q4l8sV+Vn495OVlaUbbrhBO3fuHLF9e3u7ysrK\nNGfOHEnSgQMHotr/0O07Ozuj/mNlKtvHct/S5I/d7/eruLg48n34dx4KheRyudTZ2TnpmeaamhoV\nFBSooqIiqtphrYnONwAA+xCmDfGPW+yVlJSotbXV7jJskcrHjolxvgEA56LNAwAAADBEmIYjDaan\nggAAAnJJREFUBIPBlF2mL5WPHQCAREeYhiN4PJ4pLRuYyFL52AEASHT0TBuihzG2UrlnOJWPHZPD\n+QYAnIt1pg1lZ2cPW7kCU5fKv89UPnZMLDs72+4SAABjYGYaturq6lJ7e7vWrVtndymWS+VjBwAg\nWRCmAQAAAENcgAgAAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHC\nNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAA\nAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCI\nMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0A\nAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAY\nIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwD\nAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAA\nhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjT\nAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAAgCHCNAAAAGCIMA0AAAAYIkwDAAAAhgjTAAAA\ngKH/B0JHamWXvI1KAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1128db250>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finding hyperparmeters is the unprincipled black magic of machine learning. I basically tried some numbers until they worked."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "D    = 10\n",
      "tau0 = 5\n",
      "tau1 = 5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u = pymc.Normal('u', 0, tau0, size=(D, N))\n",
      "v = pymc.Normal('v', 0, tau0, size=(D, M))\n",
      "\n",
      "z = pymc.Lambda('z', lambda u=u, v=v: np.dot(u.T, v))\n",
      "r = np.empty(nTrain, dtype=object)\n",
      "\n",
      "for n, (i, j) in enumerate(izip(users, items)):\n",
      "    r[n] = pymc.Lognormal('x_%d_%d' % (i, j),\n",
      "                           z[i-1, j-1],    # mean\n",
      "                           tau1,           # precision (confidence)\n",
      "                           observed=True,\n",
      "                           value=ratings[n])\n",
      "\n",
      "model = pymc.Model([u, v, z, r])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Posterior Inference and Markov Chain Monte Carlo"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# There was a bug in PyMC; this is just to hack around code.\n",
      "if 3 in model.__dict__:\n",
      "    del model.__dict__[3]\n",
      "    \n",
      "pymc.MAP(model).fit()\n",
      "mcmc = pymc.MCMC(model)\n",
      "mcmc.sample(10000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r",
        "[****************100%******************]  10000 of 10000 complete"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Joint probability:\n",
      "\n",
      "\n",
      "$$p\\left(\\left\\{ r_{ij},z_{ij}\\right\\} _{ij\\in I}\\left\\{ \\mathbf{u}_{i}\\right\\} _{i=1}^{N},\\left\\{ \\mathbf{v}\\right\\} _{j=1}^{M},\\tau_{0},\\tau_{1}\\right)=\\prod_{ij\\in I}\\mathcal{LN}(r_{ij}|\\mathbf{u}_{i}^{\\top}\\mathbf{v}_{i},\\tau_{1})\\prod_{i=1}^{N}\\mathcal{N}(\\mathbf{u}_{i},\\tau_{0}I_{d})\\prod_{j=1}^{M}\\mathcal{N}(\\mathbf{v}_{j},\\tau_{0}I_{d})$$\n",
      "\n",
      "Posterior probabilities:\n",
      "\n",
      "$$p\\left(\\left\\{ z_{ij}\\right\\} _{ij\\in I}\\left\\{ \\mathbf{u}_{i}\\right\\} _{i=1}^{N},\\left\\{ \\mathbf{v}\\right\\} _{j=1}^{M},\\tau_{0},\\tau_{1}\\vert\\left\\{ r_{ij}\\right\\} _{ij\\in I}\\right)=\\frac{\\prod_{ij\\in I}\\mathcal{LN}(r_{ij}|\\mathbf{u}_{i}^{\\top}\\mathbf{v}_{i},\\tau_{1})\\prod_{i=1}^{N}\\mathcal{N}(\\mathbf{u}_{i},\\tau_{0}I_{d})\\prod_{j=1}^{M}\\mathcal{N}(\\mathbf{v}_{j},\\tau_{0}I_{d})}{\\int\\prod_{ij\\in I}\\mathcal{LN}(r_{ij}|\\mathbf{u}_{i}^{\\top}\\mathbf{v}_{i},\\tau_{1})\\prod_{i=1}^{N}\\mathcal{N}(\\mathbf{u}_{i},\\tau_{0}I_{d})\\prod_{j=1}^{M}\\mathcal{N}(\\mathbf{v}_{j},\\tau_{0}I_{d})d\\mathbf{z}d\\mathbf{u}d\\mathbf{v}}$$\n",
      "\n",
      "Posterior probabilities are the **cornerstone of Bayesian statistics**. They give us information like\n",
      "\n",
      " - Predicted ratings, $r_{ij}$\n",
      " - Error bars on predictions, $\\sqrt{\\mathbf{V}[r_{ij}]}$\n",
      " - Shape of latent representations of movies and users, $p(\\mathbf{u}_i | \\ldots)$, $p(\\mathbf{v}_i | \\ldots)$.\n",
      "\n",
      "\n",
      "**Markov chain Monte Carlo:** (MCMC) Do not evaluate the posterior explicitly. Construct a _Markov chain_ whose _invariant distribution_ is the posterior.\n",
      "\n",
      "In fact, with MCMC we never know how to write down the posterior. _Eventually_, our Markov chain generates samples from it.\n",
      "\n",
      " - Markov chain: a process whose future depends only on the present, not on its own history\n",
      "     - Computationally easy to simulate.\n",
      "     - Analytically, rich theory, which justifies the method.\n",
      "     \n",
      "Gibbs sampling: The easiest MCMC method.\n",
      " \n",
      "  - Let $x_1, \\ldots, x_n$ be random variables, and $\\mathbf{y}$ be observed.\n",
      "  - For $t = 1, 2, \\ldots$:\n",
      "      - For $i = 1, \\ldots, n$:\n",
      "      - Sample $p(x_i | x_1, \\ldots, x_{i-1}, x_{i+1}, x_i, \\mathbf{y})$\n",
      "\n",
      "Properties:\n",
      "\n",
      "  - Each inner loop generates a vector $\\mathbf{x}^{(t)}_{1:n}$. These vectors form a Markov chain.\n",
      "  - At the invariant distribution (when you run then chain for \"long enough\"), the Markov chain generates samples from the conditional $p(\\mathbf{x} | \\mathbf{y})$.\n",
      "  - For any subset of variables $x_i, x_j, x_k$, keeping just those variables from the Markov chain is a sample from the marginal distribution $p(x_i, x_j, x_k | \\mathbf{y})$.\n",
      "  \n",
      " "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "\n",
      "Posterior Analysis"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "zsamples = mcmc.trace('z')[5000:]\n",
      "zmean    = np.mean(zsamples, 0)\n",
      "\n",
      "meanExpZ = np.mean(np.exp(zsamples), 0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We learned a _distribution_ over _all_ log-ratings $z_{ij}$. Let's take a peek at the first rating. The shape of the histogram captures our knowledge, including knowledge of uncertainty."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(np.exp(zsamples[:,1,1]))\n",
      "\n",
      "iFirst, jFirst, rFirst = users[1], items[1], ratings[1]\n",
      "print meanExpZ[iFirst,jFirst]\n",
      "print rFirst"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.36496704562\n",
        "1.5098059294\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEjCAYAAABAaxQzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3T1v48iaL/A/jQl272JaLw3cuC0q2awl2tkGxlgUsGlb\n9JwT3KzH8uTbLWOT7Y56JKcbrOjzBSzZ2SbTogADF7jJkanNt0X1F2iLmtkT7GJm6gY+1MjWGyWK\nJZH+/wABM2Q9rDKb5EMWX0oRQggQERHFzM6mG0BERBQGJjgiIoolJjgiIoolJjgiIoolJjgiIool\nJjgiIoolJjgiIoolJjgiIoqlr5YNKJfLyGazePPmja/yruviw4cPaLfbcBwH6XQa+Xwe+/v7ODk5\nQSKR8LUM0zThOM6DduRyOd9tCBJPREQRJHzo9XqiXq+LTCYjFEURZ2dnfsJEq9USmqaJ6+trMRwO\nhRBC2LYtKpWKUBRFKIoiTNOcu4xmsylUVRXtdns0zXVdYRiGMAxjYRuCxhMRUTTNTXD5fF6kUimh\naZq4uLgYJSY/Ca7X6wlN02bOt217YZJrtVpCURTR7Xanztd1Xei6PrOOoPFERBRdvq7gPNVq1XeC\nMwxjdNU2S61WGyU5x3EezBsMBiKZTIrT09OZ8V6SnJYgg8YTEVG0hfKQiXev69mzZ3PLvXnzZnQP\nrl6vP5jXaDQwHA5RLpdnxudyOWQyGVSr1Yl5QeOJiCjaQklwtm3j6uoKxWIRw+Fwbtn9/X0AwNXV\n1YPp1WoViqLg5cuXc+NzuRwcx0G3211rPBERRVuoV3CWZaHdbs8tu7u7CwAYDAajaa7rot/vI5lM\nLqwrk8kAAC4vL9cWT0RE0RdKgtN1HclkEqqqolAozC17d3cHAEin06NpnU4HwO/JZ55sNgsAD67A\ngsYTEVH0Lf0enB+5XG6UuBbxrvby+fzEtPGkN4t3BTj+jlvQeCIiir6NfsnEdd3RlZOu6w+mL2s8\noQaNJyKi6Ntogms0GgAAVVXx+vXr0fQvX74AgK97aJ7xpBY0noiIom+jCc570vHxY/pBkw2TFRER\nbSzBVSoV9Pt9vH37Fq9evXowb5nuwmn32YLGExFR9G0kwdm2jfPzcxiGgQ8fPkzMXybpTEtmQeOJ\niCj6QnmKch7XdXF4eAjDMGa+e+bdO1u1qzFo/CLZbBa9Xi+UZRMRxZWqqvj06ZO0+qRfwR0eHqJY\nLK79xeplHigJGt/r9SDuv+MZyd+//Mu/bLwNT7X9UW4727/5X9TbL/vCQGqC03Ud2Wx2YXLzXr72\n033oXaWNd0sGjSciouiTluDK5TLS6bSvKzfvCyR+Xr6e9qJ40HgiIoo+KQmuVqvh8+fPc5Pb+fn5\n6L/39vYA+LuH5l3yeh9tXkd83B0cHGy6CYFEuf1RbjvA9m9a1NsvW+gJ7urqCu12Gz/++OPccq1W\na/TfiUQChUIBQgj0+/25cd4VWKlUWlt83EV9J4ly+6PcdoDt37Sot1+2UBOcbdswTXNhcnNdd+LD\nyN44bpZlzY1rt9vQNA0vXrxYazwREUVbaAnOcRycnZ2hXq/Ddd2ZP9u2UalURg+GeI6OjpDP5ycG\nQh3XaDSgKAqazebEvKDxREQUccKnwWAg8vm8UBRFaJomXNedWzaTyQhFUXz/rq+vJ5bjOI5IpVKi\nVqtNzOv1ekJRFHF+fj6zHUHjZ1litRER0V/JPnbOra1UKolUKjVKQjs7O6OfNy2VSolisfgg7u3b\ntxPlF/36/f7UNjiOIzRNE+VyWdi2LRzHEfV6XaiqKi4uLhb+gUHjp2GCIyJanuxjp/LXSrdet9tF\np9OB67pQVXXi+5Vhx49TFAURWW1ERFtD9rEzMglumzDBEREtT/axc6PD5RAREYWFCY6IiGKJCY6I\niGKJCY5oRc+epaEoirTfs2f8IDjRMviQyQr4kAkB99sBIHM74HZH0caHTIiIiNaACY6IiGKJCY6I\niGKJCY6IiGKJCY6IiGKJCY6IiGLpq2UDyuUystks3rx54zvGdV2YpjkaPdtbTi6Xi0Q8ERFFkJ8h\nB3q9nqjX66Mx3s7OznwPV9BsNoWqqqLdbo+mua4rDMMQhmFsffw0PlcbxRwAAQiJP253FG2yt+G5\nXZSapiGdTuP4+Bg7OzswDGOp5GlZFo6Pj3F1dYVvvvlmND2RSKDRaMB1XRSLxa2NJyKiCFsmG1ar\nVd9XcIPBQCSTSXF6ejqzjG3bQlEUYZrm1sXPs+Rqo5gCr+CIliJ7Gw7tIZNGo4HhcIhyuTyzTC6X\nQyaTQbVa3bp4IiKKttASXLVahaIoePny5dxyuVwOjuOg2+1uVTwREUVbKAnOdV30+30kk8mFZTOZ\nDADg8vJya+KJiCj6QklwnU4HwO/JY55sNgsAD66gNh1PRETRF0qC8943S6cXj1+1u7v7IGYb4omI\nKPpC66Jc1t3d3dbEExFR9IWS4L58+QIAvu6BecaT0qbjiYgo+rbmCm6b4omIKPpCSXDLdPdNu0+2\n6XgiIoq+UBLcMkljWjLadDwREUXf0qMJ+OHd+1q1q3DT8X68e/du9N8HBwc4ODgIrS4ioii6ubnB\nzc3NxuoPJcGtYpkHQrYhfjzBERHRpMcn/+/fv5dafyhdlN7L0366/7yrrPFuxU3HExFR9IWS4Lwv\niPh5edork8/ntyaeiIiiL5QEt7e3B8DfPbBerwcA2N/f35p4IiKKvlASXCKRQKFQgBAC/X5/blnv\nCqpUKm1NPBERRV9ow+V447BZljWzjOu6aLfb0DQNL1682Kp4IiKKttAS3NHREfL5POr1+swyjUYD\niqKg2WxuXTwREUWc36G/B4OByOfzQlEUoWmacF13YYzjOCKVSolarTYxr9frCUVRxPn5+dbGz7LE\naqMYAyAAIfHH7Y6iTfY2PLe2UqkkUqmUUBRFKIoidnZ2Rj9vWiqVEsViceYyHMcRmqaJcrksbNsW\njuOIer0uVFUVFxcXCxu46fhpeKAhIZjgiJYlextW/lpp6LrdLjqdDlzXhaqqePXqVaTixymKAkmr\njbaYoigAZG4H3O4o2mQfO6UluDhhgiOACY5oWbKPnaE9ZEJERLRJTHBERBRLTHBERBRLTHBERBRL\nTHBERBRLTHBERBRLTHBERBRLTHBERBRLTHBERBRLTHBERBRLTHBERBRLX8mu0LZtmKY5+vBxJpOB\nrus4OTlBIpFYGO+6LkzTHI3EDdwPbprL5XzVHzSeiIgiQubQBaVSSRwfH4tutzua5jiOqNVqIpVK\nCdM058Y3m02hqqpot9ujaa7rCsMwhGEYC+sPGu+RvNpoS4HD5RAtRfY2LG00AcMw8Mc//nHmMDX9\nfh+qqqLZbOLo6GhivmVZKBaLsG0bL1++nJhfLBYBAB8/fpy6/KDx4ziaAAEcTYBoWdKPnTKyaLPZ\nFMfHxwvLmaYpUqnUxPTBYCCSyaQ4PT2dGWvbtlAUZepVYND4xyStNtpy4BUc0VJkb8NSHjKp1+vQ\ndX1hOcMw4Lou/uM//uPB9EajgeFwiHK5PDM2l8shk8mgWq1OzAsaT0RE0SMlwfX7fQwGg4Xlksnk\n1OnVahWKokztWhyXy+XgOA663e5a44mIKHqkJLhEIoF6vb6wnPdk43gicl0X/X5/ZvIbl8lkAACX\nl5driyciomiSkuD29/fhOM7oQY5Z6vU6DMN4MK3T6QD4PfnMk81mAeDBFVjQeCIiiiYpCa5SqQC4\nf5IxnU6j3W5PlLEsCxcXF7i4uHgw3buqS6fTC+vZ3d19ELOOeCIiiiYpCW53dxfNZhPAfZehrus4\nPT0dzbcsC2dnZ7BtG8+ePXsQ67ru0vXd3d2tLZ6IiKJJ2qe6jo6OcHt7O+oqNE0T2WwWp6en6Ha7\n6HQ6ePHixUTcly9fAMx+AGWa8aQWNJ6IiKJJ6rcoc7kcPn36hHw+D+C+K9A0zVESmiZosmGyIiJ6\nmqQmONd1YRgGdF1Hs9kcXVXVajXs7e2h3+9PxCzTXTjtPlvQeIqWZ8/SUBRFyo+Itpu0jy1bloXT\n01NYljXqiiwUCjAMA5ZlwbZtqKqK29vbBx8+XibpTEtmQeMpuGfP0vj558XvQa6PrE8BMckRbTMp\nV3C2beP4+PhBcgPu34/7+PHj6AEUANA0DcPhcPT/3lXeql2NQeMpuPvkJiT9iIjuhX4F53VL1mq1\nqQ+RAPcPoPR6Pei6Dsdx8N1336HRaKxc5zIPlKwa/+7du9F/Hxwc4ODgIFCdRERxc3Nzg5ubm43V\nH3qCM00Truvi9evXc8vt7u7i9vYWu7u7uL6+Hk33Xr72033oXaWNd0sGjZ9lPMEREdGkxyf/79+/\nl1p/6F2UrVYLe3t7vsomEglUq1UIIfD582cAv3+BxM/L114Z7ynNdcQTEVE0hZ7gBoOBr89keQ4P\nDx/8v5cc/dxD6/V6AO4/DbaueCIiiqbQE1w6nR4ljmV49+sSiQQKhQKEEFNfIxjnXYGVSqXRtKDx\nREQUTaEnOO81gPEnI+e5urqa+OCyN46bZVkz41zXRbvdhqZpEw+zBI0nIqLoUf46ymqostksMpkM\nPn78OLecbdsoFAqwbXsiyXhdjd7oAI+Zponvv/8evV5vaoIKGj9O+rDrEXf/UrTMd9PiWNd9fdzu\nKMqkHztlDBvuOI5IpVJC13Xhuu7UMs1mU6iqKtrt9txl1Gq1iXm9Xk8oiiLOz88XtmHV+HGSVlts\nABCAkPSLa13c7ij6ZG/DUr5ksru7i7u7O5yfn2N3dxd7e3vI5/NQVRWdTge3t7d4/vz5xIvgj5dx\ne3sLwzDQ6/VQLpeRTCbRarVQq9VgmubcVxGCxhMRUbRI6aJ8zBs9wHVdqKqKQqEwMUzOMvGvXr0K\nVP+y8eyiXA67KNdXH7c7ijLZx86NJLioY4JbDhPc+urjdkdRJvvYKXU0ASIiIlmkjSZA20X+F/6J\niORiF+UK4tBFyW7DqNV1X1/Utzt62thFSUREtAZMcEREFEtMcEREFEtMcEREFEtMcEREFEtMcERE\nFEtMcEREFEtMcEREFEsb+5KJZVmo1+vodrsAgEwmg3K5jKOjo7lxruvCNM3R6NvA/YCmuVzOV71B\n44mIKCKkDs4j7sdey+fzolgsPhj7zXVdUalURKVSmRk7bcw413WFYRjCMIyFdQeN92xgta0dYjtu\nWlzrisd2R0+b7G1Yam29Xk8kk8mZA4s2m02RSqWE4zgT81qtllAURXS73amxuq4LXddn1h00flwc\nDjTxTQRxrSse2x09bbFNcIPBQCSTSXF8fDx1vjeqtqIo4vr6emrs6enpzOXbti0URRGmac6se9X4\nx+JwoIlvIohrXUIAX/21zvB/X3+d2vQmSjEk+9gp7WPL5XIZf/rTnzAYDGYObrq3t4fhcIjb29sH\nZUzTxOnpKWzbxsuXL2fWkc1mAQCfPn16MD1o/GP82PLStbGuyNUX/W2ctk8sBzx1HAfZbBaGYeDy\n8nLpeFVV8fnzZ/z6669zyxmGgevra9ze3j54aCRo/GNMcEvXxroiV1/0t3HaPrEcTaBarQIAvv32\n26VjXddFv99HMplcWDaTyQDAgyQaNJ6IiKJJSoJrNBpQFAWFQmHp2E6nA+D35DOP18XovXqwjngi\nIoqm0N+DcxwHw+EQiqKM7qtdXFyg1+vBdV3c3d0hnU6jUqlgd3d3ajwApNPphXV58ePvuAWNJyKi\naAo9wVmWBQCjLsKzszOcnp7iu+++G5W5vr6Gpmm4uLiYeNHbdd2l67y7u1tbPBERRVPoCa7X6wEA\nUqkULi4u8Ic//AEvXrx4UObo6AiO48AwDLRaLRweHo7mffnyBQB83UPzjCe1oPFERBRNod+DG08W\n8x7Tf/PmDYD7JxlnxQetn4iIno7QE5zX3ec4DnRdn1u2UCjAdV1cX19PxPsx7T5b0HgiIoqm0BOc\nlzT8PEXpvXs2/pj+MklnWjILGk9ERNEUeoIbv/c16wsmnufPnwMA2u32RPyqXY1B44mIKJpCf8jE\nS1p+HvJYVzJa5oGSVePfvXs3+u+DgwMcHBwEqpOIKG5ubm5wc3OzsfpDT3CqqgaK916+9tN96CXG\n8W7JoPGzjCc4IiKa9Pjk//3791LrD72LMp/PA/B3VeaVGb+C8r5A4ufla6+MV+c64omIKJpCT3Dj\nXyfp9/u+YsY/q7W3twfAX4L03rnb399fWzwREUWTlG9RlkolCCFg2/bccn/+858BPPwocyKRQKFQ\ngBBiYYL0rsBKpdLa4omIKJqkJLhyuQxg8Vf6LcuCoig4OTmZGu999msa13XRbrehadrEl1KCxhMR\nUQTJGllV13WhKIpwXXfq/GazKRRFEefn51Pna5omNE2bufx6vS52dnZEv98PJX6cxNUWGsR25Ou4\n1iX/byNaN9nblbTaXNcVqqoKXdcn5vV6PZFMJsXx8fHMeMdxRCqVErVabWr8vOS4jvhxcdj545sI\n4lqX/L+NaN1kb1dSRvT2DIdDGIYBx3FQKpWgqio6nQ6azSb++Z//Gf/0T/80N77f78MwDOzt7aFc\nLiOZTKLVaqFWq+Hs7AyvX78ONd7DEb2Xro11Ra6+6G/jtH1kHzulJjhPv9+HZVlwXReqquLVq1dL\nxXe7XXQ6nY3FM8EtXRvrilx90d/Gafs8iQQXdUxwS9fGuiJXX/S3cdo+so+dUp6iJCIiko0JjoiI\nYokJjoiIYokJjoiIYokJjoiIYokJjoiIYokJjoiIYokJjoiIYokJjoiIYokJjoiIYokJjoiIYokJ\njoiIYumrTTdgnKqq6PV6C8u5rgvTNOE4zmhauVxGLpfzVU/QeCIiigCpo8/NcXJyIhRFWViu2WwK\nVVVFu90eTXNdVxiGIQzDCD1eiHgMBonYDgwa17rk/21E6yZ7u9qK4XJs28be3h4URcGvv/46s5xl\nWSgWi7BtGy9fvpyYXywWAQAfP34MJd7D4XKWro11Ra6+6G/jtH2kHzulptMZSqWSUBRF7OzszCwz\nGAxEMpkUp6enM8vYti0URRGmaa49ftyWrLZAENsrnbjWJf9vI1o32dvVxq/garUaNE1DqVTCTz/9\nNPMKzjRNnJ6ezrz68mSzWQDAp0+f1ho/jldwS9fGuiJXX/S3cdo+T2rAU8dx0O/3cXh4uLBstVqF\noihzkxMA5HI5OI6Dbre71ngiIoqWjSa4s7MzVKvVheVc10W/30cymVxYNpPJAAAuLy/XFk9ERNGz\nsQR3dXWFYrGIZ8+eLSzb6XQA/J585vG6GMevwILGExFR9GzsPbhGo4FGo+GrrPe+WjqdXlh2d3f3\nQcw64omIKHo2cgVXqVRQq9V8l3ddd+k67u7u1hZPRETRIz3B2baNnZ0dvHjxwnfMly9fAMDXPTTP\neFILGk9ERNEjvYvyhx9+8N016QmabJisiIieHqlXcLVaDaenp0vHLdNdOO0+W9B4IiKKHmkJznEc\nOI6Db775ZunYZZLOtGQWNJ6IiKJHWhfl2dkZ/vSnP60U6907W7WrMWj8NO/evRv998HBAQ4ODta2\nbCKiOLi5ucHNzc3G6peS4JZ5521dlnmgZJX48QRHFD9f/fVzbuH7+usUfvqJPSdx9Pjk//3791Lr\nl5LgLi8v0Ww2V473Xr72033oXaWNd0sGjSd6en6BrO9e/vyznERKT0/o9+Asy8L19TV2dnbm/obD\nIYQQD6Z9++23AJZ7+dork8/nR9O8L5isGk9ERNET+hVcoVDAb7/9trDczs7OzPHg9vf3Afi7h+aN\nCO7FAMDe3l6geCIiip6NfmzZr0QigUKhACEE+v3+3LLeFVipVFpbPBERRU8kEhwAlMtlAPddnrO4\nrot2uw1N0ya+lBI0noiIomXjA5565nVReryuRm90gMdM08T333+PXq83NUEFjfdwwNOla2NdkauP\ng6vS+j2pAU89V1dXAAAhBK6vr2eWazabcBwH5+fnE/Mcx8Hp6Smq1erM5BQ0noiIIkRsSKvVEqlU\nSiiKIhRFETs7O2JnZ2f0/6lUSrTb7Yk4x3GEpmmiXC4L27aF4ziiXq8LVVXFxcXFwnqDxgshxAZX\n29oAEICQ9GNd0atPbl30NMj+t96aLspldbtddDoduK4LVVXx6tUrafHsoly6NtYVufrYRUnrJ/vY\nGdkEt0lMcEvXxroiVx8THK3fk7wHR0REtG5McEREFEtMcEREFEtMcEREFEtMcEREFEtMcEREFEtM\ncEREFEtMcEREFEtMcEREFEtMcEREFEtMcEREFEtfyazMdV18+PAB7XYbjuMgnU4jn89jf38fJycn\nSCQSvpZhmuZo5G3gfjDTXC7nuw1B4omIKCJkDVvQarWEpmni+vpaDIdDIYQQtm2LSqUyGiLHNM25\ny2g2m0JV1QfD6LiuKwzDEIZhLGxD0HiPxNUWGsR46JV41hXnvw2b3h1IEtn/1lJq6/V6QtO0mfNt\n216Y5FqtllAURXS73anzdV0Xuq7PrCNo/Lg47JBxPljGs644/23Y9O5AksQywRmGMbpqm6VWq42S\nnOM4D+YNBgORTCbF6enpzHgvSU5LkEHjH4vDDhnng2U864rz34ZN7w4kiex/69AfMvHudT179mxu\nuTdv3ozuwdXr9QfzGo0GhsMhyuXyzPhcLodMJoNqtToxL2g8ERFFT+gJzrZtXF1doVgsYjgczi27\nv78PALi6unowvVqtQlEUvHz5cm58LpeD4zjodrtrjSciouiRdgVnWRba7fbcsru7uwCAwWAwmua6\nLvr9PpLJ5MK6MpkMAODy8nJt8UREFE2hJzhd15FMJqGqKgqFwtyyd3d3AIB0Oj2a1ul0APyefObJ\nZrMA8OAKLGg8ERFFU+jvweVyuVHiWsS72svn8xPTxpPeLN4V4Pg7bkHjiYgomrbmSyau646unHRd\nfzB9WeMJNWg8ERFF09YkuEajAQBQVRWvX78eTf/y5QsA+LqH5hlPakHjiYgomrYmwXlPOj5+TD9o\nsmGyIiJ6mrYiwVUqFfT7fbx9+xavXr16MG+Z7sJp99mCxhMRUTRtPMHZto3z83MYhoEPHz5MzF8m\n6UxLZkHjiYgomqSOJvCY67o4PDyEYRgz3z3z7p2t2tUYNH6Wd+/ejf774OAABwcHa10+EVHU3dzc\n4ObmZmP1bzTBHR4eolgsrv3F6mUeKFk1fjzBERHRpMcn/+/fv5da/8YSnK7ryGazC5Ob9/K1n+5D\n7yptvFsyaDwRhe0rKIoipaavv07hp594K+Kp2EiCK5fLSKfTvq7cvC+Q+Hn5etqL4kHjZfjll1/w\nD/9QxN3dz1Lq+5u/2eiFO9EjvwAQUmr6+Wc5iZS2g/QjXa1Ww+fPn/Hjjz/OLHN+fo43b94AAPb2\n9gD4u4fW6/UA/P7R5nXEy/Df//3f6HT+H3799f9Kqe/v/u7/SKmHiGijZI7N02w2RbFYXFju8cCj\nuq5PHSfusUKhIBRFEf1+f63xj617tf3Xf/2X+Oqr/yVt/K1EIh/rscXiWVec/zaOPfdUyF7/0l4T\nsG0bpmnOvXID7q+0Hn8Y2RvHzbKsuXHtdhuapuHFixdrjSciouiRkuAcx8HZ2Rnq9Tpc1535s20b\nlUpl9GCI5+joCPl8fmIg1HGNRgOKoqDZbE7MCxpPREQRFPYl4mAwEJlMRiiK4vt3fX09sRzHcUQq\nlRK1Wm1iXq/XE4qiiPPz85ntCBo/bt2rjV2UrGv76otvXbQ5std/6A+ZfPjwAf1+f6nHgKc9xbi7\nu4vb21sYhoFer4dyuYxkMolWq4VarQbTNB98pHnd8UREFC3KX7NqpHS7XXQ6HbiuC1VVJ75fGXa8\noihY52r7y1/+gmTyf+OXX/6ytmXOk0hoGA5tQNKj2YDCuiJXX3zriuAhLzbWfexcWF8UE9ymMcEt\nK74HSya46NXFQ97myE5wG//YMhERURiY4IiIKJaY4IiIKJaY4IiIKJaY4IiIKJaY4IiIKJaY4IiI\nKJaY4IiIKJaY4IiIKJaY4IiIKJaY4IiIKJZCH01gG7muC9M04TjOaFq5XEYul9tgq4iIaJ2eXIK7\nurrC2dkZTNPE27dvAQDD4RDfffcdgPuBT4mIKPqe1GgClmWhWCzCtm28fPlyYn6xWAQAfPz4ce5y\nOJrAsuL7ZXqOJhC9up7QIW/rcDSBkLiuC8MwUC6XpyY3AKhWq7AsCxcXF5JbR0RE6/ZkElyj0cBw\nOES5XJ5ZJpfLIZPJoFqtSmwZERGF4ckkuGq1CkVRZl69eXK5HBzHQbfbldQyIpLnKyiKIu337Fl6\n03/wk/YkEpzruuj3+0gmkwvLZjIZAMDl5WXYzSIi6X7B/f0+Ob+ffx5I+rtomieR4DqdDoDfk9c8\n2WwWAHgFR0QUcU/iNQHvfbd0enF3we7u7oMYIqLV3XeJyvD11yn89NOdlLqi4kkkONd1l465u+OG\nQkRBeV2i4fv5ZzmJNEqeRBflly9fAMDXPTjPKkkxKn755edNNyGgm003IICbTTcgoJtNN+CJu9l0\nAyLlSSS4OCerVTDBbdLNphsQ0M2mG/DE3Wy6AZHyJBLcMt2Nfu7TERHR9nsS9+CWSVqbu/f2G4DP\nUmoS4lcp9RARbdKTSHDevbd1dVWqqhrSk1G7ISxz0v/8j/dfMm9Kr7uu9xLrmmeVuua1PYz61l1X\nkPYvW1cYZD+MIW/bl/XE5qpUVZVa35NIcKuY90DKp0+fJLaEiIhW8STuwXkvb/vpfvSu8ngvjogo\n2p5EgvO+YOLn5W2vTD6fD7VNREQUrieR4Pb29gD4uwfX6/UAAPv7+6G2iYiIwvUkElwikUChUIAQ\nAv1+f25Z7wquVCrJaBoREYXkSSQ4AKNx4CzLmlnGdV20221omoYXL15IahkREYVBEU9o/Pa9vT38\n+uuv+OMf//jgfly5XEYul4Npmvj+++/R6/UCJTjXdfHhwwe02204joN0Oo18Po/9/X2cnJwgkUis\nvFzTNKe2XQZVVUdduKuQ3X7LslCv10cjQ2QyGZTLZRwdHa20PFntt20bpmmi0+nAdV1kMhnouh5o\n2/GUy2VBbOuqAAALg0lEQVRks1m8efNmTa393aa3TwpXmNuOZ937LMQT8q//+q9iZ2dHnJycjKa5\nrisMwxD/+I//KBRFEefn54HqaLVaQtM0cX19LYbDoRBCCNu2RaVSEYqiCEVRhGmaSy+32WwKVVVF\nu92eaLthGIHa7MfJyYlQFGXleJnt7/V6Ip/Pi2KxOFFfpVIRlUpl6WXKan+pVBLHx8ei2+2OpjmO\nI2q1mkilUittO71eT9TrdZHJZISiKOLs7Gxt7fXIWj8nJyeiVqutbXmewWAg3r59KzRNE6lUSqiq\nKgzDELVaTbiuu7Z6wmr/PJlMZuVYGduOV8+691khhHgyCa7VaglFUcS///u/C03TRLlcFrZtC8dx\nRL1eF3/7t38r/v7v/z5QHb1eT2iaNnO+bdsrJTmv7eMHvXG6rgtd15dur1+3t7dCURSxs7OzUrzM\n9vd6PZFMJmeeqDSbTZFKpYTjOL6XKav9pVJJXF9fz5zvOI5QFEVcXV35Wl4+nxepVEpomiYuLi5G\nJ1nrPkiFvX7CPsiGdVLqkZUkpln1xFTWtiNEOPus50kkuMFgIJLJpDg9PR1Ns21bmKYparWauL6+\nHiWfIBuyYRijHWSWWq022mn8/INNa/tj62j7PKVSaeUEJ7P9Xl3Hx8dT5/d6vdG6n5dIpi0z7PY3\nm82Z7R5nmqZIpVIr1VGtVtd+kApz/cg4yIZ1UiqE3CQxTdAT03FhbDtChLPPjnsSCa5er889w/So\nqipUVV2pjl6v57srJplMCkVRfF12y2j7PNVqVViWJZLJ5Eo7isz2n5yciJ2dnbknGZqmiWw2u/BE\nxCOr/YVCQVxcXCwsNxgMfLVnmjAOUjL/fcNofxgnpbOElSRmCXJi+lhYbQ9jnx33JJ6irFarUBQF\nL1++nFsul8vBcZzRDc5l2LaNq6srFItFDIfDuWW9d+yurq4WLldG22dxHAf9fh+Hh4crL0NW+x3H\nwcXFBUqlEp49ezazXKfTwX/+53/OLTNOVvv7/T4Gg8HCcsuMaSjDJrfPoLyHYRZtC2/evBk93FOv\n10Nv1zrUajWcnp4GfigpTGHts+Nin+Bc10W/3/d1YPC+eHJ5ebl0Pd7OYlkW2u323LK7u/cfVV50\nQJPV9lnOzs5QrVZXjpfZfq+d33777Urx08hsfyKR8HXw9LazRQlFhk1vn0GFdVK6aes4MZUhjH32\nsdgnuE6nA+D3HWwe75uVq5xl6rqOZDIJVVVRKBTmlvW+ibnoe5ey2j6Nt+Ovctbkkdn+RqMBRVEW\nrvtlyGz//v4+HMdBsVicW65er8MwjJXqWLdNbp/rEMZJ6TYIemIqSxj77GOxH03A24j9fDzZ24j9\nfLPysVwu53ssOb/fu5TV9mkajQYajUagZchqv+M4GA6HUBRllJAvLi7Q6/Xgui7u7u6QTqdRqVRG\n9WxT+wGgUqnANE1YloV0Oo1mszlxBm5ZFi4uLvD58+eV6li3TW6f66DrOn744Qc8f/58bSelm7aO\nE1MZwtpnH4v9FdwqY8CFOeip67qjs1hd1xeWXdY62l6pVFCr1QIvR1b7va/TeF1lZ2dno4PXv/3b\nv6HRaEDXdWiahuvra9/Llbn+d3d30Ww2R/Xquo7T09PRfMuycHZ2Btu2t+bgtW371rK8k1I/93ei\n8hH2RqOB169fb7oZC4W1zz4W+yu4L1++AFju5vy6BkadxrsqUlV14Ya4ibbbto2dnZ21fKpMVvu9\nr6ukUilcXFzgD3/4w0T7j46O4DgODMNAq9XydX9C9vo/OjrC7e0tDMOA4zijK7pCoQBVVUddgtti\n2/atsCxzUrpJ6zoxlSGsffYxXsFJ5j115qePfBNt/+GHH/Dhw4e1LEtW+8frsW175gMY3ieG/N7D\n2sT6z+Vy+PTp0+hKwUt0XjLZJtu2b4VlmZPSTVnniakMYe2zj8U+wS3TJRJ2/3qlUkG/38fbt2/x\n6tWrheVlt917tHhdZLXfq8dxnIVn2IVCAa7r+ur22MS247ouDMOArutoNpujq6NarYa9vb2Fo2HI\ntE37VpiWOSndlHWemMoQ1j77WOwT3DI7Vpj3B2zbxvn5OQzD8L0hymy74zhwHAfffPNNoOWMk9V+\nrx4/T2R5H/7187i67G3Hsizs7e3h/PwcHz58wNHREfr9/uhvsm0bqqpuzZOI27JvhWnZk9JNWPeJ\nqQxh7bOPxT7BeWfAm+xOcV0Xh4eHMAxjqX8kmW0/Oztbe/+9rPaP3wNa9LDA8+fPAWDhY+Hjy5Wx\n/m3bxvHxMSzLetDNlEgk8PHjx9EDKACgadrC97Zk2IZ9K0yrnJTKFsaJqQxh7bOPxT7BrWLdX4s4\nPDxEsViU8pLrKm3fpkeLV2m/twP4iQ37oLxK+71uyVqtNvMeytHREXq93uids++++y5IMzdm277E\nMsuqJ6WyhXFiKoOsfTb2Cc57wdRPF4m3Atd5v0DXdWSz2ZV2Elltv7y8DOXmuaz2q6q6dIwfstpv\nmiZc1134b7C7u4vb21skk8lAj06vy6b3rTDJPCld1TadmC4rrH32sdgnOO+M188Lput+16VcLiOd\nTq+8k8hou2VZuL6+xs7OztzfcDiEEOLBtEWf2JG17r0YP2d4XpllPi8VdvtbrRb29vZ8lU0kEqhW\nqxBCbPyF703uW2EKclIqU1gnpjKEtc8+FvsE5x04/KxI790M77tzQdRqNXz+/HnuTnJ+fj53GTLa\nXigU8Ntvvy38Afc3hMenLToAyFr341868PuUoZ/PS8lq/2Aw8NUez7Z8Y3BT+1aYgp6UyhLmiakM\nYe2zj8U+wSUSCRQKBQghFq5I7yyzVCoFqvPq6grtdhs//vjj3HKtVmvu/E20fZ1ktr9UKkEIAdu2\n55b785//DMDfB15ltT+dTo8SwDI2/c5T1LfPx9ZxUipLmCemsoSxzz4W+wQH3J+VAb9/HmYa13XR\nbrehadrMA4dpmgvvfdi2DdM0FyY313V9nZHIbHsYZLXfq2fRzmtZFhRFwcnJyYKWP1xumO03DAOW\nZfl+MvLq6kraB5f9rveg62fT1nVSSr/b1D77wPJD1EWTpmlzR+6t1+tiZ2dH9Pv9qfO9wQMVRRFX\nV1dTy/R6PaHrunAcRwwGg5m/29tbcXJyMnOI9k20fZEgAyfKar+u60JRFOG67tT5zWZTKIrie73L\nbL+qqkLX9YVtub29FalUamZd8yw7aKXf9R50/fi1yqCb9Xp94TZ/e3vra90PBgNRLpd91/1YWO2f\nZ1MDnm56n/U8mQTnOI5IpVKiVqtNzPOGRZ+3ElVVHf2DTSs3GAxEJpMZlfHz8zsEe9htX8TbyFZN\nkLLa77ruzETR6/VEMpkUx8fHW9l+rw5d1+fu7Kqqina7vfTfMBgMRD6fF4qiCE3TZtaxbLvH277q\n+vErjINsWCelstq/yDoSXJjbTlj7rOfJJDgh7ndETdNEuVwWtm0Lx3FEvV4XqqqKi4uLubGWZYlU\nKjVz6PS3b9+ONia/v2XOaMNs+zStVkukUqnRRuq12fv/VCq11IFWVvtd1xW6rgtVVUWlUhGmaYqT\nkxORSqUCHZxktb9Wq40S3Xj7NU0TxWJxqW2mVCpN/Td8/O9YLBYDtzvI+vEjjINsmCelMtq/SJAT\nU5nbTlj7rBBCKEIIsXzHZrR1u110Oh24rgtVVbf2EzzTRLntgLz29/t9WJa19npktf9xPYVCIRLv\nO61z/RiGgXa7PXpKU1GU0TzvsJVMJrG/vz/13lm73YZhGHj+/Dlub28n1l+lUsH5+fmD5S7S6/V8\n30cMu/3TWJaF4+PjiTrH67u6utrKL5+Esc8+yQRHRETx9ySeoiQioqeHCY6IiGKJCY6IiGKJCY6I\niGKJCY6IiGKJCY6IiGKJCY6IiGKJCY6IiGKJCY6IiGKJCY6IiGKJCY6IiGKJCY6IiGLp/wOqqVRP\njbwcaQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1128df550>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Accuracy Validation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's output our predicted ratings as the _mean_ rating. Let $I$ be the set of item-user pairs rated, with true ratings $r_{ij}$. Recall $z_{ui}$ is the logarithm of the rating. The standard error metric is _root-mean squared error_ (RMSE):\n",
      "\n",
      "$$\\text{RMSE}(I) = \\sqrt{\\frac{1}{|I|}\\sum_{(i,j)\\in I}\\left(\\exp z_{ij}-r_{ij}\\right)^{2}}$$\n",
      "\n",
      "Let us calculate the RMSE for both training and test sets. In general, if the RMSE is way better for training set, then we've overfit."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#trainRMSE = np.sqrt(np.mean((np.exp(rmean[users-1, items-1]) - ratings) ** 2))\n",
      "#testRMSE  = np.sqrt(np.mean((np.exp(rmean[testUsers-1, testItems-1]) - testRatings) ** 2))\n",
      "trainRMSE = np.sqrt(np.mean( (meanExpZ[users-1, items-1] - ratings) ** 2))\n",
      "testRMSE  = np.sqrt(np.mean( (meanExpZ[testUsers-1, testItems-1] - testRatings) ** 2))\n",
      "\n",
      "print trainRMSE\n",
      "print testRMSE"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.686684550069\n",
        "0.876480474553\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The Competitor, Nonnegative Matrix Factorization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Again, I jiggled the sparsity (beta0) knob; these were the best results I got.\n",
      "\n",
      "Note exactly a fair comparison, since the NMF code is far better vectorized."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R = coo_matrix((ratings, (users - 1, items - 1))).todense()\n",
      "skm = ProjectedGradientNMF(n_components=D, sparseness='data', beta=70)\n",
      "W = skm.fit_transform(R.T)\n",
      "\n",
      "reconstruct = np.dot(skm.components_.T, W.T)\n",
      "\n",
      "nmfTrainRMSE  = np.sqrt(np.mean( (reconstruct[users-1, items-1] - ratings) ** 2))\n",
      "nmfTestRMSE   = np.sqrt(np.mean( (reconstruct[testUsers-1, testItems-1] - testRatings) ** 2))\n",
      "\n",
      "print nmfTrainRMSE\n",
      "print nmfTestRMSE\n",
      "\n",
      "#print skm.reconstruction_err_\n",
      "#pint np.linalg.norm(reconstruct - R)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.907722163866\n",
        "1.39713801489\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}