8 - Anomaly Detection & Recommender Systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomaly Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(307, 2)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = loadmat('machine-learning-ex8/ex8/ex8data1.mat')\n",
    "X = data['X']\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f6052eca2e8>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAHVCAYAAADl4K3UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X9w3PV95/HX2/I6XTtcZIpKsYCYy2TMhFDsZA/oOdeLaWITkxCVtBCOpqTp1c1N0gmZnHp20ynkx0x05xJyvXSaIQmT5EKpSzAqKSTGFzOTC1dSZGRjHOxCEgNeU1BqZAJegiy/7w99v/JXq+/3s9/VrnZX0vMxo9Hud7+7+9mv5dVLn31/3x9zdwEAAABIt6jdAwAAAAA6GYEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEDA4nYPIM0ZZ5zhK1eubPcwAAAAMI/t3r37Z+7eU2u/jgzMK1eu1NDQULuHAQAAgHnMzJ7Ksx8lGQAAAEAAgRkAAAAIqBmYzewcM3vAzH5kZvvN7GPR9pvMrGxme6KvjRn3v9zMDprZk2a2udkvAAAAAJhNeWqYT0j6hLs/YmanSdptZjuj225x97/IuqOZdUn6K0nvlHRY0sNmdo+7/6jRgQMAAACtUHOG2d2fdfdHoss/l/S4pN6cj3+xpCfd/Sfu/qqkv5X03pkOFgAAAGi1umqYzWylpDWSfhht+qiZPWpmt5nZ8pS79Ep6JnH9sDLCtpltMrMhMxsaGRmpZ1gAAADArMkdmM3stZLuknSDu78o6a8lvUHSaknPSrq5kYG4+63uXnL3Uk9PzXZ4AAAAQEvkCsxmVtBEWL7d3bdLkrs/5+7j7n5S0pc1UX5RrSzpnMT1s6NtAAAAwJyQp0uGSfqqpMfd/fOJ7WcldvstSY+l3P1hSW80s/PMbImk90u6p7EhAwAAAK2Tp0vGWkkfkLTPzPZE2/5U0rVmtlqSSzok6Y8kycxWSPqKu2909xNm9lFJOyR1SbrN3fc3+TUAAAAAs6ZmYHb3H0iylJvuy9j/iKSNiev3Ze0LAAAAdDpW+gMAAAACCMwAAABAAIEZAAAACCAwAwAAAAF5umTMe4PDZW3dcVBHRita0V1U/4ZV6luTd/VvAAAAzGcLPjAPDpe1Zfs+VcbGJUnl0Yq2bN8nSYRmAAAAUJKxdcfBybAcq4yNa+uOg20aEQAAADrJgg/MR0YrdW0HAADAwrLgA/OK7mJd2wEAALCwLPjA3L9hlYqFrinbioUu9W9Y1aYRAQAAoJMs+JP+4hP76JIBAACANAs+MEsToZmADAAAgDQLviQDAAAACCEwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAE1AzMZnaOmT1gZj8ys/1m9rFo+1YzO2Bmj5rZ3WbWnXH/Q2a2z8z2mNlQs18AAAAAMJvyzDCfkPQJd3+TpEslfcTM3iRpp6Q3u/uvSfpnSVsCj7HO3Ve7e6nhEQMAAAAtVDMwu/uz7v5IdPnnkh6X1Ovu97v7iWi3hySdPXvDBAAAANqjrhpmM1spaY2kH1bd9CFJ38m4m0u638x2m9mmwGNvMrMhMxsaGRmpZ1gAAADArMkdmM3stZLuknSDu7+Y2P5JTZRt3J5x17e5+1skvUsT5Ry/kbaTu9/q7iV3L/X09OR+AQAAAMBsyhWYzaygibB8u7tvT2z/oKR3S7rO3T3tvu5ejr4/L+luSRc3OGYAAACgZfJ0yTBJX5X0uLt/PrH9ckl/IulKdz+ecd9lZnZafFnSekmPNWPgAAAAQCvkmWFeK+kDki6LWsPtMbONkr4o6TRJO6NtX5IkM1thZvdF9z1T0g/MbK+kf5J0r7t/t/kvAwAAAJgdi2vt4O4/kGQpN92Xsk3ufkTSxujyTyRd1MgAAQAAgHZipT8AAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAATUDs5mdY2YPmNmPzGy/mX0s2n66me00syei78sz7n99tM8TZnZ9s18AAAAAMJvyzDCfkPQJd3+TpEslfcTM3iRps6TvufsbJX0vuj6FmZ0u6UZJl0i6WNKNWcEaAAAA6EQ1A7O7P+vuj0SXfy7pcUm9kt4r6evRbl+X1Jdy9w2Sdrr7UXd/QdJOSZc3Y+AAAABAK9RVw2xmKyWtkfRDSWe6+7PRTf8i6cyUu/RKeiZx/XC0Le2xN5nZkJkNjYyM1DMsAAAAYNbkDsxm9lpJd0m6wd1fTN7m7i7JGxmIu9/q7iV3L/X09DTyUAAAAEDT5ArMZlbQRFi+3d23R5ufM7OzotvPkvR8yl3Lks5JXD872gYAAADMCXm6ZJikr0p63N0/n7jpHklx14vrJf19yt13SFpvZsujk/3WR9sAAACAOSHPDPNaSR+QdJmZ7Ym+NkoakPROM3tC0jui6zKzkpl9RZLc/aikz0h6OPr6dLQNAAAAmBNsovy4s5RKJR8aGmr3MAAAADCPmdludy/V2o+V/gAAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABi9s9AAAAMLsGh8vauuOgjoxWtKK7qP4Nq9S3prfdwwLmDAIzAADz2OBwWVu271NlbFySVB6taMv2fZJEaAZyoiQDAIB5bOuOg5NhOVYZG9fWHQfbNCJg7mGGGQDmGT5+R9KR0Upd2wFMxwwzAMwj8cfv5dGKXKc+fh8cLrd7aGiTFd3FurYDmI7ADADzCB+/o1r/hlUqFrqmbCsWutS/YVWbRgTMPZRkAMA8wsfvqBaX41CmA8wcgRkA5pEV3UWVU8IxH78vbH1regnIQAMoyQCAeYSP3wGg+WrOMJvZbZLeLel5d39ztG2bpPjdt1vSqLuvTrnvIUk/lzQu6YS7l5o0bgBACj5+B4Dmy1OS8TVJX5T0jXiDu18TXzazmyUdC9x/nbv/bKYDBADUh4/fAaC5agZmd/++ma1Mu83MTNLVki5r7rAAAACAztBoDfN/kPScuz+RcbtLut/MdpvZpgafCwAAAGi5RrtkXCvpjsDtb3P3spn9iqSdZnbA3b+ftmMUqDdJ0rnnntvgsAAAAIDmmPEMs5ktlnSVpG1Z+7h7Ofr+vKS7JV0c2PdWdy+5e6mnp2emwwIAAACaqpEZ5ndIOuDuh9NuNLNlkha5+8+jy+slfbqB5wMApBgcLtMVAwBmUZ62cndIerukM8zssKQb3f2rkt6vqnIMM1sh6SvuvlHSmZLunjgvUIsl/Y27f7e5wweAhW1wuKwt2/dNLoddHq1oy/Z9ktSW0Ex4BzAfmbu3ewzTlEolHxoaavcwAKDjrR3YlbqyX293UQ9ubm0Do+rwLk0smvK5qy4kNAPoSGa2O886Iaz0BwBz2JGUsBzaPpu27jg4JSxLUmVsXFt3HGz5WACgmQjMADCHregu1rV9NnVSeAeAZiIwA8Ac1r9hlYqFrinbioUu9W9Y1fKxdFJ4B4BmIjADQJMMDpe1dmCXztt8r9YO7NLgcHnWn7NvTa8+d9WF6u0uyjRRu9yumuFOCu8A0EyNLlwCAFB7u1X0rentiJPq4jHQJQPAfENgBoAmCJ3wtpACY6eEdwBoJgIzAKjx/sHtOOGt3jF3ao/kTh0XAMQIzAAWvGaUU6zoLqb2Q56tE94Gh8vqv3Ovxk5O9NIvj1bUf+deSeljbmXJSD0BuNMWXgGANJz0B2DBa0b/4Faf8HbTPfsnw3Js7KTrpnv2p+7fqh7JcQAuj1bkOhWAs06ApHczgLmAwAxgwWtGOUWru1WMVsbq2t6qkpF6AzC9mwHMBZRkAFjwmlVO0SknvK0d2DWtHKJVJSNZQbc8WmnruACgEcwwA1jwOrV/cKiv8/Klhcz7pZVDzMZrTBtfVtC1OsYV79uqXtYAUIu5e+29WqxUKvnQ0FC7hwFgAemUTg3xOMqjFZmk5Dt0sdA1WeYxOFxW/7f2amy89nt4b3dRD26+rO6T8UL7Vp+sF4/vfW/t1V27y1O2V7+OtHHddM/+zHKS5UsLuvE9F3TE7D2A+cXMdrt7qdZ+lGQAgNpfTjE4XNanvr1fLxw/FRqrQ2ayr3PaIiFppQ3SqTKJvK8xT+eKrFrlBw6M6HNXXVjXuCTpFydOZo7nheNjdM4A0FYEZgBos7TZ2izJkFkdgNcO7JpRPXD1DK+ZVP3hY/LEvXgGPGt89Y4rLXxXS/6xkGemvFM+MQAwP1DDDACzIFR/XC1PYIyFwu9M6pTjfs7JcoisSr14pjkrLGeNr9a48nbEODJaydW2rt7WdgDap573ynaihhkAGlQ9m7nu/J5pdbzJ+uNq522+N7XGt1roMbLGEofSrNnWrNnfNF1mGg/8zgiNLzmu7qUFuUvHKmNa0V3U8VdPTClFydIbhfG08XaZ6aR78PHimmkAnSHrXIjZbMlZLW8NM4EZABqQ9oZf6yS3aqHQGj9W7wzLCmr9QqonrIdmwbvMdPPVF9UcX9p4CotMMgVPYCwsMm39nYv08W17co03jUn66cAVM7w3gGbLeu9r5R+3eQMzJRkA0IC0coqsQJdVepBWsiBJ3cWCbrlmtb5wzWpJ0se37Zn8yLLWx5jx7Tds2xNcSCRvv+N4UZYsJ91zhfm04zV20rVsyeLJRV+WLy1okVXdMbreSH/m7kArPgCtN5cWLiIwA0AD6nlj715aSA25aasEfuGa1dpz43pJmlaP23/nXvV/a29mjW6yhrfWuPs3rJqY4Q3o7S6qb01vsBY6TxgdHC5njulYZUwPbr5MPx24QkuXLFbVqt8aG3dt3XEw84+LPDrwA1VgQcv6A7gTFy6iSwYANCCrbVp1WUahy/TSK6dqa6tbtWW1fMuaka1WGRvXJ/5urz6+bY8W1ag1lqRFZhocLk8+Z1Yf5OTJeX1rerVl+6OqjE1vAVcrjP7Z4D7d/tDTmbcnf0GGZp2q2+nlea2xYxl9ngG0R/+GVaklY+1eNCoNM8wA0ICsDhDXXXrulBnjZUsWTwu6ydKILPXMYI+7y6PvefaNZ6X71vRqz43rdWjgCn3hmtVTxl198s0rKWFZCofRweGybn/o6cxSlepfkFmz1a6JEyTjmeafDlyhm6++KHWlwDSdOGsFLGRpn6618oS/ejDDDAANSFtAJO3kvPM235t6/7RAnOwoUc8Mar2SvY1jWTPd8ZiyRhIKo6H7SZr2CzL0cpMlKPF44+eo1aWkE2etgIWu3YtG5UVgBoAG5XnDzyrdqA6a1V0k0sJynq4SeeWZwa61sEoyjKa1tQs9R1wfnZSndKJ61cPqxyi9/nQWLgHQNARmAGiBvLV6WYuYJGuiX/tLi/Wms07T//vx0Vwt1haZ9EuLF+l4SjlFHNhDK+OFFlZJtrvLWlK7e2khs8/yuvN7pm0L7Z8UCuJzZdYKwNxAYAaAhJkuqVzrflmlG9JEL9J4W1YXiWQwfuH4mB788dHcr+mkS78YdxUW2ZQ66jiwpwXdj2/bo6GnjuqzfRdmBlOTpvRKTQvWlbFxvWbxosw+znftLqv0+tNzl2QkUZMMoFUIzAAQyZohlRQMzXnvVz3rOThcVv+39k6WVuRdcW8mxk+6/s3SgpYuWTwtsH/i7/ZOK/1wSd986Gl986GnM1f4qw6sWcH6WGVMt1yzOvV50uqo85RkmNJnpwFgNhCYASCSNUMat2vLmnHOul91EKz2qW/vb0odcl4vHB/Tje+5YHJMcdCvdVJh2u1p5SS16rSznqc6aIdm2mOuU7PTUu2TLgGgEQRmAPNSvaUVoUU14qCXNXOcNbNaHq1MKbdYd36PHjgwMnk9VKfb212c3O/4qydy1fTmkRz/p769P7jcdbUuM510zzyeWXXa687vmXzeNMXC1A6naY+Tttx4ZWxcfxr1hY5vy/upAADUg8AMYN6pt7Qi3j+PtJnj0OIl8fbyaEXfTCzcUWsGNVkbfN2X/7GumuWQZO/nekP4uPtkkN+646CGnjo65Q+A/g2r9LmrLpz2h0ropEFJqpyYejJiWr131vFKO5Exz+w+ANTDvAPXCi2VSj40NNTuYQCYo9YO7EoNWL2JAFcd6OqpHzZJPx24YvJ6Wtu1tBnRvLqLhcllsf9scN+UoJ3H8qUFuSt15b54bHnKHtLuF3pNxUKXPnfVhZLyhd2k5Ix62ux11r9paKzJfyMASGNmu929VGs/VvoDMO+ESiS2bN+n8mhlygIY9QbH6pPd0laramQq4oIVp01evuOHz9R9/9HjY9pz43r1ZnSRWGQ2oxMMa72myti4Pnn3PvXfuXfKMc6j+t9kcLg85fa0FRVD6KABoJkIzADmnayw1GWWenJel6Uvprx8aSF12es8K8YtW5I/3FV78MdHdcGff1fnbb53Rqv8reguanC4rJd/cSL19jyPWViUvcR0yMuvjk9bArxeaUuGp/1RkvHPJkms6gegqahhBjDvZJ18llVHO+4+7fZioUs3vucCSbU7MKTVTDfq5Vdrn4zXXSzoFydOThv30iWLdMO2PQ09/7jPvKQkS1Z7ujRpnxJUt+VbmbHceLwvADQLgRnAvJO1SEhWrXJWbXP8OLXCV62T2mbLTVdODfTdSwt66ZUxPfH8yw0/doOTxKluvvqi4BLbSXlKKrICeNYnBgAwUwRmAPNS1tLIWctTN7KUcmiJ5tlS6JoIhfG441nulKYRM1JrNtiiMbyas4/08qUF9a3p1dBTR2uexJi37CVrfDMpYwGAEGqYASwYaXWwn7vqwoY/vq/nBLPuYmHyub9wzWr97qXnzug5x8Z9Sp1vs2e5r73knOBJdi5p6ZLFuU7E61pkk+UtDxwYCe5bz79J1kmNWduzDA6XtXZgl87bfK/WDuyadsIhADDDDGBBSVueOrm4SJ5V4qoXRVl3fo/u2l3OFViPVcamPE/fml79dOSlGfVZTpaXNHNZ7WJhkT7bN9Ee7vaHns6sZY6XvI6PxaKMWenTXrN48piGZuO/cM3q3H+8ZJ3UmHd2Ovk4M1kOHcDCQmAG0HT1rrLXrvGkhaWPb9ujG7btmaxrznOC3127y3rfW3snF/Ewy64BTrZOiz3y9LEZva64VndwuNxQ3+dqr0R1HQ8cGAk+5oru4pSSkKwTDZP9oLP6MsclG3mk9b2OHyO59HceM13WHMDCQmAG0FSdNmMXGk9aWMpaYjkO3WlhrzI2rgcOjEyuzjc4XFb/t/ZqLFDfWxkb1w3b9tTVOaJafL+tOw42taNFXGISmg1OzuTWWinRon361vRmdjCJSzbyyCo/Wbpkcd0/Y1mvsR116QA6V80aZjO7zcyeN7PHEttuMrOyme2JvjZm3PdyMztoZk+a2eZmDhxAZwrN2HXaeGqFoni/OBCGyh6Sj9W3plfX/LtzcnVraOQEtd4cwbZehUWm46+e0Hmb79WijPF3mU2pM65VP+3RPlJz6sibGXKz6s9Z+ARAUp4Z5q9J+qKkb1Rtv8Xd/yLrTmbWJemvJL1T0mFJD5vZPe7+oxmOFcAc0GkzdqHx5Fm2+choJdcJdcmANThc1l27y7PereG5F1/R6k/d37TZZZN0UtILxydKKNLGHy9/nQy4ef5tq/+gaOTThqx/t5mE3KwZbxY+AZBUc4bZ3b8vqf6zUaSLJT3p7j9x91cl/a2k987gcQDMIZ02YxcaT57llld0F2sGwuqA1aq+zCdO+pT64Ea5pPGU4mszZc4GDw6XM2eik15XLDRtnGn/bjMNubPVOQXA/NJIDfNHzez3JA1J+oS7v1B1e6+kZxLXD0u6JOvBzGyTpE2SdO65M2uzBKD9Om3GLjSeZElB1kzzuvN79MCBkczbe6MuGVt3HNTHt+3JNWs917hLhwaumLY9LlXJM5P+8qsnJuuYG5W1MM1MH7vRGW8A8595jjc6M1sp6R/c/c3R9TMl/UwTExKfkXSWu3+o6j6/Lelyd//P0fUPSLrE3T9a6/lKpZIPDQ3V90oAdIy50iUjae3AruAqgHlXqJOU2bFiJo/VKdICc9Yxy9LbXZw8MTKk035+AMxfZrbb3Uu19pvRDLO7P5d4oi9L+oeU3cqSzklcPzvaBmCe65QZu+rgdUugz2+o1jnPTHSSa3poTpvVPjJaaWp3i9nSnVFOUW9dep79O63LCgBIMwzMZnaWuz8bXf0tSY+l7PawpDea2XmaCMrvl/SfZjRKAKjD4HBZn/r2/smT16TawavWiWTxHwF5Z1VdEzOqyVlSSVMWSbnlmtWZvYs7RWGR6aYrJ1q+Jf8A6V5af01ynjr2evsiMxsNoBVqBmYzu0PS2yWdYWaHJd0o6e1mtloTvxMOSfqjaN8Vkr7i7hvd/YSZfVTSDkldkm5z9/2z8ioAIJK1qIV0qvfx1h0HpwWrvLXXeWdVly8tTCk/GBwuq//OvRqLTqorj1bUf+deLVvSpZdf7czyjOTiLdXHNfnHSFJWOYpJuerY6+mywmw0gFapGZjd/dqUzV/N2PeIpI2J6/dJum/GowOAhDyziXk6VKQFq7wnkr2uWMjVmeLY8TGt+fT9Gj0+sRT2Cy//YjIsx8ZOuhY3ufVcvBBKb3dRx189kRpsu4sFLXvNYpVHK6kB1yRdd+m5k8tjS/mOa2gRFld2iE3+u2Ytr502O80qfQBahZX+AHSk6nC87vwe3bW7XHM2Me8McFqwCtVex2Ueedu4JfsZh0o4KmMn9buXnqtvPvR0rsetZdxdhS5T/4ZVwaWqb7rygsyZeJd01+6ySq8/ffJ45DmuJ6OgnrX0dbIcJWvmOqv3c9rsdKf1/AYwf9XswwwArZZcWc81EThvf+jpXCsI1tPvOW+wiseTVYbQqM/2XagvXLN6cuW+Ro2Ne7A22qz2jHH1sc1zXLuXFlJ7JBe6TC+9cmLKv+eW7fsm/yhKG0eXWc2+yJ3W8xvA/EVgBhA0OFzW2oFdOm/zvVo7sEuDw7Pf7CYtRGUVLlSH3jyLkSQfM89rasVCJH1retW/YVVmR4pmcg/PeseSx7Z/wyoVusILlLinLwSybMniaeUotZYnP+munw5coQc3X5Y569/MBUwAIISSDACZ2nVSVT0fqVfPJqbVIleXcyTleU2z+RF/d7GQ2tVjtoXqjWPTZmprlFsfi8pVqktbztt8b+r+oeXJ88wSN3sBEwDIQmAGkKldJ1Vlhais3sbV0mqRS68/PbOPcq3XFFq9r7DItGTxoimdLgqLpLGTqbtP8+6LzgouZJLVdaJR4+4qFroynzdtue/qWeJqoRKJrFDc6MqQndLzG8D8RkkGgEztOqkq66P26y49d8pH/Wm1rVklJH1revXg5suUVVQQek1ZZR7dxYKuufgcVefIvGFZkrb90/Ta7CSX9LuXnpv/AXOKj198PJcvLahYOPUr4ZcSlweHyzVLOEIhN1Q6kVbCkVWzDADtwgwzgEyNfFzeiJl+1J6nhKTWawq1rqvePvTU0Ya7W+QJ183qoBFLW3Ww+pi8cHxMW7bv09BTR3XX7nCNt0l631uzZ3rj7Tfdc6rLSDKQM0sMoNOZN7kHaDOUSiUfGhpq9zCABS9tEZBioatjZwCzVuHr7S5OLiIyOFxW/7f2amz81Htfocu09bcvkqRpr7fQZVq2ZLGOVcamBOg/G9zX9CDbCl1muvnqiyZbuiUXU8nav1atszT1GKeZaz9LABYGM9vt7qVa+zHDDCDTXDupKncJSXX+i66n1WyPjfvkrGhyxvqOHz7T8HhjoVriZjvpPmXGt1Zdcp6wLE09xmmz9CwyAmAuIzADCJpLH5fnKSFJO3lt7KQHW5wlxSEvb5CsZe0bTtfvlM6dLIuo5yS/rkWm8RqBt1ryWORdhKWex80qi8n6g4BFRgDMBZz0B2DeyNOXNzQLnbc2+8hoRV0W7kmcx9o3nK7b//DXJ6+bJpbeXr60IJNqPsciTZysF993ym0pd52tHsXJx82aSc56LSwyAmAuIDADmDfydFwItT7r37BKhbSkmbLvtZecE9wn62F6u4s6NHCFDg1codv/8NenrWo4WhnTK2Mndcs1q3Xz1RcFn2PspGvpksU6NHCFbolWCoxf9+evXq0vXLN6ykIoyRPtpFNhuxHVxzjrD5K4jV0Si4wAmCsIzADmlbh9XNYqcTVnoWvk5Xjfz/ZdGGz3llUpUR0oG63tPTJamawZLo9WtMhM5dGKtu44qKGnjuoXJ0614Yg7X8St9q74tbNqPn5IdbcNaWJ57DTLlxZoHwdgzqKGGcCCEjqRce3ArindM2JdZjrpPu2kx8/2XajbH3q6roVFqme4a52o2F0sBGuNX1csTKkRjmury6OV1LElw/gDB0bqGPl0acE+q7Q7XjabgAxgLiIwA1hwsoJbVng96a6l3FgMAAAU2UlEQVSfDlyReltoFcBqJk0rQah1ouJNV16Q2fqtWOiSmTJPqMsK8uXRSmYLvnpVH7NjGeE+azsAzAWUZABAJFTfnCVrFcA0LtVdItK3pldbf+ci9UZjiE+ei0saRo/PLIjGHTnSmCZWF4zLJ7qLhcya7OpjM5NjCACdjhlmAIisO78ndTGSdef3ZN4nrcTj+Ksn9EJKkF2+tKC1A7umlIJIEyfjxbPE3cWCbrryginBOlTKkLZKX14uTWtjZ5Kuu/RcfbbvQkmn2sSl1WSnnbTXv2FV6gIlnNwHYC4jMANAJKumt1atbzLQDg6X9alv75+2T6HL9NIrp4J0ebSi/jv3SqYpddPJk/TyyAr5SaHV+lwTs9VZC9OknZQYP2baSXtzbbEbAMiDwAwAkdwrBWZIW/5Z0mRrt+qT99Lqkitj47ph2x5t3XFQ687v0QMHRoLBM8+Jeyfd1ZtRK11rSetQXXdWCObkPgDzDTXMABBptP42azZ22WsW133SW3m0om8+9PRkf+Z4xby4JVwsT5iPw/ZM+iBTkwwABGYAmDTTUBlrxiqCIXEbt6Raj5vslTyTPsiNHhMAmA8oyQCASKP1t6EWcWknw81EdShPe9z4RL7eqvHPpFSCmmQAIDADyCFeSW4hBKZG6m9rdYhIdsOYqeoZ5VYEWmqSASx0BGYAQdUnssW1tNL0nsILXVZ4lRScXS4WuvS+t/bqrt3lmoF63fk901rTEWgBYHaZZ61j2kalUsmHhobaPQwAUuaKcLW6K7TCXJn5Dq2qlyybiF9PPX2Vi4WuXLXI1ZLPFbedqy7hAID5zsx2u3up1n7MMAMIarTV2myZSzPfWcfKpCl/dMQzxVnt6dJUxsb1ib/bO3n/PKofP+7R3MnHEADaiS4ZAII6ta1YWgu3tC4SnWAmx/A1i0+9PVvWGtaRcffUlnNZstrfSZ17DAGgnQjMAII6ta1Yq2e+B4fLWjuwSys336s3bLlPKzffq7UDu3KF1HqOYTz7m1zkJE/lXD1Bt9YxavenBwDQaSjJABDUqW3FQi3cmq3REoZ6jmFo9reWvLXPWccueTsA4BQCM4CaOrELQ1b/4XXn9zT9ufKUMNQ6PnmPYSOzu121ajcioZ7QnfDpAQB0GkoyAMxJfWt69b639ioZEV3SXbvLuWt582plCUPW7O7ypYXJVfqyjOfsepRc9U86FbTzrv4HAAsNM8wA5qwHDoyoOiLmnfGtRytLGLIWP7nxPRdMvqZQq7+8OvFTAwDoVMwwA5izWnXiX9pJe7FmlzAkZ39N6bO+WWUns1GOAgBghhnAHNaqE/+SJ+21YqGPWrO/DxwYqWs7AKAxBGYAc1ZW+cJsnLTWSSUMnbqYDADMV5RkAJiz8pQvzEedupgMAMxXzDADmNM6aea3VVo5sw4AIDADmCcGh8sdt7jKbOnUxWQAYL4iMAOY86pX4su7At9Mn6sTgupCnFkHgHahhhnAnJe2El/cj7mZ4mBeHq3IdSqYN3uhFABAZ6kZmM3sNjN73sweS2zbamYHzOxRM7vbzLoz7nvIzPaZ2R4zG2rmwAEg1qquEa0K5gCAzpJnhvlrki6v2rZT0pvd/dck/bOkLYH7r3P31e5emtkQASCsVV0jaOcGAAtTzcDs7t+XdLRq2/3ufiK6+pCks2dhbACQS9pKfLPRNSIrgHcvLTT1eQAAnaUZNcwfkvSdjNtc0v1mttvMNoUexMw2mdmQmQ2NjLBaFYD8WtWPuX/DKhW6bNr2l145QR0zAMxj5u61dzJbKekf3P3NVds/Kakk6SpPeSAz63X3spn9iibKOP44mrEOKpVKPjREyTOAzrP6U/drtDI2bXtvd1EPbr6sDSMCAMyUme3OUzY84xlmM/ugpHdLui4tLEuSu5ej789LulvSxTN9PgDoBMdSwrJEHTMAzGcz6sNsZpdL+hNJ/9Hdj2fss0zSInf/eXR5vaRPz3ikANABVnQXVU4Jx42eYNgp/Z0BANPlaSt3h6R/lLTKzA6b2R9I+qKk0yTtjFrGfSnad4WZ3Rfd9UxJPzCzvZL+SdK97v7dWXkVANAis3GCIf2dAaCz5aphbjVqmAF0smbPBq8d2JU6a01dNADMrrw1zCyNDQB1avay1PR3BoDOxtLYANBmrVp4BQAwMwRmAGizVi28AgCYGUoyAKDN4vIOumQAQGciMANAB2h2XTQAoHkoyQAAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACFjc7gEAAAAkDQ6XtXXHQR0ZrWhFd1H9G1apb01vu4eFBYzADAAAOsbgcFlbtu9TZWxcklQerWjL9n2SRGhG21CSAQAAOsbWHQcnw3KsMjaurTsOtmlEAIEZAAB0kCOjlbq2A61AYAYAAB1jRXexru1AKxCYAQBAx+jfsErFQteUbcVCl/o3rGrTiABO+gMAAB0kPrGPLhnoJARmAADQUfrW9BKQ0VFylWSY2W1m9ryZPZbYdrqZ7TSzJ6LvyzPue320zxNmdn2zBg4AAAC0Qt4a5q9Jurxq22ZJ33P3N0r6XnR9CjM7XdKNki6RdLGkG7OCNQAAANCJcgVmd/++pKNVm98r6evR5a9L6ku56wZJO939qLu/IGmnpgdvAAAAoGM10iXjTHd/Nrr8L5LOTNmnV9IzieuHo23TmNkmMxsys6GRkZEGhgUAAAA0T1Payrm7S/IGH+NWdy+5e6mnp6cZwwIAAAAa1khgfs7MzpKk6PvzKfuUJZ2TuH52tA0AAACYExoJzPdIirteXC/p71P22SFpvZktj072Wx9tAwAAAOaEvG3l7pD0j5JWmdlhM/sDSQOS3mlmT0h6R3RdZlYys69IkrsflfQZSQ9HX5+OtgEAAABzgk2UH3eWUqnkQ0ND7R4GAAAA5jEz2+3upVr7NeWkPwAAAGC+IjADAAAAAQRmAAAAIIDADAAAAAQsbvcAAGCuGxwua+uOgzoyWtGK7qL6N6xS35rURU0BAHMQgRkAGjA4XNaW7ftUGRuXJJVHK9qyfZ8kEZoBYJ6gJAMAGrB1x8HJsByrjI1r646DbRoRAKDZCMwA0IAjo5W6tgMA5h4CMwA0YEV3sa7tAIC5h8AMAA3o37BKxULXlG3FQpf6N6xq04gAAM3GSX8A0ID4xD66ZADA/EVgBoAG9a3pJSADwDxGSQYAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAiYcWA2s1Vmtifx9aKZ3VC1z9vN7Fhinz9vfMgAAABA6yye6R3d/aCk1ZJkZl2SypLuTtn1/7r7u2f6PAAAAEA7Nask4zcl/djdn2rS4wEAAAAdoVmB+f2S7si47dfNbK+ZfcfMLsh6ADPbZGZDZjY0MjLSpGEBAAAAjWk4MJvZEklXSroz5eZHJL3e3S+S9L8kDWY9jrvf6u4ldy/19PQ0OiwAAACgKZoxw/wuSY+4+3PVN7j7i+7+UnT5PkkFMzujCc8JAAAAtEQzAvO1yijHMLNfNTOLLl8cPd+/NuE5AQAAgJaYcZcMSTKzZZLeKemPEts+LEnu/iVJvy3pv5jZCUkVSe93d2/kOQEAAIBWaigwu/vLkn65atuXEpe/KOmLjTwHAAAA0E6s9AcAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAggMAMAAAABBCYAQAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQACBGQAAAAggMAMAAAABBGYAAAAgoOHAbGaHzGyfme0xs6GU283M/tLMnjSzR83sLY0+JwAAANAqi5v0OOvc/WcZt71L0hujr0sk/XX0HQAAAOh4rSjJeK+kb/iEhyR1m9lZLXheAAAAoGHNCMwu6X4z221mm1Ju75X0TOL64WjbFGa2ycyGzGxoZGSkCcMCAAAAGteMwPw2d3+LJkovPmJmvzGTB3H3W9295O6lnp6eJgwLAAAAaFzDgdndy9H35yXdLeniql3Kks5JXD872gYAAAB0vIYCs5ktM7PT4suS1kt6rGq3eyT9XtQt41JJx9z92UaeFwAAAGiVRrtknCnpbjOLH+tv3P27ZvZhSXL3L0m6T9JGSU9KOi7p9xt8TgAAAKBlGgrM7v4TSRelbP9S4rJL+kgjzwMAAAC0Cyv9AQAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIaHSlPwAAALTA4HBZW3cc1JHRilZ0F9W/YZX61vS2e1gLAoEZAACgww0Ol7Vl+z5VxsYlSeXRirZs3ydJhOYWoCQDAACgw23dcXAyLMcqY+PauuNgm0a0sBCYAQAAOtyR0Upd29FcBGYAAIAOt6K7WNd2NBeBGQAAoMP1b1ilYqFryrZioUv9G1a1aUQLCyf9AQAAdLj4xD66ZLQHgRkAAGAO6FvTS0BuE0oyAAAAgAACMwAAABBAYAYAAAACCMwAAABAAIEZAAAACCAwAwAAAAEEZgAAACCAwAwAAAAEEJgBAACAAAIzAAAAEEBgBgAAAAIIzAAAAEAAgRkAAAAIIDADAAAAAQRmAAAAIMDcvd1jmMbMRiQ91e5xpDhD0s/aPYh5gOPYHBzH5uA4NgfHsXEcw+bgODbHQjmOr3f3nlo7dWRg7lRmNuTupXaPY67jODYHx7E5OI7NwXFsHMewOTiOzcFxnIqSDAAAACCAwAwAAAAEEJjrc2u7BzBPcBybg+PYHBzH5uA4No5j2Bwcx+bgOCZQwwwAAAAEMMMMAAAABBCYAQAAgAACcwozO2Rm+8xsj5kNpdxuZvaXZvakmT1qZm9pxzg7mZmtio5f/PWimd1Qtc/bzexYYp8/b9d4O4mZ3WZmz5vZY4ltp5vZTjN7Ivq+POO+10f7PGFm17du1J0n4zhuNbMD0f/bu82sO+O+wfeAhSTjON5kZuXE/92NGfe93MwORu+Vm1s36s6ScQy3JY7fITPbk3FffhYjZnaOmT1gZj8ys/1m9rFoO++PdQgcR94fA6hhTmFmhySV3D21YXf0y+GPJW2UdImk/+nul7RuhHOLmXVJKku6xN2fSmx/u6T/6u7vbtfYOpGZ/YaklyR9w93fHG37H5KOuvtAFDyWu/t/q7rf6ZKGJJUkuaTdkt7q7i+09AV0iIzjuF7SLnc/YWb/XZKqj2O03yEF3gMWkozjeJOkl9z9LwL365L0z5LeKemwpIclXevuP5r1QXeYtGNYdfvNko65+6dTbjskfhYlSWZ2lqSz3P0RMztNE+9xfZI+KN4fcwscx7PF+2MmZphn5r2aeONzd39IUnf0A4h0vynpx8mwjGzu/n1JR6s2v1fS16PLX9fEm1u1DZJ2uvvR6JfATkmXz9pAO1zacXT3+939RHT1IU38gkBAxs9jHhdLetLdf+Lur0r6W038HC84oWNoZibpakl3tHRQc5C7P+vuj0SXfy7pcUm94v2xLlnHkffHMAJzOpd0v5ntNrNNKbf3Snomcf1wtA3p3q/sXwa/bmZ7zew7ZnZBKwc1x5zp7s9Gl/9F0pkp+/BzWZ8PSfpOxm213gMgfTT66Pa2jI/A+XnM5z9Ies7dn8i4nZ/FFGa2UtIaST8U748zVnUck3h/rLK43QPoUG9z97KZ/YqknWZ2IJohQJ3MbImkKyVtSbn5EU2s4f5SVOYyKOmNrRzfXOTubmbUUjXAzD4p6YSk2zN24T0g7K8lfUYTvzg/I+lmTfyCRf2uVXh2mZ/FKmb2Wkl3SbrB3V+cmKSfwPtjftXHMbGd98cUzDCncPdy9P15SXdr4qPFpLKkcxLXz462Ybp3SXrE3Z+rvsHdX3T3l6LL90kqmNkZrR7gHPFcXPYTfX8+ZR9+LnMwsw9Kerek6zzjJI4c7wELmrs/5+7j7n5S0peVfnz4eazBzBZLukrStqx9+FmcyswKmgh5t7v79mgz7491yjiOvD8GEJirmNmyqAheZrZM0npJj1Xtdo+k37MJl2riZI1nhTSZsydm9qtR/Z7M7GJN/Dz+awvHNpfcIyk+q/t6SX+fss8OSevNbHn0Efn6aBsiZna5pD+RdKW7H8/YJ897wIJWdc7Gbyn9+Dws6Y1mdl70SdP7NfFzjFPeIemAux9Ou5Gfxami3xdflfS4u38+cRPvj3XIOo68P9bg7nwlviT9W0l7o6/9kj4Zbf+wpA9Hl03SX0n6saR9mjhbtO1j77QvScs0EYBfl9iWPI4fjY7xXk2cYPDv2z3mTvjSxB8Yz0oa00Sd3R9I+mVJ35P0hKT/I+n0aN+SpK8k7vshSU9GX7/f7tfSgcfxSU3UMe6Jvr4U7btC0n3R5dT3gIX6lXEc/3f03veoJsLKWdXHMbq+UROdMn68kI9j2jGMtn8tfj9M7MvPYvZxfJsmyoAeTfwf3sj7Y9OOI++PgS/aygEAAAABlGQAAAAAAQRmAAAAIIDADAAAAAQQmAEAAIAAAjMAAAAQQGAGAAAAAgjMAAAAQMD/B+9qJR3RJP9CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f60539104a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter(X[:,0], X[:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian distribution:\n",
    "### Estimating parameters for a Gaussian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimateGaussian(X):\n",
    "    mu = X.mean(axis=0)\n",
    "    sigma = X.var(axis=0)\n",
    "    \n",
    "    return mu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([14.11222578, 14.99771051]), array([1.83263141, 1.70974533]))"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu, sigma = estimateGaussian(X)\n",
    "mu, sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selecting the threshold, ε:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((307, 2), (307, 1))"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xval = data['Xval']\n",
    "yval = data['yval']\n",
    "\n",
    "Xval.shape, yval.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.183842  , 0.20221694, 0.21746136, 0.19778763, 0.20858956,\n",
       "       0.21652359, 0.16991291, 0.15123542, 0.1163989 , 0.1594734 ])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "dist = stats.norm(mu[0], sigma[0])\n",
    "dist.pdf(X[:,0])[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# probability distribution of the dataset\n",
    "p = np.zeros(X.shape)\n",
    "p[:,0] = stats.norm(mu[0], sigma[0]).pdf(X[:,0])\n",
    "p[:,1] = stats.norm(mu[1], sigma[1]).pdf(X[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# probability distribution of the validation dataset\n",
    "pval = np.zeros(Xval.shape)\n",
    "pval[:,0] = stats.norm(mu[0], sigma[0]).pdf(Xval[:,0])\n",
    "pval[:,1] = stats.norm(mu[1], sigma[1]).pdf(Xval[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def selectThreshold(pval, yval):\n",
    "    bestEpsilon = 0;\n",
    "    bestF1 = 0\n",
    "    F1 = 0;\n",
    "    \n",
    "    step = (pval.max() - pval.min()) / 1000\n",
    "    \n",
    "    for epsilon in np.arange(pval.min(), pval.max(), step):\n",
    "        preds = pval < epsilon\n",
    "        \n",
    "        tp = np.sum(np.logical_and(preds == 1, yval == 1)).astype(float).astype(float).astype(float)\n",
    "        fp = np.sum(np.logical_and(preds == 1, yval == 0)).astype(float).astype(float)\n",
    "        fn = np.sum(np.logical_and(preds == 0, yval == 1)).astype(float)\n",
    "        \n",
    "        P = tp / (tp + fp)\n",
    "        R = tp / (tp + fn)\n",
    "        F1 = (2 * P * R) / (P + R)\n",
    "        \n",
    "        if F1 > bestF1:\n",
    "            bestF1 = F1\n",
    "            bestEpsilon = epsilon\n",
    "    \n",
    "    return bestEpsilon, bestF1        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.009566706005956842, 0.7142857142857143)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epsilon, F1 = selectThreshold(pval, yval)\n",
    "epsilon, F1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([300, 301, 301, 303, 303, 304, 306, 306]),\n",
       " array([1, 0, 1, 0, 1, 0, 0, 1]))"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# indexes of the outliers\n",
    "outliers = np.where(p < epsilon)\n",
    "outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f6052e17320>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAHVCAYAAAAU6/ZZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X+Q3Gd9J/jPM6OWGDkpjwg6sCdszFEpm/W6bIEKcuvcVkx2MYsDUWyMl8rtsrepYrdqcxUo0J0cUrGdDWXdCYe9q7vKFXuhwh4cZ2w5iomzNhx2VWq5M7uSZWEc20ey/hEGG5SzR/zQxB5Jz/0x3eOenu+3f8x09/Toeb2qpjT69q+nv+ppvefpz/N5Us45AACgRFObPQAAANgswjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWNvG+WCve93r8iWXXDLOhwQAoEDHjh3765zz7l7XG2sYvuSSS+Lo0aPjfEgAAAqUUnq2n+spkwAAoFjCMAAAxRKGAQAoljAMAECxhGEAAIrVMwynlN6YUnoopfTnKaXHU0q/0Tx+a0ppPqX0aPPrPaMfLgAADE8/rdXORMTHcs6PpJR+MiKOpZS+2rzs0znnT41ueAAAMDo9w3DO+fmIeL75/Q9TSk9ExNyoBwYAAKM2UM1wSumSiNgTEd9oHvr1lNI3U0qfTSntqrnNh1NKR1NKR0+ePLmhwQIAwDD1HYZTSj8REYcj4iM55x9ExO9HxJsj4qpYnjm+o+p2OefP5Jz35pz37t7dc0c8AAAYm77CcEqpEctB+As553siInLO38s5n805n4uIfxMRbx/dMAEAYPj66SaRIuIPIuKJnPPvtR2/qO1qvxIR3xr+8AAAYHT66SZxdUT844h4LKX0aPPYb0bEB1NKV0VEjohnIuKfj2SEAAAwIv10k/j3EZEqLvrT4Q8HAADGxw50AAAUSxgGAKBYwjAAAMUShgEAKFY/3STOC/ff+Eg0Dp+KqRxxLkUs3XBhvPuut272sAAA2ERFzAzff+MjsePuUzGdU6RIMZ1T7Lj7VNx/4yObPTQAADZREWG4cfhUpI7ucClSNA6f2qQRAQAwCYoIw1N5sOMAAJShiDB8rmrLkC7HAQAoQxFheOmGCyPH6mngHDmWbrhwk0YEAMAkKCIMv/uut8bL778wzqYcOXKcTTlefr9uEgAApSumtZrgCwBApyJmhgEAoIowDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACK1TMMp5TemFJ6KKX05ymlx1NKv9E8/tqU0ldTSt9u/rlr9MMFAIDh6Wdm+ExEfCzn/Lcj4uci4l+mlP52RByIiK/lnH82Ir7W/DsAAGwZPcNwzvn5nPMjze9/GBFPRMRcRPxyRHyuebXPRcS+UQ0SAABGYaCa4ZTSJRGxJyK+ERGvzzk/37zohYh4fc1tPpxSOppSOnry5MkNDBUAAIar7zCcUvqJiDgcER/JOf+g/bKcc46IXHW7nPNncs57c857d+/evaHBAgDAMPUVhlNKjVgOwl/IOd/TPPy9lNJFzcsviojvj2aIAAAwGv10k0gR8QcR8UTO+ffaLro3Ij7U/P5DEfHHwx8eAACMzrY+rnN1RPzjiHgspfRo89hvRsTBiPhSSunXIuLZiPjAaIYIAACj0TMM55z/fUSkmot/cbjDAQCA8bEDHQAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYm3b7AEAAIM7cnw+Dj3wVHx3YTEunp2J/ddeGvv2zG32sGDLEYYBYIs5cnw+br7nsVhcOhsREfMLi3HzPY9FRAjEMCBlEgCwxRx64KmVINyyuHQ2Dj3w1CaNCLYuM8MAW4SPxWn57sLiQMeBemaGAbaA1sfi8wuLkePVj8WPHJ/f7KGxCS6enRnoOFBPGAbYAnwsTrv9114aM43pVcdmGtOx/9pLN2lEsHUpkwDYAnwsTrtWeYyyGdg4YRhgC7h4dibmK4Kvj8XLtW/PnPALQ6BMAmAL8LE4wGiYGQbYAnwsDjAawjDAFuFjcYDhUyYBAECxhGEAAIolDAMAUCxhGACAYllABzDBjhyf10ECYISEYYAJdeT4fNx8z2Mr2zDPLyzGzfc8FhEx9kAslAPnK2USABPq0ANPrQThlsWls3HogafGOo5WKJ9fWIwcr4byI8fnxzoOgFEQhgEm1Hcrtl/udnxUJiWUA4yCMAwwoS6enRno+KhMSigHGAVhGGBC7b/20phpTK86NtOYjv3XXjrWcUxKKAcYBWEYoA9Hjs/H1QcfjDcduC+uPvjgWOpl9+2Zi9uvvyLmZmciRcTc7Ezcfv0VY1+4NimhHGAUdJMA6GEzuzrs2zO36V0bWo+vmwRwPhKGAXrotoBsVIFw0lqZrSeUT9pzAKgiDAPntWEEsnEvIPutI4/FFx5+LnLz771moicxdE5Sj2SAbtQMA+etYfXHHecCsiPH51cF4Za6Vmbj7AE8SN20dmzAViEMA+etYQWycS4gO/TAU2uCcEvVTPS4QuegoVs7NmCrUCYBnLeGFcjGuYCs29imUoo3Hbhv1eOPK3TWhe5b73288rxcPDsT8xVj0I4NmDTCMHDeGmYgG1dXh7oxR0Sczctzxu31t6MInVU1yHXhemFxKRYWl9aMa/+1l66qGY7Qjg2YTMokgPPWJPfHrau/rRpzlVYpxLCfY105xOzORl+3b++y0d4jedfORuzYNhUfufPRePPNfxqXjLFfM0A3ZoaB89ak9Mc9cnw+bvvy4/HS6eUZ1JnGVJw5l2Pp7NqZ3qox180Uf3dhceDn2KvzRF05xI5tUzHTmF5zWd242uWIWDi9tFILXTXDrcMEsFlSznVLNYZv7969+ejRo2N7PIDNduT4fOy/+8RK8O1mbnYmvn7gnWuOX33wwcpAXHf99se+9d7HV8oYLtg+Ha+cORdL514dy0xjOm6//oqIWA7CdcE7RcSnb7pqVZA+/cqZlYDfOa6qMok6vZ5H53Pa7F9ugK0hpXQs57y31/XMDAMMaJBAduiBp/oKwhH1i97WU3975Ph87L/rxKrg++NX1gbT1iK4l8+c6xpcL56dWVM33dlLuH1cVTPMdVrPu9d51bsYGAVhGGAAgwayQbo61C16W0+5x6EHnloVhLtpzRzXqQve3cb10Tsf7euxI5afd9V53X/3ibj13sfj1OLSykz0uHcCBPp33yefiFc+9ULMLkQszEZs//gb4rpPvGWzh9WTMAzQReds5aCBrFvNb7teM71V3Sy6zaQOq7Xa7Ewjbn3f5bVhs31crfF89M5HYyqlldrgbhpTqXYmeelsXtWpoo7exbD57vvkE7HtthfigqUUERG7FiJevu2FuC9i4gOxMAxQo2q2sk63EoeqmuGpiLhwZyNeOr0U0ymtBOqjz74YDz15sucMcK8Z6kFC+GsaU5W1vxERF+zY1tesa+d4+gnCEbFcjBwbC7R6F8Pme+VTrwbhlh1LKU5/6oWICQ/DWqsB1Bik7nV2Z6OyVdq+PXNx6P1Xxq621mSzM434vZuuilvee3nMNKZXdVf4/MPPdd3lrdWS7SN3Ptp157n9114ajanV/zFVuf36K+KW915ee3k/gfrI8fn42JdOVJ6r6ZQixfIiudmZte3Zls7mOPTAU+sOtI3pNBGt8qB0swuDHZ8kZoYBavQ7W9mYTvGjv3m1s0LnLG3dhh1XH3ywZ9huL8GoWrBWN+bW47V3k+g011wUFxHx0TsfrdwGejp1D9StMdXNBJ/LOZ4+eF1ERLzpwH21Y/70TVf13X2i3QXb+5u5BkZrYXa5NKLq+KQzMwxQo262cnamsbKZxNzsTFywfduaxWrts7R1+g3b8wuL8aYD99XOvrbLESsz0/v2zMWjt7wr/vVNV/XcmKOuqKFXucOt9z7eswtFS+3GHWk5jO/YNhW7djZWNunoZ2b7VI/Ff8B4bP/4G+Llxur3i5cbObZ//A2bNKL+mRkGqFHX0qxzQVm3Gc9u+q3rjVgOq/3W4VbNTEdUd31oLXqrM9elfOHI8fmunSjWBO6a4beOLywuxUxjOj5901WrxvbdhcXaBXnqhWEyXPeJt8R9EXFaNwmA80e/Lc3qQm1nUOvs/nDNZbvj8LH5gUsD+tHZ4aKuG0W30oT2MFvVuaJbiJ5OKW6//opVj9mrhVvnuDs7VQzaaxkYr+s+8ZaJXyxXRRgG6KKu3rddP5tiVHV/OHxsPm5421z8yYnnV4Lijm1T8cqZc7VlC4PoNTPdbYHgXMfscVXnim4h/o4PXLnmvE332W6tatyTsrU2cP4RhoFirHcr31636yeoVQXPxaWz8Scnno+Xz5xbOdb+/Ua1NrMYtBdxili1PXLd2OvCbV2pb79lHt02HxF+gWEThoEirHcr335vV7VV8dUHH1wJoXW1wf2UDvQjxepFcDON6bjkp2ZWdYmYX1iMj975aBx99sX43X1X9F3eUReaz+YcM43pNUH5XI7Kc9TPzLDSB2DcdJMAilA3u9mr48N6btcK0O39gkctR6zqcHHD2+bi//7LF9eUW+SI+PzDz8UlB+6LH798JhrTq6dxq8Jo3Uxt63GqVJ2jfmaGd2x79b+l1i8Unb2bAYbJzDBw3jtyfL42kLbaltWVTdTNis4vLK6a+W2/bb+bdfTa/W0Qu3Y2VkobWptg9IqeC4tL0ZhKsWtnIxZOL9Weg7qa6NYCwDqd53yuZia6fVZ7YXEpbr7nsTj67IurFhf2O5MPMKiU+90ycwj27t2bjx49OrbHA85Pg9T+9rNRRctMY3pNB4SrDz7YM8BFLG+8ccH2bT3LHqZTinM5x8WzM3HJT81Uzt6uR2MqxaEbr4yIGHjzitmZRlywY9uqLhedW0JHrK2JPvTAUz1nvZ9pbrgRUf1v0XkeW6bScrlFp7nZmVX1zAB1UkrHcs57e15PGAa2kroWW7dff0VErC+wtesMW4MEuH6kiHj64HXxW0cei88//NxAt52dacSt77u8dre4Vk/gYZdlVP2SELHcX7nXeWgPwxFrf5EZdKyt8wfQS79hWJkEsKXU1fDeeu/j8fKZcwO1/6rSWRZR1SliI2GztQvbF7/xVwPf9tTiUuzbMxcfvfPRystHVZu8uHQ2bvvy43H02Rfji9/4qzibc0ynFDONqTi9VN/9IkVUlpL0M/NexyYbwLBZQAdsKXU1vAuLS7Xtv6rUHe8Vtk6/cqaPUdZ76fRSHDk+33ebsXatsW1GIHzp9FJ8/uHnVsZ9Nuc4vXSuto1aRKwsHmz9efM9j61ZBLf/2kvXbBXdjU4TwLAJw8CWMmgQbLX/ajfTmI4PvuONlcc7w1ZnZ4hhLHb7SM3MbsvsTKNybNdctjv2/M5XNjQDfMH26Zidaaz79p0GifRVHSb27ZmL26+/YlUnjG4sngOGTRgGtpSqmcSZxnTs2lkd8OZmZ9aErduvvyJ+d98Vlcc7w1a/nSGGZSoibn3f5avGtmtnI1Lk+PzDz204jJ9+5WxcsGN4FXI59w6w7ep2l/v6gXfG0wevszgOGDs1w8CWUrfbW8TaLgqtmd66ncv62dGs15bGw3ZhM9S3xvbqAr7h7Ex38exM1+eUYrkzxitn+5vznU4p9l97aey/60QsVbV/qHh8gEkiDANbTrcQu57tlrvpd8FcaraYaH/cPb/zlYFncl86vbSqn+6wZ6Z7ddjIEbFz+7aYbluM2M0H3/HG2LdnLm778uM9n2u/u8vt2tmovK+62f86691+GyiLMgngvFH1cftGdzDrd4FX1Xq4W957eUx3W2FWo722dpgz042p5XN0zWW7u17v1OLSqjKNusWGOxtT8bv7llvaLfQIwnVlKFVuee/la3bGa0ynuOW9l/e8bUvVLoBVC/gAzAwDA5m02ba68XT2B+53B7Oq+7v9+iv66lfcHrpapiJiPfO6rRB84Uyj50Ye/TrTrLR46MmTXa938ezMqjKNugV/7W3V6mbQd+1sxPHffldf42s/9xfONCKl6LozXjfdttE2Owy0E4aBvq03YG7GeOrC0Me+dCIiVo+3FcLmFxZXbajRur/br78ivn7gnXHk+Hzsv/tELPWop209zk++ZltfdbRVWrW1NZOyG7rPbrPN7aUMrfNbJzWvs2/PXO2Wzf3O5nb+Wy4sLsVMYzo+fdNV63pt1T3HcdeAA5OvZ5lESumzKaXvp5S+1Xbs1pTSfErp0ebXe0Y7TGASdJttm7Tx1IWeszmv+ri8/eP0iLWtwtqf3749c3Ho/Vf2Vbt6NucNzei2Ammv8oN+NaZSnH7lTLzpwH0x1aX3cnspQ6965dy8TkR1i7R+yyLqHmsjr626hXoW8AGd+pkZ/sOI+J8j4t92HP90zvlTQx8RMLEmbbat23i6LXxr/7i8nwVq7Y/T2c1iKqV1baDRy2/e88346J2PDuX+U0Sci1d7JFfdX9WWy/38u3aem/V+QjDs11bdTLVNO4BOPWeGc85/FhEvjmEswISbtNm2buPptfCtFbL6CVvtj9O5MGsUQThiuR53WPefI+JsRblGSlE7i3vk+HztDHK7qZSGsiht2K+tjc5UA+XYSM3wr6eU/klEHI2Ij+WcXxrSmIAJNWmzbd3G0wo9H/vSicpAOdssdejVOq29vODi2Zn48ctnxroJxyjlHPHMwevWHG8F/n6CeKvsJGJjdeOjeG1tZKYaKMd6W6v9fkS8OSKuiojnI+KOuiumlD6cUjqaUjp68mT3FczAZJu02bZe49m3Zy7u+MCVa9p0RUT86G/OLC+I6zGDvHQux0unl1Y6RXSrA+6nBdtWUFc60ppJ7tRvbe+R4/O1re4m7bUFlCPlPn7zTyldEhF/knP+O4Nc1mnv3r356NGjAw8SoNMgLd6uuu0rlSF2bnZmpUtEP5tGdDPXHEP7mPrZrGMzzc404tFb1rY9e9OB+9YsJOwlRcTTFbPMLZ3dIiKq65QBhiWldCznvLfX9dZVJpFSuijn/Hzzr78SEd/qdn2AYakKrr1avJ2qmc1t1Qu3FtKtNwy3Sik+euejcfHszEo7sEsO3Leu+xuHxlSKW9+33Pas/ReL2S6dMtrbznXqVdur7y8wqXqG4ZTSFyPiFyLidSml70TELRHxCymlq2L5ffGZiPjnIxwjQERUzy62dAtWdbO07QFukK4FF2yfjtmd21c2h/jB3yytBOn5hcXYf9eJvu9rvboF017mZus3J6n7haDb4/VT2ztot4hJ29wFOH/1DMM55w9WHP6DEYwFKFg/4adXG7T5hcW4+uCDa27bz+KsQcoaFl85G9u3nYmI5VnnzpC4dC7Hzfd8M2aHuHtcy3RKsWNbWrX7W+flZ3OuDK8pIn715/7WyhbKEb3Paft91qkrdWj/N61rEVc1ozxpm7sA57f1LqAD2JD2xVRX3faV2H/3iZV2Za3w09myq5/Z26rb9lqcdeT4fCycfqXvsbd69uaony1dXDoXv3TlRdGYGt4Wco3pFHd84MraIBwRcccHroyZxnTluHJEHD42v+rc9HNOz+UcczVlELt2NuLQA0+tWRTXTwu6uhnlSdvcBTi/CcPA2HUGpYXFpTVbHFeFn357zlbddt+eufj6gXfG0wevi68feOeqIHzzPY/Fj18Zfru0h548GYduvHIlhE9vcG/lpbM5PnLno12vc9uXH+8609t5bvo5p7M7G5VdNxrTKX70N2cqf4mpm3GeTqlnt4hJ29wFOL8Jw1Cwbq2uRqmfj+Yj1oafXm3Q2rVKJno9p37Hsh7fXViMfXvmYv+1l8bFszMr5Quj1M8iwPbzuv/aSytbz7XLuXp2/YLt22LpXPUvMXXB9VzOa34h6TRpm7sA57eNbLoBbGGbWZfZ7wxfZ/jp3Aq5tQlGXV1uP89plLONF8/OxG8deSy+8PBzK2ULo9mvbjBrQmWPQbW6cXRuYvGmmm4Z3VrL9RNoJ21zF+D8ZmYYCrWZdZn9BKK68NNZ7nDr+y7vOlvc6zl1G0trvrS97LfR57vmVERcc9nuVUG43c7GVOzq0sZsVDrP66EHnlozu9tp0Jnauu2w+w20NuAAxsnMMBRqM+syq2b+GlMpfuI122Lh9FJtN4luHScOPfBUbTeIbs+paiwRyxtS/NKVF8XhY/OrLuuydm2VcxFxz7Hv1E66nl46F9e/7afX3P8wzM404oId21b6Bue8XJc9ndKaXw56ddBoTKXaANvPdtjrbY9mK2VgXIRhKNRGPsbeqPUEpV5lHfv2zMXVBx/s+px6hen240effTE+//BzG3qe3bo+RMSG77/KTGM6bn3f5Ss9hG+99/GVMpJWR4f5hcXYf/eJ/mo2upQTt85d+2O8pm3qXKAFtoK+tmMeFtsxw+TYatvj1gXd1pbKEbGmPjfi1ecUEWtno6dTXLB9W5xaXD0b/VtHHhtJUB216bTcem3Yz6H9HHfaaq8joBwj3Y4Z2Po2+jH2uPUq6zhyfD4OH5tfFYRTRNzwtldnjTvLEZbO5pUZzfaZ5i9+46+GMuaUIl6zbXpk3So6nct5ZUb4C0MM8+3nuPP1YptlYKsThqFgW+lj7F5lHVWhLMdyr9+I/mqhWyGu225rg/jVd/yt2Pszr43bvvx4Xy3PWq5+82vjPzz9Us+FbZ3az8UwP/O7eHamtkylLujrCQxsFbpJAFtCr+4EvWaO+62F/u7CYs/NMXr1Ck4p4r9q2/Z45/bleYfW/fa6/0eeOxU3vf2Nld0mZhrTcfWbX7vmePtCt2EG0dY5rpsBrnsuegIDW4UwDGwJvdpt9Wr/tf/aS/vaGvni2Zn44DveWHv5VKpedzY3OxPPHLwunjl4XTx9+3Xxu/uuWLXTXsTyAraZxnTX+49YDpkPPXkyjv/2u+Jf33TVylbIrW4Qj3/3hzHd+Vza/rqRIJoi4oLt02vOcV3Abj2ndnoCA1uJMglgy+hW1tHXRg09snB7W7A7/8NzlW3U6ioXqsJi3Wxqq3Sjm+8uLK7U6M4vLEaKV7tBVG0ysnQ2r9TpXnPZ7nUvnsux/Bw/fdNVq8717M5GZanHrp2NuOW9l2+Z2nOATsIwcF7otSDw0ANPxdLZtUl2OqU4l/Oa65/ps59wS9VsbLfSjdmZRu3OeRERF840VoX7fmqAW1tQ//jlM32NuU7VAri6MurWVs3CL7BVCcPAeaNbKKsLpudyjqcPXrfmeN2CvSoporIsoNuiv/3XXhr77zpRuUhupjEdKcW6ulD0O+ZeOs/XqZrgXnccYKtQMwwUYdAthasW7NXJEZUhvNuiv3175uLQjVeuqgeOeLVOd2GA7hP9SrG8sK/1mN2qRjrPy6DnD2CrMDMMFKGujvaay3ZXXr+q7OL0K2cq62bnmq3HOks0IpZ3ZGvN8M7ONFZ2h2s9Rt1MdrftpVt2NqZ67nLXkiLiV9s6XNRtYhJRvQCur5psgC1IGAaKULdordtits6wWrfD3TWX7V7Tg3f/XSciUqyqU355gELkfhbB7bpgR+yK6tKIXTsbsXP7ttpFbd3ar1XtHrfVNmkB6JcwDBShVx/iXrrtcPfQkyfX7m5XUQu8uHQ2PvalE3H02RfjoSdPdg2V/Xac+PRNV1XO2N7y3su7BtW6eua52Zna21koB5yP1AwDRdhozWu3He4G2eTibM7x+Yefi/mFxcjx6k5uR47Pr7peP/d5cTO4duu/XKfXJiYApTAzDBRhozWv3WaWB+k8UaWqlVmv+2wf+3pmbJU9ACwThqFwVQu/zsdAtNHw16tNWmfQHlRn2K66zxTLs9FzQ/p3UvYAIAxD0VrbBbcv/Lr5nsciorpV2Fa3kfBXN7N8zWW7K0soBnVxRUeKVj3y+f6LCsBmEoahYHXbBXd+ZE/1zPI1l+2Ow8fmuwbhxnSKyNUL6tpVdaQ4fGy+r/rfOqXM+gNshDAMBdtoh4VRmsQg1zmzfPXBByuDcOcWzxHd+waniMo2auv5xaR13uYXFlfKKiLO/1l/gPUShqFg3epgN9NWKd8YZIvnfXvm1jyvlm5zxvMLi3Hk+Hxfz7vz/jvv16w/wFpaq0HBJrW9VrfyjVE4cnw+9vzOV+KSA/fFJQfui6tu+8qaVmdVBm3Xtm/PXNzwtrmVrZf7VdV6rUo/tcuTMOsPMEmEYSjYenvUjto4yzeOHJ+P/XefWLXN8sLiUuy/60TPADroLxOtjTvO5u71w536/UWg397EALxKmQQUbhLba42zfOPQA0+t2jK5Zelc7llSMGi7to10neg36PbbmxiAZcIwMHHqeuxec9nuoT9Wt5DZTwAd5JeJjcxs9/OLwDh6EwOcb4RhYOLs2zMXR599Mb7w8HMri8ByRBw+Nh97f+a1Qw103WZThz0TXfdYre4TUylVllCkiL5mdO0qBzA4NcPARHroyZO13RCGaf+1ly73Au7XtH1VAAAOVklEQVTQmEpDLymoqzG+4wNXxtMHr4s7PnDlmssjIv7um/v/BWDfnrn4+oF3xtMHr4uvH3inIAzQgzAMTKRxLaLbt2cuDr3/yti1s7FybHamEYduvHLoQbLXgsVWt4nOaP7Ic6f66iYBwOCUSQATaZyL6Ma5iLDXY3WbETfLCzB8ZoaBiTSpPZBHbZJ3BQQ4HwnDwESa1B7IozboRh4AbIwyCWBiTWIP5FGrao9Wwow4wGYRhoGJd+T4/FjahY3rcbrRHg1gvIRhYKIdOT6/aqZ0fmExbr7nsYiIoQbEcT1OP0qcEQfYLGqGgYlWtYXxKPoNj+txAJgswjAw0cbVXUEXB4AyCcPARBtXdwVdHADKJAwDE21c/YavuWz3QMcBOD9YQAdMtHF1V3joyZMDHQfg/CAMAxNvHN0V1AwDlEmZBECoGQYolTAMEKOrTT5yfD6uPvhgvOnAfXH1wQfjyPH5Dd0fAMOlTAIgRlObPEkbeQBQTRgGaBp2bXK3jTyEYYDJoEwCYEQsygOYfMIwwIhYlAcw+YRhgBEZ14YhAKyfmmGAERnXhiEArJ8wDDBC49gwBID1UyYBAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIq1bbMHAACc3+6/8ZFoHD4VUzniXIpYuuHCePddb93sYUFEmBkGAEbo/hsfiR13n4rpnCJFiumcYsfdp+L+Gx/Z7KFBRAjDAMAINQ6fihRp1bEUKRqHT23SiGA1YRgAGJmpPNhxGDdhGAAYmXNpsOMwbsIwADAySzdcGDlWTwPnyLF0w4WbNCJYTRgGAEbm3Xe9NV5+/4VxNuXIkeNsyvHy+3WTYHJorQYAjJTgyyTrOTOcUvpsSun7KaVvtR17bUrpqymlbzf/3DXaYQIAwPD1UybxhxHx7o5jByLiaznnn42IrzX/DgAAW0rPMJxz/rOIeLHj8C9HxOea338uIvYNeVwAADBy611A9/qc8/PN71+IiNfXXTGl9OGU0tGU0tGTJ0+u8+EAAGD4NtxNIuecI6K2dXbO+TM557055727d+/e6MMBAMDQrDcMfy+ldFFERPPP7w9vSAAAMB7rDcP3RsSHmt9/KCL+eDjDAQCA8emntdoXI+L/iYhLU0rfSSn9WkQcjIh/kFL6dkT8/ebfAQBgS+m56UbO+YM1F/3ikMcCAABjZTtmAACKJQwDAFAsYRgAgGIJwwAAFKvnAjqAUt33ySfilU+9ELMLEQuzEds//oa47hNv2exhATBEZoYBKtz3ySdi220vxK6FFClS7FpIse22F+K+Tz6x2UMDYIiEYYAKr3zqhdixlFYd27GU4pVPvbBJIwJgFIRhgAqzC4MdB2BrEoYBKizMDnYcgK1JGAaosP3jb4iXG3nVsZcbObZ//A2bNCIARkEYBqhw3SfeEmdueUO8NJsjR46XZnOcuUU3CYDzjdZqADWu+8RbIoRfgPOamWEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxRKGAQAoljAMAECxhGEAAIolDAMAUCxhGACAYm3byI1TSs9ExA8j4mxEnMk57x3GoAAAYBw2FIabrsk5//UQ7gcAAMZKmQQAAMXaaBjOEfGVlNKxlNKHq66QUvpwSuloSunoyZMnN/hwAAAwPBsNwz+fc35rRPzDiPiXKaW/13mFnPNncs57c857d+/evcGHAwCA4dlQGM45zzf//H5E/FFEvH0YgwIAgHFYdxhOKV2QUvrJ1vcR8a6I+NawBgYAAKO2kW4Sr4+IP0opte7n/8g53z+UUQEAwBisOwznnP9TRFw5xLEAAMBYaa0GAECxhGEAAIolDAMAUCxhGACAYgnDAAAUSxgGAKBYwjAAAMUShgEAKNZGdqADAGCd7r/xkWgcPhVTOeJcili64cJ4911v3exhFcfMMADAmN1/4yOx4+5TMZ1TpEgxnVPsuPtU3H/jI5s9tOIIwwAAY9Y4fCpSpFXHUqRoHD61SSMqlzAMADBmU3mw44yOMAwAMGbn0mDHGR1hGABgzJZuuDByrJ4GzpFj6YYLN2lE5RKGAQDG7N13vTVefv+FcTblyJHjbMrx8vt1k9gMWqsBAGwCwXcymBkGAKBYwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFAsYRgAgGIJwwAAFEsYBgCgWMIwAADFEoYBACiWMAwAQLGEYQAAiiUMAwBQLGEYAIBiCcMAABRLGAYAoFjCMAAAxUo55/E9WEonI+LZsT3g+rwuIv56swexxTmHw+E8bpxzuHHO4cY5hxvnHG5ciefwZ3LOu3tdaaxheCtIKR3NOe/d7HFsZc7hcDiPG+ccbpxzuHHO4cY5hxvnHNZTJgEAQLGEYQAAiiUMr/WZzR7AecA5HA7nceOcw41zDjfOOdw453DjnMMaaoYBACiWmWEAAIolDAMAUKxiw3BK6ZmU0mMppUdTSkcrLk8ppf8ppfQXKaVvppTeuhnjnFQppUub56719YOU0kc6rvMLKaVTbdf57c0a7yRJKX02pfT9lNK32o69NqX01ZTSt5t/7qq57Yea1/l2SulD4xv1ZKk5h4dSSk82f17/KKU0W3Pbrj/7pag5h7emlObbfmbfU3Pbd6eUnmq+Px4Y36gnS805vLPt/D2TUnq05rZehxGRUnpjSumhlNKfp5QeTyn9RvO498Q+dDl/3g8HUGzNcErpmYjYm3OubEDd/E/gv4mI90TEOyLif8w5v2N8I9w6UkrTETEfEe/IOT/bdvwXIuLjOedf2qyxTaKU0t+LiB9FxL/NOf+d5rH/ISJezDkfbIaLXTnn/67jdq+NiKMRsTcickQci4i35ZxfGusTmAA15/BdEfFgzvlMSum/j4joPIfN6z0TXX72S1FzDm+NiB/lnD/V5XbTEfH/RsQ/iIjvRMR/jIgP5pz/fOSDnjBV57Dj8jsi4lTO+XcqLnsmvA4jpXRRRFyUc34kpfSTsfy+ti8i/ml4T+ypy/n76fB+2LdiZ4b78Mux/AaXc84PR8Rs80XHWr8YEX/ZHoSpl3P+s4h4sePwL0fE55rffy6W38w6XRsRX805v9h8s/9qRLx7ZAOdYFXnMOf8lZzzmeZfH47l/wyoUfM67MfbI+Ivcs7/Kef8SkT8n7H8+i1Ot3OYUkoR8YGI+OJYB7XF5Jyfzzk/0vz+hxHxRETMhffEvtSdP++Hgyk5DOeI+EpK6VhK6cMVl89FxF+1/f07zWOs9Y+i/g3/v0gpnUgp/buU0uXjHNQW8/qc8/PN71+IiNdXXMdrsn//LCL+Xc1lvX72S/frzY9WP1vz0bTXYX/+y4j4Xs752zWXex12SCldEhF7IuIb4T1xYB3nr533wx62bfYANtHP55znU0r/WUR8NaX0ZPO3fAaQUtoeEe+LiJsrLn4klvcF/1Gz7ORIRPzsOMe3FeWcc0qpzPqlIUgpfSIizkTEF2qu4me/3u9HxL+K5f8g/1VE3BHL/5EyuA9G91lhr8M2KaWfiIjDEfGRnPMPlifWl3lP7K3z/LUd937Yh2JnhnPO880/vx8RfxTLH/21m4+IN7b9/aebx1jtH0bEIznn73VekHP+Qc75R83v/zQiGiml1417gFvE91plOM0/v19xHa/JHlJK/zQifikifjXXLIjo42e/WDnn7+Wcz+acz0XEv4nqc+N12ENKaVtEXB8Rd9Zdx+vwVSmlRiwHuS/knO9pHvae2Kea8+f9cABFhuGU0gXNQvNIKV0QEe+KiG91XO3eiPgnadnPxfIiiOeDTrWzHymlNzTr5iKl9PZYfr39f2Mc21Zyb0S0VkJ/KCL+uOI6D0TEu1JKu5ofX7+reYxY7nAQEf9tRLwv53y65jr9/OwXq2NdxK9E9bn5jxHxsymlNzU/GfpHsfz65VV/PyKezDl/p+pCr8NXNf+P+IOIeCLn/HttF3lP7EPd+fN+OKCcc3FfEfGfR8SJ5tfjEfGJ5vF/ERH/ovl9ioj/JSL+MiIei+XVlps+9kn6iogLYjncXth2rP0c/nrz/J6I5QL+v7vZY56Er1j+5eH5iFiK5Rq3X4uIn4qIr0XEtyPi/4qI1zavuzci/re22/6ziPiL5td/vdnPZcLO4V/Ecv3go82v/7V53Ysj4k+b31f+7Jf4VXMO//fm+903YzmMXNR5Dpt/f08sd5T4S+dw9TlsHv/D1vtg23W9DqvP4c/HclnON9t+dt/jPXHD58/74QBfxbZWAwCAIsskAAAgQhgGAKBgwjAAAMUShgEAKJYwDABAsYRhAACKJQwDAFCs/x+29Umrm6NxpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6052ed3358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter(X[:,0], X[:,1])\n",
    "ax.scatter(X[outliers[0],0], X[outliers[0],1], color='m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recommender Systems:\n",
    "### Movie ratings dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'R': array([[1, 1, 0, ..., 1, 0, 0],\n",
       "        [1, 0, 0, ..., 0, 0, 1],\n",
       "        [1, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),\n",
       " 'Y': array([[5, 4, 0, ..., 5, 0, 0],\n",
       "        [3, 0, 0, ..., 0, 0, 5],\n",
       "        [4, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),\n",
       " '__globals__': [],\n",
       " '__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Thu Dec  1 17:19:26 2011',\n",
       " '__version__': '1.0'}"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = loadmat('machine-learning-ex8/ex8/ex8_movies.mat')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 943), (1682, 943))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = data['Y']\n",
    "R = data['R']\n",
    "Y.shape, R.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.2061068702290076"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# average rating of 2nd movie is\n",
    "Y[1, np.where(R[1,:] == 1)[0]].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Movies')"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAAHkCAYAAABSY8jJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztvX+031V55/t+eoJB8AdgDZKAhdRoVwQ5mhyEtHSlQQ063sLMOBOgKtzCYt2htkzx3ganzOqvdMbM7cXBNdVbFCdaG8BxLLJcjmcQpDo3KIfUA8FQTRpdkBMkVvxVmUlJuu8f388+eb777J+fz/78/D6vtbJyvp9f++fn+Tz72c9+NimlIAiC0Hd+pu0MCIIg5ECEmSAIg0CEmSAIg0CEmSAIg0CEmSAIg0CEmSAIg6A3woyILiWibxLRfiK6ue38CILQLagPfmZENAXgWwDeDOAggDkAVyql9raaMUEQOkNfNLMLAOxXSh1QSv0DgLsAXNZyngRB6BB9EWarADzFfh8sjgmCIAAAlrWdgVwQ0fUArgeAKUytOwkvGTv//CtOxgnf/an1XvXik0A/eW7x95FzXojl3/6fi/9rTl6r8NO95H2WeW1V9HPMvHDM/KSkHXOtrbyvft1z+NZjJwWfqfPtSsd2PHTMbC+T0Hmd/4WjLxxLJ6ZdY+BtpfNt/m9iHve1d9U8hc7bro3Nz8lrFX747IsW67Hss3h9/AQ/+Dul1MtDaffFZnYRgN9XSm0ufr8PAJRS/952/UvoNPVGumTx9+yheWxeOV0q7Sr3Cs0xM38Mc9NTwetytOfQ+8TC1g1YtX1X29lY5Ivq07uVUutD1/VlmDkHYA0RnUNELwBwBYB7XRc//4qTx36HOt6BnePnZw/NR9/bNgd2Ti/Jfwqzh+axsHXDkmOcmfljzt/mtanY8l71mT54e5rlMnGdT+0TVcsTyidgr8eyabgEWdU06qYXwkwpdRTAewDMAngCwKeUUt9wXX/KaX+/2IEO7Jwe+9vG6qvGO1uos6Z0Tlearg6acnxm/tiSvKeyeeU0Lr/yK0n3cA3olsPnjZ3j5TXrySZ4t6zdjdlD82Pl27xy2nqveYz/3rZiT1Te+T0hTc513sxvqD9sXjm9eL354YghRuPct3GH9/zM/DFvPnkarj67/FG7WUE/v8y5HNdrejHMTIUPMw/snK78wneBrqn+sQx9SCbUz9CGmcnwL7nti+T7Svm0I/1VDmkOLi1Cf5ldX2jzOTovR85/bvE8f56pJbjKcHTTuiXX6Gfwa8379PP1cVf5+f22Ybr5DDMPrrzrv/m9tnzz5/nandefq41sdWXer6/n5TKv421j/m9ilpmXgefHp/ny+/Tz9PVltR0zrdBQ0+y7tjaPubcMg9fMgHjjcCwp2t4QNBOzvGW0XX1PTk15Zv4Ytq3YU6p+bfnI3U+q0NV+c3TTOix7YLf1XJm2janzWM1sIoRZCl3tRD5yvoShZ6V02CrCRshDbgGd63kp79nEDzOBcqq1r4K1im1TtW3DsxzEGI5zdlbbs3h5tSCLKd/c9NRifabMhIWuLTsLGSLXUMg17ATcQ9i6ZnBz9Q2d/3vuvLjUfSZ1fOBEMxNqpUtDN6GfiGbWMm355JSZ+q8TLchc2lLVeuqi71NVY/uQ8E2m5GYihFlZFd43axYiZTjmwze01cQ4PPqwzZ7GXKt/xwgUl3Zms7/F1tnsofna3W7KtH1fNdHQjGsZXJMFrrSrpD8Rwiw0PndpM/o+3Tljp6eB4w1SxmbE0S+r76Wt+vLw+rl7r/9Latbl5pXuCQFerynOwKYjbmxeNLF1bV5nuy+3YOJ10jWtUpe1bJkXtm6IqlNf2lXSnwhhBvgrNdYDXr88rpeXaxSmUPB5T9vuN6mqrsfeH9J0tC+Tz29Ns2r7LhzdtA4LWzc4O6jNoGwTqHUMV8yyxmh5Ia3x6KZ1Xr8/rjmb7i4mOUwGqQLTpSHF+OCt2r4LW9aOa2Lm7zqZiAmAkBE6p5G6qqd+F1xDQnng581r68r/zPwxPHTTBdHDFmE4yAQAyqm4oUXWIaouOao6dKrKgZ3TQWHkO1+XIJ6bnloiyMpoMyk2zDoXvE86dUySDFaYzcwfW1RxYxYh6+Uttxw+b2wJxrYVexbP8aU8oY7uEoqh+1xLhfZt3IEDO6eXdALbcibbMfO3Xm5jpmUOC3RUDXM5kWkT5OnY0jaX1/jK6lriZD5338Yd1qGtj3Nvu8F5znwWXyDuwlx6ZD7Pt7TKVn7bsi3befPv0LlcgjlmOVPMcjXzncwxATARw0xg6VCSD4dcSzRcy3hihqVlhlu2e/Sx2OGrK91Yz/2Ysh3dtA5PXueO2GHLQ8pQnl9rexZfGvXKjy7V2MqaDWLqOObZ/Brdt/Sx2L7WNKH0U/o+r8eU98D1bBlmMg7snF5SSbyCn5lZPnZOfxlcjRvjBV1muLV55dIQOXpmz/eS8Xtc6cZMQLgwh273f/IO7/U24RP77NCzgLCh3jcja9MQNY/f+CHvc4G4tufXaMGl+99Ftz5svcdWJq4Rt43OX6j8M/PHxvrqmgevWTzuw5xYKGNWmRjNrE+kahYxGkWuL/+BndOLQ9G6/KlceTXrpYwG5runLo27TF7axNdX2ghFJZqZge/LEDvt39RXMrWD887l0nRyDmHmpqeCi9HN364vrS2/Li3STNOVB99Xnd9jpv3QTRc479P4bG6a2H4Saucm/dBs629txAgyXq8p7iVVyyuaWU109avbFimaYVNff7ONfOFtukjbdrYUqrwPopkxUiS+76uasrzJtC3EziaZM5Zm1AxbWWLKV0X7ND35+cyi7ctrmzmNOaeJFWS5o2e47Fmp+LQRnqfQCglX4EeNayjOKeNsXGUE4mpfmyDj74RZZ2XyMBHCLIQZF95VkdpWEuPqwV9Il53FlY6t4fXzzA58YOf0kvjvNsGZonGYAoaXZduKPZibnlocktkEj5n/1VcdX0NpO2cSO03vClekX4zUGP8xw0wTW10fOf85Z6RXniavO1teQ8N5G+b1T14XJxR8eyPYXIJcbFuxJ2qTGvOdWLV911gaZbS4wQozLulDSyq4HSRmDaGeofFhCx1tYmswM68uwamfv2/jjrH85LCz2ASmRpfFJxxtoZtdWoYtv3w20hZ1w6f57Nu4Y3F5mi89mx9XaJbWhq1tt6zdvVgG/aHR5Qi1T0jzDt3vqhtXiHANXw9r8z2M3TAGsNs8zefb6s02C53iHzdYYbZq+67Figg1BF+b6foa6mMLWzfUut4sxocJOP4SrXnwmiWamc3Fowq8vCnq/7YVe7CwdQP2bdyBu/eus2q9trq0aWv8Xq3R2PKyeeW01zWDP9u2aD4V1+J0ly0rZOOyRVvh94T6nmuIbi7eN8vK+53tXMwHHBj1R1se+DvIn8/rj5dNXDMKujABkIuyxvC+7ObUJyN2Dvo2ydAFZAKgZpqaNjd3ZYqFC7IUI3BV95PU+0P7Pfoo0wbm7lZV8UXIsJFrksFMpyvOtT5CQ11N2XdrkMLs5LVqbPYvFLtfX8uD05nXcpvPzPwxvPKjYedK/syQ3cNcp6e9v11aC++83CZlM9Y+M7N8UaDZ1vPx+rENc1NeFPN+HePKteWZLXZZrI1oy9rd1ufa2k+jhSffAk/j2rbPl7+56aklaZnx7zimyaNM/C/exr7ZQmB89tk3q+zLg2mvjLmXvyvAqC74agBfOcsGNpVh5oTTFX84PvyqkqeulEcYp8rKCRlmGphfA5s2xrFFzTS/ULGhrAEsSds166Sfa+bPdp3t+TYNzPclNY3lMZpBrFe3zovva66fxe1I3EfP1T6+Mvm0ENfzgHHNzOVvZ/PjCmlXMWHJYyKFxAyReZ3z/3154teFMK9zRf2wjUjMdghpwaKZwa2ZuYyvdRjLTcN224bf2PRzaza6HlKiRVRdLXB00zpcdOvDrWlovv7UdD8ok17VtZn8fv13jGZmPlvnXTQzhpbwLuOrrsBUw6Pvy2F2htgOlWu1gvmc2PRjXUNi0fXgqnvbBIBrkxObdmZ7sZY9sNtrQ+KUMTabm7iYMctsedLp5hJksZpVmfSqrs3k9+u/N6+cDmpaPGwQkD5ZMlGa2VBp8ms/M38Md+9dN1HuFCZta9nAUi2mL644ZZh4zSwU5dOFObtlnotZu2n7O9ZW4sqrz3YU47nuW3Np5tVXX9tW7BmL8OqzZfEIvQtbNwQjp4awtYfLdhj7jDLX+gSZOZuXGmnWlrZNCzQFl6nVlKVq/aWmH7PiJhbRzAz69oVz2ea6oD1wJs05NhdV661v/dnGxGtmnBQpH9PwXXJQ3LJ295h2ogWYucC4yZ2lTb8+HtDRpEt1mYucDtW+dbIxtCHIbH3N1f9ytr9oZsIgGIIGUoUua75V8yaamdAIXdmObZIFGdDsZrupNCVkByvMzOVHscZnX4jnMhMKKWgHTFs6PkOxmV/bNL2rTOakhjbWm9fwv3WaB3ZOL4mmYBr6tQOrTj8mCF/sshtXuWKGLrbJHD4sjnVW9k0KpW4VZz5fTxy53C588b9Sh28xdR5ahqT7RsrepfpvvhO8a7IoxKCHmdpRL1XNdTmONqXK+9Kx5a3shia2Z80emscl77zWOXmgO1xZ9wwzH7Y88GNNLk8qk1aVPhGTnu+auvpjFafZ2DIBfp9Gno4MM3F8EbPe74/DJT83Tmo/KvMa4PhGvDFG2LKGTZfm4uoAs4fmFyNrAG4NzLYw3taZ1jx4jXfncB2XzIdtaZg+FrPlXWxwRpemoNvNp7XFLAK3Yd5ne+ltC+Vty3t4OctoIjERR3JNRujnnD53JHhtzIYmZh8yl4XxGIOxDFYz++HTLxuLjNCXxcdl8lplEa+Zng7DXWXj1pzbsfWVLhvkmyamLng/MvvPxGtmPLRMSsjfOohxjNTYBEOIKoLDrBvb1963uNgX1ZVfm2LbSqkvG65F5GWPlcG1L0Po+VU0qdjgAjZ89tlUx+YY7dXEF+k2lsEKM/2Szswfc64LM1V/bXCN8XJOMUy70jcbzezwXIj4oh7YPMRNbH4+thdn88qRUZ+Xj6v83A7pM97rfNxy+LxFbc9Mc9uKPdYJDBu6HW1pmNfptjfbQJ/X99leGn4sNDnhy8OaB68JrqawDZN1PbmES6itXSsi+MSSjVsOn+esl9SRAlckbG3Lh94x5YxlsMJs88rRy+bacYertXwxLHC8MWzjfR1kzqZB8cbgX+bQEiUzzzp9bmg1A9bZggvyvJud4ZmZ5UvSt03nz8wf807z606vrzG/ulvW7sbmldOLZbnnzosXj5nXb145jVsOn+f9cvMJB1++NdtW7FmsA9t1c9NTS4SimRZ/Vohzb7thiR3u7r3rFsusZ3vNvqXR/UT3Nb63g34mv8e3h8Hc9NRin9F513XA28CGrheXALZ9YGOwCWedp7npqbGy6TKXHaIPVpjpL/SyB3aPfdW5Id0lZLSgMo2Q/Dkm+rjt66Zfblc65jHb5g/6xTNdIVy7fJudgU8SaHQ6vB7MF9jUFvV5veWcLcoqj356+ZVf8X7Zbed4Xen0+CSO73m3HD5v8Xwoaq5rJyx9ja3NzHrlHxx9/Za1u5e0g64//nHlgtecIdRCgJfHzL+tbPq8fq6vLmzoD4ztOAA8fuOHgvfb4M+8e++6xfrg+dJ1VdbWONgJAHMFQIpRuuyEAU+jrUkHrY2WXZdpy7ftS2lb+6nvDblbhI6H2ipUt74vu2/Nasx6VtM1wbY29snrjo0JMa3x8P9deW6z35htZ4tLFtKazKgqKeVxPbuzEwBEdBYRfYmI9hLRN4joxuL4aUR0HxHtK/4/tThORPRBItpPRI8R0Rti0jGnx1M24XWdt02v101Mmlz1n5uewrIHdi+5PnZ4YNO2+JBZf/1tsaZ8WoA+ZtrubNfaNAMOd53xudzY4MLKLKeOPpLiHG2LW8e1SF0WlwDTaJeVMhNAtuemTgDwttP5sm1xF1ppEJoQAtxOsVVnf9sYZh4F8F6l1FoAFwL4DSJaC+BmAPcrpdYAuL/4DQBvBbCm+Hc9gA/HJLJl7e7FBnVVkuuLwRvWnLmbmbdvMmJ+cavEsI/JI8c2hDTzGOu3M3to3tup9PAytnymgTdGYwwFVtT5M4dgwGgYHmOD488xcdm3gPCyKV1eHfrbzB/vI7qPzh6aX/JclzHeJ9z0udVXzUcJFRcuW+PsofmgHZHvXq/vM6lr6Vnjwkwp9bRS6q+Lv38C4AkAqwBcBuDjxWUfB3B58fdlAD6hRnwVwClEdEYoHf7CubQa3xIf/T/vTPs27lhUxW3Ge06ZL6pvxlV3sJAmor/K5nW2iRBberbNXvlLfcvh88Zm4XxOq+bMp21WN3aZmW24sm3FniXphzSHGCFss+e4MNPX/WDV9l3B5T+2NEx7nVk/vJ+Z5/QMahUO7Jx2CiyX7ZcTqxGaH4qYWeIQrU4AENHZAF4P4GsATldKPV2c+i6A04u/VwF4it12sDjmJTQVDWDMHqDhL42pmemGvOfOi60dPTQNH2LzyqVbvZlDN9uXng8DV181+nqaX2bbsNCmeYS8yret2IN77rzYaqzl2oS5lEznO8ZuyV8m3Ta2F8ymhdmGybH4hoQasx3N9PXQ1zZUA5aaBMz0zPY124inb+Zz9VVuW2OsK9GWtbsX03R9VELY1nmaq2zMcmp7nMY1W++jNWFGRC8C8F8B/Gul1I/5OTWalUiamSCi64noESJ65HmMllxwTcAnaPRSH5uNzPwSzswfw+M3fsj61Q1NfXNsu+aYS194vl2a0Mz8MZx72w1jzzj3thuWaCh8xyMzT7xutOblcnZd8+A1i1qHSzvlQ+5XfnRck4354nLtUJdD54vDTQkal68Stz1q9wPXi+IThuZHwrRF8g/SthV7lpSbzwb62pn3L5e9s4pW64LPptvucbk66XS2rdgztmRNP8uMhmz2ddOWGjNzbdLKbCYRnQDgcwBmlVK3Fse+CWCjUurpYhj5oFLqNUT0Z8Xfd5rXuZ5fJZ5ZjI9Lnct1ysxmVYnlZStLzvLVOTtXZVazCfq0jK4str6Su967PJtJAO4A8IQWZAX3Ari6+PtqAJ9lx99dzGpeCOBHPkGm4V/fWC9zwP61589c2LphzLkvREhNNr+u2ibnus6mifAJANfsp+uZ2rmVX2dqQabXua5X1/6f/F49aWI751pt4breJMYYbRJjQM9FyM8t1H+0FpnqqFonvtFIyvDe1ERjQkOFaGOY+YsA3gVgExHNF//eBuD9AN5MRPsAvKn4DQCfB3AAwH4AHwFwg+WZS9DGyrnpqeAMl4npPMltQau273J6v9vucX2ZTY9vfr/NSVIPmbnDoe1+PZy02ZJM9L38nI6KwYegW9Yed/XQ99y9dx1Wbd8V1AiXP3rSmIDmabnqhg+Jc2s2ekhuG2aGXEKAuCVOvut52WyTFfxZq6+aXxzWx+aljBBIwfceaZeUmHfN7Lertu8aK7vNBBJiYpxmfaRu/hEzfDDV7zqHpnyYqdMpu6FJbNwz29CWl1E/J5SPMk6zoXyUJSbdGKdRX99wnW97SKzTDzl+x5TvnjsvLtUmrr7S2WFmU/AhTIojpYlNxXc1pk9jinXajTFKm75bpp/ZzPyxKEFmm/Hdsnb32Po87hnO07R1Vl5G/WLqfNgmDAD7UNHmfuCql5w+S76ZQ01o6Mr7hq3fuWbBteDn9+TUvEJDQHN9siuvIbat2BPVJra64X22d64ZdcIXFJt2IR/mLKWtYVw2DJ/9wJa+z1Newx0h9W9z3Sd/wbQjaUx5bVP5ut7MIaj+f2HrhtIaZmgdpobnXQ8LQ/Ynjm4fWzvxwI7mC2MKVluatmf61uua/cAVqQQYvcwuYZI6u5fqf8fx2UFts+KcWw6fNzYzGfOumGm4zocYrDAzMRvBFV4nNDsGlNMGynZCU1sINbKeOjddCGK3mosx4oa+mvwZup5d5belF6orXZa56SlrubSmamsnfYx/7DQxL5BttYWJbZWBrjNXpN1U59oQpgYZWhHA82wro74/1Pfv3hs2b7j6YlX76KCFGZ/Zc+0ADYw3/Mz8sSXe9twBks8u2Wb5NK7gfK488rz4tCqbLxw/ptdlxnyJXbOJ5uwvN5ZffuVXcMk7r7V+cXUn1fUEjMIs+8pjc9Ll1+tlWLyc+mWxLY+aPTTvDc3NJ2hsdcmJEezmi3lg53E/M67l26JO6PQWtm5Y0l6uYXlI49bnbX03Bj3hw+/l94c89U1vAFuUjYtufThYDq5BxzJYYcbtZfs27nC+tMB4hdlm23hDbl45veg0a8aWMod/PA3XTJnPwMr/1o1rapgP3XTBkljyW9buXvIlNp0WddrmEHzfxh1Y8+A1i0JGu2potq3Yg4tufXjxY2AKUv0S6vLf/8k7vC4UtnqxtYFrQbrLqdeHbkfb8ihgXCibmILymZnlVk1DO5Caz+AOwTr9I+c/Z135YVuuFdJezMgXsUJM9xeb3ZNjlt/UZnlbm31Hc8+dFy9e5xrScg06lkHOZq4//0T1lj9/22LDbluxp7aZRKDc7GYMfIbLl0ZV58yY+13XlJ2l5RMLNuHNnxszM8vvdV0fmqlzYStjHbOSqe0Y2z/qoqk0J3o281uPnbT4peG2kZw+OCE7W6oPUygkiq/TpCyeTkVrWq40YnyzbPgmP2YPjUdXtQkm84vONVGX4EudqeM+fnyHKd/9VdwrUgSDGd2kjZUGdfa7MgxSmAH2BcA+I+jC1g3RHsx6iOkTjmZaXP22DUvKLO4Ned9zXNrTgZ3T3hmquempsfPaNsijyfqwrTcNYQvtc2Dn9Fi9mUNJvr7WNoNoet7H5MvMw9171y15tqu9Uo4DaS/97KHR3qYcs0/VIUR0vflWJNiG/Sau+137F8QyyGFmlbWZPlLU6rLOnFWdQF1DrFgn2pxOqDHp2tJr24HUl4fYPuC6vwtlC5Ezj7p9bfbf2M2EJ3qYeeScF1q3R3N9jfUQgq8Xs0UxsDlE8q8tP3b5lV8ZS0drDDwfpnMqd4Dlx7lmoZ9h+vLoZ3IjPJ8R4s6rek2lmW9gaQebmT+2mJbeqJXXjU7LNTPI0zXrAjg+28XbS8+mmetr+f08TR251IyGcXTTOm/ezPrl+TJDG+lnmX3A1md0GfQzeVuaM+e8L7gio+hn67bTz7T1C543rUmZx/X9Zr2aM/K8rnhezft4Guak2uVXfgWzh+bHom08fuOHFgUcMD5zqY/p8FspGr1oZsi/1GhmfjwOusanpegXpc5lT22Qw0gco1Xqa+rQfLrYJq4+Vjc+zT1HW9vab6I1M44vDpjpAwSM+x/xrwL/8nCbi8YMumez2Wmtybxe36PzYn4tbfnn511xu2Kw+VaZ9hfTwVhrB66t+GyG8hRbkVkeXm/8PL9XC7t9G3csaT/XEiGbFmLa1wB/dFcXtueY6fPfMbY4ftzWx1LyZ8uvL239XN/QUH+MQ33P1xfMII0p5RmkMDt5rVpUWfVaQ22wD3WGzSunlwSn035J+jyPc65Vc9fwkDN7aH7RHcHXKfQyGJ4ud6o014DaZrV8jqM6fzo//Hmrr5rH/Z+8Y6xTHjn/uUUhtW/jjsXt4/hQgWPzK7MtzQHGZ0N1mtyvb256aonQ1GXjs9S6vrn/GP9A6OfzPHAfO52WDo3O4W3K07SZLbiPmun/x/OkhZ2uf+5ozB1vTfhmLrYX3Sd4zfxy4W32IfMjGrO2WC91M4fSZj5cZdMf8i1rdy+2OQ88GmKQw8wzXnuqOveJN3nVXt9wxKZKm2q9OfQIqdj8+hS/Jd95XYbQ0CqUNkeXXT9zZt4eBSF2SMGft2Xt7iXPsj1n9tA8zr3tBmeatiEWz68m1Q/MZhMrE7hSPydmMoWXKWZmOJT3HBMPZd8bYFQ/egbTjMKRmh99b+wwc9DCrI/ECijbPWVJeUHNv6sSa48y0zR/x9aBvq4JO1gdabjqvmzIpzLkmvGOrZ+Jtpk9+z/9Qywbtlkv2zUpzytzr0uQ+ewQMS+xb6F5SJU3hz38b9tzdV758DCl7mzPMv32ygpTc4iaap9J8X1yaSSxz3MtwrfZ456ZWQ7A7cOV6rPle0aKIHPZKH39wWVXDjFIYXbei7+/OF3ssi2YzE2H94JMeYFMw3fMcif9vxlpAfBHtY0RFHptpk3I8sXA5hQ74A8fY9uoQud11fZdi/XvKv+5t91gbR99jO8CZU4C8PtCu0ppzIgQ5vraUF/hti1bfkOY7eh6nu1afty0S+oF+S5Bk2PWky+dsmErhysclS1qCT8HpEenGaQw2/OTl2H1VaMlMb69/lJW5nMjcOiLwV+02Dhc2s7ADcw+FfzopnXBeFv8bx6fjM/wzk1PjXlp84kO/QyfILZtd2aeN7Ugrs0dOf+5RSHC64U/V/u33XPnxWNCju/FsObBa7zxywAsTgjZSDF0266JGS7ZtNijm9ZZ7acm5oyyucjeNfS0zeia2GZUXdcf3bTOGTjA9Pbnq0RibLVmHmLDVmkmymbGjcRtrGWri6rl8XW0MhMTZQgZnYHjWoxtuVOsBmxLL2Q8d+Vpy9rdjdSNLw9VNC5zIqVsfnOWs4qf2SCF2fLVq9Qvf6c7O9oIQt3kXobWJSZ6AuC0F9qjgbqcIm0B9kL36t+uWFam86DN0Rawq/OuYS9ftmQ+2+Xo6FqiovPIy2PLC99OzubM6ysHz5+ZvunzpJck6WE8zxefnLHVK8+/+WyXMyj3ozJtjza7os1v0GWo1m1tPsdmqrC1m5lnnrZpC9bnbW4sOg+6j/Lr+RIum2FeH3OZVMy6NMvL6zc0XHT1ZdMuGmKQmtkZrz1V/e93/op1KBLrC5ZCjJpddkG19ttxXRc6bzqh2twSUnx+UuB5c9WRfm7MMNM3rONDUR/cUVXnzVb+VE0nVD883VzuGmWGmeaQ2hyelzHFuHwRY/GVY/bQPKbO2D/Zw8wz/90NnYpO4IoUUMdawhA5fZKasD+a+XXlv2lbaGp6sdc32SeqpNVUfU/0MBNwf6FDan0sqTHKj5z/HE6fO7LkeBsCN6dzJQ9eX3clAAAgAElEQVSoGIM5TA6tCQSW5teV/9DOQRxz1jfG78kktKO6Scg/L8WvquwaTHO4HOp/vjy5JmFiKVsGF4MUZue9+PvW4641b2UEii9GuWnnAUZ+UFX2BYyZWvcR4xRchZB7gc6/9gfLIcTNOrFtnuFC+2UBo7bctmJPcptUGS4+ed2xJf1E96mYzWhi0nY53abUPa8njq0v2YKf+uo0NCx3peNikMPMtkIANbFEJjd157mNYXSb9LEPpNCGW9PEDzNtuIaGtqUtNnT4Gw4fLpjPMdPhM4OxuGZBfc/Rs542Z0jzt2/jYhPfV9LceFf/rwWZa59S/beeEePneHvZlkjZAnC6ZnVd2NrUdZ2JLa2QtmHOxtrSdrWDqa342kvPYprD+ZCmw2dcbfCJAl9+be2Uas5x9X0XopkZmF8eU7MInXcdc0WC8N2TC9ezY4JF5nKsBI7PEC5s3YDHb/yQdRaTz0ryPLiu4Qvz9RA2VTMKzVyWbRvffWaaui20EHHtKMbLu/zRk7z5dkVnMZ/Nryvjr8bzZM6E5ujXopmVxHxRzYbQRl+fK4AZYA5w7+dorkGsghmDiueHo88ve2C3U2vhMcFy7QKky6+XMJn5mZkfxdXSQkmnZRra923cMbZmE1i6gDwE1zzMWHQmerOUnJh9QX9U9AYytnLwdbtAeHd1285Z21bsWTyuNR+9OcrM/LGxTZd98LrSda/zxvtIkyYG0cxqYtJsRXVSh52mr+3DtemqZehLHYhm5iFHOJQQq6+ad9qJchA7y5N7+jsH3M4YG9EkhpQ6jo2y0TW4WaCqIAoFWTSJCbDQJhMlzLRabU59p+wDGCtEzCUmoZcndWIgZsUBEB52pQrcKh3WVkaeP3PZVCj+lZmXfRt3eCc7bOmGlhb58p7aZlUwjehV3GzMJUsmB3YujTQze2g+aEvzhTMqQ2oZJ0qY6dj9piAwbQ8+p8Kywx3bZqgclz+PjZjQKLaOZ+sc/LqYJU0uHyjfi2GuH9S7KJl50nVghiEC7PXusleGruHp8r0MbDN1+retPkM2K/MZtvRtx3ztZLNLuXAJA+0PxvdA4Lj2xgDCM6h8IgGI07x9dtuJX5vZts1sYeuGxU0/gHZ8c3KQYlPJ4V8VW0/8uibDRfvykRPbjGKZcsbkr676y2mPE5tZzfhCQmuPco0r2mYZ6vTiN0npjKEVAK588+MuzciEL1t68rpy9WHLX+oQOnVdZiw2jc8ncFzai7m8yxbJo64PQcwqhtxMjDDzhTTxDdti7CK2Ts2nxavYVmK26UqNyBmL6+XmEUxd2GL22z4AZnkWtm5YfBH4Tuo8bdfQmIe48eXfFj3Xt8wnx14GKYIvtE2gzpN2inVpxKbpwmZiMbE59Nrq0WwTM2+5AnamIMPMBujLFHidSB3UR5Xhbtfb5cDOaRy48hYZZpYhdtF2iNBmJH3Hp23a6sc17LBd2+QsIaeLbiwxVNGC6u6bvraMqe+U/IlmJghCJ0ndBFg0s5po0lDfB9p2qBS6ic/emzo7PhHCLOVFqjrMdBm2y6YfUsVtkSJsM3OxdZBD6PDIC0c3rSvlnxZDyvrB3Glr6h4Sd+0jEJMfW5246jrnbOpECLOUaeKqbhQ5/Y64E6KLzSunrTOHsc6jJr7Nhm3YOinfo3PZA7uxeaU9KCa/Nja9mPT1CxdTd1XxecXn0M6bsrfGCmXfu6TrndeJDnvkqmvbao+yDFqYcS3J5/fEVV1bDDJ9PKXCyzYOf+m5B7TPTYJP5Vd107C9gKYg4rvmuMLUmPlxuUvY6skWjtoVXjsmWmnMc3hefEIoRjPR2jLvd7aX1rfbFc9jEyaL0FIlnT9bNFmNLXYdsHSlgS2e3uaVabHLbLQ2AUBEUwAeAbCglHo7EZ0D4C4ALwOwG8C7lFL/QETLAXwCwDoA3wewRSn1Hd+zq04A+Ka6Te9+F21vStHmlHvXp/s1OVYtDHm/yjbhfagPEwA3AniC/d4O4ANKqVcB+AGAa4vj1wL4QXH8A8V1Sdg8n134BJleXxizkYW5sDxWU+P55F9lV/51kD7+20asxhazoDuEbSF/KHqqeX1sHnxrHG3wyLK2xdSu9F35sAkyX39zaSiua33XlIlsUSe2dvYtjPe9EzGOwyataGZEdCaAjwP4YwA3AfjfAHwPwCuUUkeJ6CIAv6+U2kxEs8XfDxHRMgDfBfBy5cm4uGYIddEXrTM3bZa765rZfwTwOwD+sfj9MgA/VEodLX4fBLCq+HsVgKcAoDj/o+L6MYjoeiJ6hIgeeR5Lt3TzEfs1a+KrZ37ZYtLk17g0sCZnxcyvsC/tFC2sSVyaVcwLXdfyshjqqjuXPcwkpe1z07gwI6K3AzislMq6ElUpdbtSar1Sav0JWL6kQx3dtM65xk4PFfSwTsc9Mzdt0Gvd9DAiFA6FE7sBhKleu8LM6PLooa9Oj0ch5ejQz2VcNswOHCo7H6YvbN0wtnO2+Szb9nBb1u5ecj0futmGLvx6c5hnW2/Ir+XX+Vxr+MSHhm8aorn/k3dED4v1dSkvvbkpCrB0Axlbur7hN39ffHkJhT3atmLP2NpdbXIw8+zLj67nVEHY+DCTiP49gHcBOArgRAAvAfCXADYj4zDzPY9uHJsNrGrorWIs7nIIIDNvs4fmccvh85YIJG4bSqkL37V8xi+mfnzP4sOgHM9LwdbH2m7zKmUz79W/Y8qU633jdHaYqZR6n1LqTKXU2QCuAPCAUurXAHwJwDuKy64G8Nni73uL3yjOP+ATZBpd6ZtXuqMKpJD6DDO0TewERJmhbJVhptk5bT5hppE7pS70lLstfR1Jw+b64XqWia5XrpE8dNMFwXwtbN3gLUdKO9j6mOulb2oImtpfbe4S5u8Y4Wzeq59rarR1mGy65Ge2FcBNRLQfI5vYHcXxOwC8rDh+E4CbUx+c4vOVa6G5LcRMDGWm+XmoFx3fy8xjm0brLWt3O9O3RdiNzatLU3hmZrn3HgBjQ3Mbse2QYmoAwh7vbdkLc/YPLtD0c8+97Yaxa+pwZ2lVmCmlHlRKvb34+4BS6gKl1KuUUv9CKXWkOP6/it+vKs4fiH2+7kyhr1QdRstcERhiluTYXuhcy6lyPOPuvUvtSvxc2Ty4yuiz6+SO/utbMVHm+SGhkto+uRxuq74TNtsoEBe0MxaJmiF0jjZDYVelz3kHuukE3FmbWdPEhENuI46Va1lNzN8mqfHVQteZs4H8b/M5obRNe1lMXm3CoC9RSHTezXqKccg1sUXa9RFysI4hJMhy2LpsfSu1X9kYrDDTU79z01PZ45HnGIbybbnM9Zj8b/OaUKPrvJmdzra7tY1bDp+HfRt3OPNkDtltw0S+DlHviGV7li2fvAyclCGbr31863BzwQMEmIveXcMtW56WPTByUeHC3fdR4rvQm88Lfcx8dcY/SKFdxBa2bghGzeD54ztF8WtiVtmYDFaYzU1PWRvXhmtZiwvXHpj8vpgF1TEzTq6ZJc3c9NSYQHHZXHwLhFPy5ZtYMKflY2yVNk0g9uPjegF9dqeLbn148e+UOtHYfO7MPPH+lmpYD00cpc5S8ln9GGz+efyDZG6SYrJq+y6rzdKVPq8ffk1oa0YbYjMTBKHTiM2MUTY4o+s+2/FQUMUUj3CXjcq2K5LLJlh2KJyyKN9k9tDxobMOmRTycPfFI6uDlIXeqc+z/Tbhmp1Pk3cdq5tQ/cQsjo85boZbylHWQQszm0NlCO7g6rrPtqI/NBxwDU1DbhXc1rXsgd1LrjcdRPnQwHY8xJa1S9OInYgAxh0sbf5GJiFbni/qQkx+TLQtZmb+2JI8lZkISnW/4DYzPpSK2Z09J656Cw0jXf1YY24pyI9zTEfc1VfNV3YoHrQw45vKpnR6V4PqZ4TWp9mItVmE7EWmUOI2IOC4oDU7lF6bGSL0cpp2JpvnOBcKr/zo1NgaSDNfNkMvrwP9QrvyZbvfN+Pms8Vc8s5rnec0oX40Nz3lnWTg/YALBpuGVyYwo65nngebkHbVp2tPUk2ojsz8hmY/+Xnel8XPrIDbzPSM5uyheax58BprA9XhG1TFKZPfy9cc5gzDwtPg5bfl2+Z7FJuXMvUQCjbpS9u2tjQ2f6F+YDtvK5+vzGbedd2m1tPC1g14/MYPZV0D6cpj7DmNbV0oP+Z6hq4Ds55jbWaDFWZn3vxvO+f81xeqesfbFpDXvfC6i86eTdL2wvYYyi5+n/gJgMuv/EppJ8LQ9am2lbJG01hy+kvNzB8L+viEbGixLjGu+2Ov4+3QRUGW2i5VnLe7LMhMl526/PsGK8yAtJfKdh93OOQNYPu6+DpiSC3vUkfctmIPNq+cjo5vFZt313VV7i/zlW9yJUFM2XzRKvqIrX5jI4pUZbDC7O69S90YUuEOh6EG0FEBXIbusjs7uaauedn4hEXVMuuOx20WpjCem55KnkpPmY10BdE0KTOdz1dSuJYc+cit5edeWN4FbBMgqUvayjBYYaZninwzQrbOG+rgrpnRVdt3LS7f0b/5/SlfXe6uwKeuAXtEBr7EhC+B4qR0IJ+PnPYf872E/HotNFz5sn0k+Gwxnz01yxByE+CYfk229ohZCVJVy+dl8M30lRly2j6CuYnJF5+VNPsvMN6mvi32UhnsBMB7Ht24OGRKMTy2sXGDL0pqTN5j8pxiIDbT5Peawjkmf/yaKjOjJrlmSl2kGqxzRa41qdJ2XSeU39lD85g6Y/9kTwAA8YurOal+XlUxX64ya/Fi1jK61iGawy3Xl5Jvz8bztHnltFfDMOurSjBGkzJrK7U2lytSSsiWGnOvTwuZPTQfdGTl5NhMt0lC72hKnQ5amGlcjpK+mTgT3fnLvnhmB9N2hRxaYJmXWnPL4fPGOoytc81NT3kdhX0zicsfPanURyWGMtFQ+CL42NUMPnz+ZLH4tK7NK6eDkSpM+rQVXs7JgMEKs20r9ljH6xxbRbo6eFXV3cxDTifdmJc6pQ5S7g+xavuuMWFbRfCalAkTw4VMqOxV7E45nZtz1tmQGaww49pAn9TutrEJiJTAkLYJkn0bd5RaJxuizAcmRZsrIyxzMzc91StNq00GK8yA4y9aSgfO7QOTapupaxYqVqDn2M3KnK2qa5hZxu6VYn8qQ4qDdGxfG+rHOHeE58EKMz7MDHWa0Aa9nNiwwWX3D/RFrLBhy69tA9WUHY9CefKd44H8dMTVbSv2LAq4GP+7FMFrQ6cRcn0IpVNGqK++yv/cqpucxEbRjbk+hdjnuLYVdDnTxobGimHQrhlAeueJmQZvw33DR848H9g5jX0bd/Rier+sG4KrvlIXqZehzHrVSV93OvFrMzXmUqQQ3LvddZ9v2JoarC/2ft9zYl6I2KF26LpQEErbOd9aztx7msbgqq+qgiymbObeCq5AnJw6BFmZjUli2sqm8drarQ5zyuCFWeqXm3u3p3it63tNUo3IrhjwvpcsppPFvqR6OOh6ZsrLbgbgs+3xaW5kUTcxHwWbEI6hjI+gbbOYJnb/LiMgY/essC1/M6lD+x20MNNfvhSjbxVjtW0NZ+rzXC+Qa5t7IG7zh9gXc8va3UsmAXwvU8zmJ74lRKFn+KhrYkHnp8psZoyNySXEzdUXZQRPXZMGvr5QRduSsNkBtE0ipTOY23TxjlVlYXNV9PCvruUymrv3rlvykvH602szY9BD9rr8pOpyndAvbOwOQbb6iLWphvZGKFt3VWy6vn7uc+BN6etmneWwQQ92AuDld753cW/KspFHu0RK9NKq5FzfV3ZWt02amuDpal+rSu76m/gJgByaTIxx1nV9WVyb/JrGYk6Mcd9lhDXXZsbs/m5GxTDh583lUjntYzFRT0xi2ih1mZSrDqqUVd9bZZettogNpx5DSh0OXjOrg5gvat+iF7RFrq94qpbTNfcawc3E7wHwnkc3DlKF7wJDHR4J3WTih5mhl62O6e466dKSljKCrMn8t922XdnQNye2Ley6xmCFWYiqjoi5wyeHqGNIFONPVXfY5pDgCbmq2K7TW7elkqusKRv65kqzbiGTu//V8cEZ7DBT75s5qUz6EhihW7j2JI2xXU78MLNOYr6CZTbLsGELsWNi+8p1TZA1OfTr4wxg28S2T9k6csXvy6nxiWZWE0M2ktexA7wguBDNLAHzq1T26+OLVVX2y5dqU6mqXYS2fgOAJ6/z73Zl5sEVFsZ2rXmsTHlitqcr46MWc83RTUtXUHB4Xdq2ZGsLW9iknHXfRPlEM2sI8WuaDEJaa9f9D6v20zr6uWhmLWNqM00Jsqoe9m1rCEC/bVfPzCz35r+rgkz31zKbxHBWXzXfmvuGCLOaaNIAzwVn1ZclRyjnlIintuPLHz0pKg9dZNX2Xb3UwHV/zWHn5ZppkxM/IswS6KrGYBOcdefV98LahlkuAWV7TtUPQVfbCeh23urAt0VhbkSYJdCnL26u/T1zcHTTukY11S63U5fzlouyw8yqfW/Qwkzbf1JDM/Oggr6ZOBuxablmslLDbvP8+u6pYgsz88frxBffTB/nmlrZlRMLWzcseUlS2tU20xraXCW23X3tV5WYmdnQsSbywuHtrcNwmdhm7asK+lZmM4noFAAfBXAuAAXg1wF8E8DdAM4G8B0A/1Ip9QMiIgC3AXgbgOcAXKOU+mvf8+uezWzSu34o/mpdm80tW69V275r9RBLm/2w67OZtwH4glLqFwCcD+AJADcDuF8ptQbA/cVvAHgrgDXFv+sBfDh3ZlLV4jKduazqPQRBBjQ7vIoxOpet11Dbh9q5T4KM12Mf+mHjmhkRvRTAPIDViiVORN8EsFEp9TQRnQHgQaXUa4joz4q/7zSvc6WRqpn19WspCJNAlzWzcwB8D8B/JqKvE9FHiehkAKczAfVdAKcXf68C8BS7/2BxLBuxgqzt0DJ9YdJm7IRu0IYwWwbgDQA+rJR6PYCf4viQEgBQaGxJKiMRXU9EjxDRI8/jCIByoYd9Rs4yw8smQgXF3BO7zMQ2oeDDVrfmx6EuR9wuC82YvHXBQbmrlGnbNoTZQQAHlVJfK35/GiPh9kwxvETx/+Hi/AKAs9j9ZxbHxlBK3a6UWq+UWn8Clo+da3MIGRNemzdcmR2HYu5x5cO2T2foeXx2at/GHd5ruzaBkSJAyggb3ZYxfS6m3Wbmj03EiMCs6zLvbOPCTCn1XQBPEdFrikOXANgL4F4AVxfHrgbw2eLvewG8m0ZcCOBHPnsZR++rmLqjeZOYm6aW8eC37R9Z9gWIqSuex9DeleZGym1rIyntW6YvpLyEsVvZdS2cUx3keO/ams38TQB/QUSPAZgG8O8AvB/Am4loH4A3Fb8B4PMADgDYD+AjAG6ITURXUIq20/bLloumXoDY6CB8l6lclPl6162ZpRCTf/NjILgZdNSMMhEKcg+LmhhmTepsbN3l7sIQuQt5aANe7i7PZjZGmSFb7o5Td0dc2LrB+kJ32Tiei7oFeBeESJU8tBW9IoeNr0y5By3MYrEFE/QRM/xoamiwavuu7IKrjgCPrjqz1VPVl6FtQZ57eJo7VHXOPNg2reYmjlBd5HxPBi3MdEWGOgP/wvuGLil2n7KheA7snE5+GcwYVDPzx5ZEqYgVEKlr5Gx1y+tnYesGLGzd4KwzWz3xl0E/P+WFbjuEkK9/pKzj1JpV0yaEGHcbDW+/AzuncffecW0w9K643pMyH4RB28w0Kbazmflj2LZiT2eD6JmYNpWqdqQc5bfZeVxrGqvYhKrc21TE19ztUzbdusob+1zbdbH3TvyO5mUXmtdlcO2LkT6UT1/95Hph+mz0DtWfbAGYjkwAJMANpTEvURkbhtnBQy4MuYm1TWxZu9s7JNVuLrY6sH15u4ZvE5HU8Es2dDu7+kgZQdaE02xMWWPDYZntbk5EpEQiTkE0Mw8pGkIuzasurSRFazKv5WUrq32VqR/znj5obH3IY98QzcygzFfW7JS+tZ65hpBlXoSYr1qKzdAkdYVC6BmxhO6JCVzpI8estYlpAK+SPjDqc0Nx5E6hL2szG4NXiE9ImGqxS63XL3KZF7PskCtmiMHzUyYdc79PU2D5XibbOV3Xs4fmxz4Atg6a4sYRI+h5mqHrc8xau54ZI4BWXxUWVJtXTk+kplfmHRv0MLMLKwCaIGeeD+ycxr6NO7IZ8rWNrS+zw1Wpo//wvTib3HfTZhqo6/3wPVeGmTj+AuU2RPueV8WQaUbQ4M+rOtRI8Qb3LR7XewBojdGnbd1z58XYvHLa+/LFto2t/KF7ed70tfp/2/P09Sl1ZZbfdIeIuTcU5587v/K6jImtXwWbdhQSZFXCUdnaK4VBa2ZVvyJtuFPEpFnm61zFH6jK9Qd2TmPL2t2Nf837QB/zX9dEl++3+JnVuKGJILjoo4DqOjLMLEmsv02XyJkf39ZxIVzDnrr8zbrQDqnD/9Q8d6GMnC76DmpEmJXAtfavrYbmazPrnsb3Pd9ma9u3cUd0EEJN7Atctwkgpi5TtTDed2JnPJuA17lvPWwoGGcuyvRjEWYGc9NTi41oukWEXrI6ZplSX6h77rwYwNK8pxi1feXwvbz6HE/bjKSbg1A7+M7rvJX13wqlbdaPmcaR85+zXtu2BsbbSP/Nj+l6i/Gjy1GWMkP1ibGZHd20Dk9eV33X5LIGbZfhtC9rNkNw9wEb3JY0M38Md+9dFx2tJPRszVDWPdrqwlcHdSzi7xJiM7OQQ2isvmq+VAdxDU1T8hS7Ro9rYbFD36rr/0LChn/RTUEG2LUA/YV/ZmZ8gxoXR85/LmodYIg6taSY/PiEOkfn0yXA61pNYqOtQJCcidHMBCE3KZpgV7WkLuQrlAfRzBD/lWnbXiH0k5QhbdsCw0UX8pUrD4MWZmYEVhdtRyYVusEkLugeEoMWZtpOE9K8+AyTj65qcHXHu1rYumFs6Q1/6btaJy58AqsLWkoX6FubagZtM+uCPUAQhGqIzawEff0ixVBltik1BFDKdSkhgGzkdFRuw+nZ7HO5hrptrTTgWjyH122u4Akmg9bMNKKh1YfUrVA3opkxUl82MQSP0wc705C16iHQhB/aRAizVFwvaJtCrs2XtarA0ntnxlJ1w5iY++usz1B8sljK5LFtoe5aAphzQ2IXEzHMFISuMJTlaz5yLyuTYSYj5WuVS/sq8xxbpNlcxGpGObSaFKO+zeiesvt6l8m1sYuLtrUwGzPzx6IFWe4JF9HMOkST8d1zMQmahjCirYX8opkxQl+wVOOk73kubSFG20iN+9UGZtljV1m4yl9lu7eua2ZDo05BlkPLnAhhFtIcLrr14eTnuTYecRnLQx2h7Ca5Ji7VPXXjENceoWYeYycHXOUPlTkmfhowns+UaLltD9VEII/Iod3LMNOgj0O9Lvt6dTlvQjpttKcMM0vSF0HGv+gpu2iXoYr2ktLxY7YgK3N/XffloErd5tAqU8ve5Q/TRGhmoh1Uw1d/OTTZ2GdIO6YzhOi7opkxfC8AN0z7vnSpX0Hzi1c2soVvPwLbJhSuvFTxwPbVn00ImXlZ2LrBqwGYz3CVxcxHbq/ymOdVTbNKhJMyGqQWZGZ/McltOzTzattzwfWO+DZU8TERmpkmFEu+zq9Y2WdPojbC26nJ8ocmYcq6ofTJDluXq02VOpBNgCMnAJp6WfoaYln8yIS6iO3nMsyMxFaZoaFAmaFG2RDLMUMLW35zDRvqFmS2fLbtLjGJVNm+z3ed7pu2fpz7gz3xmlmIWM3E/MqkaFexW6m1ge6cug66pDUKk4EMM3u4nEkQ6qBPNjsbMsw0yDHzxVXluqKS9sUjPDWkT0q5qoYLKjtjGNOmMeXwXZPSD3P1hSYFmSsEUBMka2ZEdCqAs5RSj9WTpeqIZiYIwyGrZkZEDxLRS4joNAB/DeAjRHRr2cwR0W8T0TeI6HEiupOITiSic4joa0S0n4juJqIXFNcuL37vL86fXTbdXNS9G5KNvhvF24ivL0wWscPMlyqlfgzgnwH4hFLqjQDeVCZBIloF4LcArFdKnQtgCsAVALYD+IBS6lUAfgDg2uKWawH8oDj+geK6ZHIMM7UQi92aLid9d49wDXXa+DAIzdFk+8YKs2VEdAaAfwngcxnSXQbghUS0DMBJAJ4GsAnAp4vzHwdwefH3ZcVvFOcvISJKTjDDbKF2r+i7YBk6MR+uOl+yJuLdC0uJFWZ/CGAWwN8qpeaIaDWAfWUSVEotAPgTAE9iJMR+BGA3gB8qpY4Wlx0EsKr4exWAp4p7jxbXv6xM2jH0xQA/FOpYcRHz4apzvWJX3WzaoMl1oVHCTCn1X5RSr1NK/avi9wGl1D8vk2AxgXAZgHMArARwMoBLyzzLeO71RPQIET3yPI6Ufo74UDVL322BQneInQB4NRHdT0SPF79fR0S3lEzzTQC+rZT6nlLqeQCfAfCLAE4php0AcCaAheLvBQBnFekuA/BSAN83H6qUul0ptV4ptf4ELC+ZNaFpYiPVCkKI2GHmRwC8D8DzAFC4ZVxRMs0nAVxIRCcVtq9LAOwF8CUA7yiuuRrAZ4u/7y1+ozj/gGrQ03fShp2uzUhSNiQJPY/TZU140tq+S5Sp+yg/MyKaU0rNENHXlVKvL47NK6VKjRGI6A8AbAFwFMDXAVyHkW3sLgCnFcfeqZQ6QkQnAvhzAK8H8CyAK5RSB3zPFz8zQRgOuVcA/B0R/TwABQBE9A6MjPelUEr9nlLqF5RS5yql3qWUOlLY4S5QSr1KKfUvlFJHimv/V/H7VcV5ryAT7PTJNlVmpjF1IbS4hDRD2XouMyMcq5mtBnA7gA0Y+YB9GyPN6TvJKTaAaGblkYXkQtfIqpkVWtObALwcwC8opX6pq4IsFf5Fr/NrXfbZ2nZg3l+XL1Nd+wmkaIZV26EuLTTHc3PlrUwd+fpMbL5SN5PW/zfhe+fVzIjonUqpTxLRTY/xpWsAAB4LSURBVLbzSqnSS5rqRDQzQRgOuTSzk4v/X+z4JwiVSNn7Mye5ZirF2z8fVesy1mb2cqXU9yql1CCimQl9Q8KTu8k9m/n/EdF/J6JrCw/+wdJH3yJbnsvaZmLvS7W/uOq1rfrummYmgqw6sRMArwZwC4DXAthNRJ8jonfWmrNM8C22fC+ZDiDom8njL0ATRmZfeuYmwOZ58+XwvXR8WGd7qWx58q0/NJ+RKjhcW+iFttNzPYMze2geB3ZOZ5uxDa3DDJW9CZcZvc1bromVmPY03znznpiPYeqHokxwxp8FcCuAX1NKdXIOv65hZpvhh2OHIXXsJ5C6TV7KkKmO4VXZbf1i2reu4WDXh5llXXZylCvrHgBE9BIA/xSjJUw/D+AvAXxKKdXJhXVt2cy6EGu9ah6q3B+6V3zYhDLkFmbfBnAPRgLsoQz5qxWZABCE4ZB7AmC1Uuq3AewhohdVy5ogCH2ly8viYoXZa4no6wC+AWAvEe0monNrzJfQIXx+X6HO3QUDuA1Zm1kOn/2rbU+AWGF2O4CblFI/p5R6JYD3FseESNpu6Cr47GAh427dNrKyQqnJCKiTQmpb53aOjhVmJyulvqR/KKUexPHVAb2jDW2gzpe66j6TQHu7J8XMdPnKp4XS7KH5XmpbVfpi1X7c9gc292RZrDA7QET/lojOLv7dAqC3oXiamAKv4lvlwvWypmgZqWXPJSCqLOi//MqvBJ+tX4wuLC/S5euyfSlErrw3KTBjhdmvYxQx4zPFv5cXx3pJE50sVROLETK5h0b8xbd9JWfmj2VLk5cv9ZkP3XRB8Nlas+zCZiK6fDEhwVdfNb+kPzYpBF1ppXz0fFp9k644yU6zfaBt14wu+JsNiTocgYX+EOuascx3koju9Z1XSv1qasYmgc0rpzvv0d0nRJAtpe0PZtvp2wgNMy/CaKekr2C01+X/Y/ybSMoEsqsanDEVnce2DPsmfL2dqy6q2lds7dIFG1odtC1I1jx4zdhvcy1mrqFyynNCwuwVAP4NgHMB3AbgzQD+Tin1V0qpvyqdw55jaly2l3P1VfNLbENlXtayNgedR97pywoLn0D07dzE71191Wix87IHdpe2w4UEs00TzqHV1TFLaitL7g9PnbY309an6978PwZf/aY8xyvMlFLHlFJfUEpdDeBCAPsBPEhE74lOYQKIeTlXbd/V+rpEMyR27MsT0gJ85brl8HlR18Wcb0sbCc2mlsFWFl1XubRJ1wRE18KK55pkCk4AENFyAP8EwJUAzsZoH8uPKaUWfPe1SdsTAIIQomxkjy7RVOCAXBMAn8BoiPl5AH+glHo8U/4EYaLpuyADureBc8hm9k4AawDcCGAXEf24+PcTIvpx/dkThMmjz862bRKymf2MUurFxb+XsH8vVkq9pKlMCsOlrdnWtpfy+Fj+6EltZ6GXxK4A6C0xndZ8ocqGx67yYuaaMUuNtW/Ls1lm85qUcobq32YIj3l+Va/5mCGSmfcq7eu615bvnEPQ2DzHCveYflrHrHkMg18BoCuobb+cOokxxPquyekAaUunLgdLiVxbDz6H7zJt6XuebkNfW2aNNNs3+jibWWXFgLzUwpDJHWl2EHAV2abe53SONJ8fejYXZDHDXP48LsiaCINTxg/KVQ7b8ToN4FWfbavfHCGY6qKt9PT/TYZlmihhxm0RWnjwys5pqzAdFk+fO7L4d6iBuXByGYNXbd+15Dmzh5auOgDyd+hYr/qYjmwrX2i7uypUXS9r1u/RTeuSnGqbXq9bR3oxHvv6f11fM/PHgtvPVUWGmRNI2w6bfFgcM7yuYxg9hEAAOcrQh3qQYWZPKaOBpN7jE2RVY/qbLGzdsOSZXDDFxPzSBuIc6Od04QWuqm3G1F2ILtRDLkSYIX1cH9MJy758ZTpXzg6ZI6Y/n2IPaYC259naw5Wu7Vpf3evndCGaSNV2S9FW64hI4sJmU23Cr28ihFmoIVKHXDGdsOywqEyjVzWy5g7xbU7dp760Me2h88BtMhpf3evrbO4FfdxDQBNqw6rDdN6GoXqy2VSbmG2fCJtZH+wCQ0W/ZOI6IpRFbGaMfRt3tJ2F3uMalsVoBCLIhBA5ZqsnQpiZUTG7Sp12HN+zY9LVw7IDO6ejh3TmtbEMITpsl9d+dhFz5FSm/iZCmOWY9WmCskt+Yr5qvmenRKNdfdV8tKalO2iqkK4z5n9TQiaXNjopQjFHOSdCmJkRVjWuCiyj8jatTWgBMTN/bMlXrYqGl3tIODc9NRZtlsPrOab+cnT4vgx5dRv2Jb+alIAG/LhZTv07pc1lAkAQhE4jEwCMLgiy3EuK2nYjSC1Pyhe27bLVxRCGjLFtE+ofdYxkJkIzmySGtmFu20uvhPYRzWxCcQmyroZizu3Q3Dfa2F9yqNQmzIjoY0R0mIgeZ8dOI6L7iGhf8f+pxXEiog8S0X4ieoyI3sDuubq4fh8RXV02P7kMzGWHQE1OENhCCHVhqG2rA54v1wbBobrT5a0afTYEz0eZfmDrXzHt4iq/bd/Ktigz/Ay1T2r71TbMJKJfBvD3AD6hlDq3OPYfADyrlHo/Ed0M4FSl1FYiehuA3wTwNgBvBHCbUuqNRHQagEcArAegAOwGsE4p9QNf2pM8zBSELpEj4knrw0yl1JcBPGscvgzAx4u/Pw7gcnb8E2rEVwGcQkRnANgM4D6l1LOFALsPwKVl8pNDMxJVvjx11J20R/dp0rWkaZvZ6Uqpp4u/vwvg9OLvVQCeYtcdLI65jifz5HX9n0kSxnENrYY6G9pHmvzgtDYBoEbj22xjXCK6nogeIaJHnseRJedz2BTatkv0mSbrbuiTBn2iyXZvWpg9UwwfUfx/uDi+AOAsdt2ZxTHX8SUopW5XSq1XSq0/Acutidf5lai6rlIbh1N8kcr4LeWog7rzKAyfOt7FpoXZvQD0jOTVAD7Ljr+7mNW8EMCPiuHoLIC3ENGpxcznW4pjpahzjWZoXWXsS+1aemVr/Jz2CN9+m+awLSVdV3lC6XaBsgEJuxD40STXByi1vfT1ps06RmNLrcc6XTPuBPAQgNcQ0UEiuhbA+wG8mYj2AXhT8RsAPg/gAID9AD4C4AYAUEo9C+CPAMwV//6wOFaKqi8WPx+qaPNZXADY7tXnXY1cZf1lzBS+S0DNTU9VGraZi81NAbFtxZ4l98RuwszPaYE7M38s6oWL0YRdO2bZ4P59txw+b8kHIKcm4mt7Vzpz01PRtsSHbrqg1DnAXU+8fnQeQ304NdrNxK0ASJ0qLruus697WdrynWtta9/WyJbZ8Lbudm+7X9W1obMP2QS4Y35mVZcZ9U0QCM3S1WVsOZajiTALCDMRDoLQD1p3mu06IsiESaANn7u2IgVPrDAbEn3whHflsUt55y9+lVnWnDO0VZ/Vts+dzn/MhFXVvjDRwqzOF6kul4PUWVQfC1s3WK/3pTEzfwyzh+YX77Xdn5IHm4acGq3URcjFwHze5Vd+ZfFvvuN6aho5DfTms1Jm1NuC2+50/l2TBrw8VUdLE2szE5pFbJRCWcRm1jLmF7tqeKGmnEvrSqesICujPcf4ppUtZ5eGxSl02Tk5FyLMDHLFWHJt0ODDZ99w3X9g57Q1T7H5tDn38uFLjqFMlRfJJwRdz/Xds/zRkwAsLWeO/MSQe/f4WPro8wik9T8ZZjZEV/2A2qZrw8+68iPtv5TYupZhZsuYX5RJ78gujaRru83XJVgntf19ppPcdS2amSAEaGMJj3Ac0cwSkGB+eeFa2BAMzyLI+oEIMwBHzn+u7Sx0mlSBxI3NtqgYTeVDyE8u37Y62lKEGY7PcIWImVnq6gsXu8Qkd9y0zSunsbB1g7PubC+HS1POMSM3dC3cLF/uAJ65tNQ6ZlfFZoZmwqp0bdauK9Qxy9d2mBwhL2IzSyAUtDEHIsjs1DHLNzc91VvnVqE8IswggqYKXRUak9amXTVvNIkIMwc2G5P54vrsUG3ZZpoOv1JVaAzdhtUUMqwWYTYGF1a24Y/54vqGSG2FXumbc2aueorVTJrQYERLOs5E7JvZRdocmnT5BaiSt7rC1JgviamZuNJtQoMRLek4Q943cxDUMTTSL0DVqXTXS1xVWLruD315czqc2mJfufIVEz8rJp0mSInkUVfecj035jmzh+ZLRZUJIa4ZgiB0GnHNEARhohBhJiTTVXcMoX3anJ0WYVaCLhvrNVUETujeSfPh8tGHvtAkrtnpJvYrEGFWgq7NVtmETxWBU4ewqluba0tb7FpfqJOu+wSKMDPo45c2t/Cpow7qEJBldvaRIXJ5qvgENhFGSYSZwSR9aV3krIMywiNWAyiTTxkiDxcRZiWpsjFF37WDqvticmyCi2sAfa+rpujbiKKOdhVhVpJUrYC/1HVqB661mbHaTuilWNi6AXPTU97reEetula0SU2qjJG6K3YkX39sUtDFCqlQu5YRdiLMLDS9WNtHaqM+M7PcejzW3hES0joqr+863lH1WlFXOboU5feWw+cl39PWGtwUmjSdmEKqrCAt8xGTFQA9pg9BCJvKowS/bIY2+pysAMhAE74xsdi+cHV0qhxDEv6M0JA01xBIBFkzdPnjKcLMQ5Xp5NwGzqY6UY50UnZzt53r0jBfqIaEAOoQZQ28oimUp28x2QQ3EgKoQ3TJQN1VfMPxvrkMCP1FhFkA0bDC+IbjXbWxdMWlQsiHCDNhotBaZB9cKoQ0RJj1BJ8mUbfBfEhDxSbWCArtIMKsJ/g0iboN5l0YKsqyJiGECLME5IUqT9W6E9ulEEKEmSAIg6A2YUZEHyOiw0T0ODv2fxPR3xDRY0T0l0R0Cjv3PiLaT0TfJKLN7PilxbH9RHRzXfntM32waYlmJdRNnZrZDgCXGsfuA3CuUup1AL4F4H0AQERrAVwB4LXFPR8ioikimgLwpwDeCmAtgCuLa1thy9puOnPOTU/JEFiYeGoTZkqpLwN41jj235VSR4ufXwVwZvH3ZQDuUkodUUp9G8B+ABcU//YrpQ4opf4BwF3Fta3QBUO4C9F8hLrp+gigTZvZrwP4b8XfqwA8xc4dLI65jk8cZkdK0cS4W0eoQ1YJOum7pqrmmDvfQ6FJjbzMxzxmg+pc1BoCiIjOBvA5pdS5xvHfBbAewD9TSiki+k8AvqqU+mRx/g4cF3SXKqWuK46/C8AblVLvsaR1PYDrAeBEnLTul+ht9RRKEIRG6WwIICK6BsDbAfyaOi5JFwCcxS47szjmOr4EpdTtSqn1Sqn1J8AeoHCoSJSJ4SK20HgaFWZEdCmA3wHwq0opvoL7XgBXENFyIjoHwBoADwOYA7CGiM4hohdgNElwb5N57gMSZWK4iC00njpdM+4E8BCA1xDRQSK6FsB/AvBiAPcR0TwR/b8AoJT6BoBPAdgL4AsAfkMpdayYLHgPgFkATwD4VHHtoJFF0IKQTp2zmVcqpc5QSp2glDpTKXWHUupVSqmzlFLTxb//g13/x0qpn1dKvUYp9d/Y8c8rpV5dnPvjuvKbQozqbxqkUwQUX7rUVcN2rnzZntPE0MqWbpvD9a62cypm2zU5TJY9AGpi9tB8K4ua+xYLv0v5bavNJpXY+o6dABBhJjRCl4SW0Bw52r2zs5lDINa3Kgepw4+mhkq+fFW1+TVxf53DH9uzdX31dXYyNt/mdU1+wEQz6xhVv2RHN61rZHazqXQEQTSznlL1S9aUgBFBJnQNEWZCEk3sJVp1Zi8mj5Pi/pJjlvTAzumkdm9rv1kRZjXh60R9tZsAwC2Hzyt1X0qZqy7oj5kha3oPAF3+pts+R3CE1VelzfLya33lzS30RJjVxN17xw3xvFH7PKtnliuWPpc5B7r8Q6+HlAmA3G4wMgEgiH+V0GlkAkCIRgRZN5CAAdUQYSYIHUFmiKshwqyDpBqJbdfLV34Y9HmyqGlEmGVGz2JWmRJPNRLbrje/8rmE24Gd04NZFG3SdLliBFXfJwyarFOZABCEBpmZP9bpvSS6iEwAtExbjoNCtxFBVh8izIRBoT8iYmvqFgtbN/R7Q5O2kGGmMHQmabgqw0yhFLk1mtC6vipf60lZX2ljUgRZCqKZNYSEzBGEcohm1jFEkAlNU6fdsIsTXCLMBkDuTpu7o/bVKN/FFzaFOn3UurgEToaZLSALu/MjewwMFxlmdhgRZPkRQSaIMBMEYRCIMBM6Q+w6vqZtb0Ndizo0RJgJnSHWd6rpIaX4dPUDEWaCIAwCEWaCIAwCEWaCIAwCEWY9oW8Op4LQNCLMesJQ/agkvLeQCxFmQqvImlUhFyLMBEEYBCLMBEEYBCLMBEEYBCLMBEEYBCLMBEEYBCLMBEEYBCLMhFrpe7RWoT+IMBNqRQJRCk0hwkwYBJO+3GvSyw+IMBOEQTDU5W4piDATBoG8zN2miWi9tQkzIvoYER0mosct595LRIqIfrb4TUT0QSLaT0SPEdEb2LVXE9G+4t/VdeVXEIT6aCJab52a2Q4Al5oHiegsAG8B8CQ7/FYAa4p/1wP4cHHtaQB+D8AbAVwA4PeI6NQa8ywIQk+pTZgppb4M4FnLqQ8A+B0AfMPOywB8Qo34KoBTiOgMAJsB3KeUelYp9QMA98EiIIVhI2GCuknXNnpp1GZGRJcBWFBKPWqcWgXgKfb7YHHMdXywiF/WUiRMUDfp2kYvy5pKiIhOAvBvMBpi1vH86zEaouJEnFRHEo0gflmCUI4mNbOfB3AOgEeJ6DsAzgTw10T0CgALAM5i155ZHHMdX4JS6nal1Hql1PoTsLyG7AuC0GUaE2ZKqT1KqRVKqbOVUmdjNGR8g1LquwDuBfDuYlbzQgA/Uko9DWAWwFuI6NTC8P+W4pggCMIYdbpm3AngIQCvIaKDRHSt5/LPAzgAYD+AjwC4AQCUUs8C+CMAc8W/PyyOCYIgjEFKqfBVPeMldJp6I13SdjYEQcjAF9Wndyul1oeukxUAgiAMAhFmgiBkY2HrhizPKeNbKMJMEIRsrNq+K8tzyvgWijAThAlmSKGDRJgJwgQzpGgjIsyERnBpAEPSDIR2EWE2ALq24NeGSwMYkmYgtIsIM6FW+iBohe5RJuCCCDOhVroWWUHoB2UCLogwE8bIrUl1TTMTG91wEWE2AO7emy94YW5NqmuamdjohosIswEgL6ggiDATGqJrw02hGZpsdxFmQiN0bbgpNEOT7S7CTBCEQSDCbEKQHY6Epmirr4kwmxC6usNRW7a0Ltvw+u4+0lZfE2EmtEpbtrQu2/BkdrocIswEoQOIGaA6IsyEWgkNmfo+pMpFV80AfUKEmVAr+zbu8J6XIZWQCxFmQq2EFgx32RAv9AsRZkKrdNkQL/QLEWaCIAwCEWbCIJDhaj76OikjwkyonVx7KZrwly7HcFUE4oickzJ1tb0NUko1llhTvIROU2+kS9rOhiAIGfii+vRupdT60HWimQlJiPYidBURZkISMcM5c2jR5FBDmFxEmAnZWbV9l/d3E4gGOXmIMBMGifivTR4izARBGAQizARBAND/obkIM0EQAPR/aC7CTJgo+urdLoQRYSZMFBJyaLiIMBMEYRCIMBMEYRCIMBMEYRCIMBMEYRCIMBMEYRCIMBMEYRCIMBMEYRDUJsyI6GNEdJiIHjeO/yYR/Q0RfYOI/gM7/j4i2k9E3ySizez4pcWx/UR0c135FQSh39Spme0AcCk/QES/AuAyAOcrpV4L4E+K42sBXAHgtcU9HyKiKSKaAvCnAN4KYC2AK4trBUFoma6tplhW14OVUl8morONw/8KwPuVUkeKaw4Xxy8DcFdx/NtEtB/ABcW5/UqpAwBARHcV1+6tK9+CIMTRtdUUTdvMXg3gYiL6GhH9FRHNFMdXAXiKXXewOOY6LgiCMEZtmpknvdMAXAhgBsCniGh1jgcT0fUArgeAE3FSjkcKgtAjmtbMDgL4jBrxMIB/BPCzABYAnMWuO7M45jq+BKXU7Uqp9Uqp9SdgeS2ZF7rB7KFuDW+EbtC0MLsHwK8AABG9GsALAPwdgHsBXEFEy4noHABrADwMYA7AGiI6h4hegNEkwb0N51noGJtXLjU8d80YLTRPna4ZdwJ4CMBriOggEV0L4GMAVhfuGncBuLrQ0r4B4FMYGfa/AOA3lFLHlFJHAbwHwCyAJwB8qrhWEMbgxugmBVvfo7MOCdkEWBCETiObAAuCMFGIMBOEBjm6aV3bWRgsIswEoUGWPbC77SwMFhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMAhFmgiAMgkEGZySi7wH4KUYhuYfAz2I4ZQGkPF2mi2X5OaXUy0MXDVKYAQARPRITnbIPDKksgJSny/S5LDLMFARhEIgwEwRhEAxZmN3edgYyMqSyAFKeLtPbsgzWZiYIwmQxZM1MEIQJYnDCjIguJaJvEtF+Irq57fzEQERnEdGXiGgvEX2DiG4sjp9GRPcR0b7i/1OL40REHyzK+BgRvaHdEiyFiKaI6OtE9Lni9zlE9LUiz3cXO9Sj2MX+7uL414jo7DbzbYOITiGiTxPR3xDRE0R0Uc/b5reLfvY4Ed1JRCf2uX00gxJmRDQF4E8BvBXAWgBXEtHadnMVxVEA71VKrQVwIYDfKPJ9M4D7lVJrANxf/AZG5VtT/LsewIebz3KQGzHahV6zHcAHlFKvAvADANcWx68F8IPi+AeK67rGbQC+oJT6BQDnY1SuXrYNEa0C8FsA1iulzgUwBeAK9Lt9RiilBvMPwEUAZtnv9wF4X9v5KlGOzwJ4M4BvAjijOHYGgG8Wf/8ZgCvZ9YvXdeEfgDMxesE3AfgcAMLIEXOZ2U4AZgFcVPy9rLiO2i4DK8tLAXzbzFOP22YVgKcAnFbU9+cAbO5r+/B/g9LMcLyhNAeLY72hUONfD+BrAE5XSj1dnPougNOLv7tezv8I4HcA/GPx+2UAfqiUOlr85vldLEtx/kfF9V3hHADfA/Cfi2HzR4noZPS0bZRSCwD+BMCTAJ7GqL53o7/ts8jQhFmvIaIXAfivAP61UurH/JwafRo7P/VMRG8HcFgpNZTdbpcBeAOADyulXo/RMrkxW2xf2gYACtveZRgJ6ZUATgZwaauZysTQhNkCgLPY7zOLY52HiE7ASJD9hVLqM8XhZ4jojOL8GQAOF8e7XM5fBPCrRPQdAHdhNNS8DcApRLSsuIbnd7EsxfmXAvh+kxkOcBDAQaXU14rfn8ZIuPWxbQDgTQC+rZT6nlLqeQCfwajN+to+iwxNmM0BWFPMzLwAI8PmvS3nKQgREYA7ADyhlLqVnboXwNXF31djZEvTx99dzJxdCOBHbMjTKkqp9ymlzlRKnY1R/T+glPo1AF8C8I7iMrMsuozvKK7vjJajlPougKeI6DXFoUsA7EUP26bgSQAXEtFJRb/T5ell+4zRttGuBgPn2wB8C8DfAvjdtvMTmedfwmiY8hiA+eLf2zCyTdwPYB+ALwI4rbieMJq1/VsAezCamWq9HJZybQTwueLv1QAeBrAfwH8BsLw4fmLxe39xfnXb+baUYxrAI0X73APg1D63DYA/APA3AB4H8OcAlve5ffQ/WQEgCMIgGNowUxCECUWEmSAIg0CEmSAIg0CEmSAIg0CEmSAIg0CEmdAJiOhsInrcOPb7RPR/tpUnoV+IMBMGC/NoFyYAEWZC5yGi3ypivT1GRHcVx04moo8R0cPFAvDLiuPXENG9RPQAgPuJ6Awi+jIRzRfxuy5utTBCbciXS+gDNwM4Ryl1hIhOKY79LkZLa369OPYwEX2xOPcGAK9TSj1LRO/FKJzNHxfx7k5qPvtCE4gwE7qCaymKXub1F0R0D0bLiQDgLRgtaNc2tRMBvLL4+z6l1LPF33MAPlYs5L9HKTWfP+tCF5BhptAVvo/RmkfOaRgFA/wnGK13fAOAucIWRgD+uVJquvj3SqWUjmz7U/0ApdSXAfwyRtEfdhDRu2suh9ASIsyETqCU+nsATxPRJmC0/wFGcbb+B4CzlFJfArAVoxA0L8IoAupvFpEfQESvtz2XiH4OwDNKqY8A+ChGAlEYIDLMFLrEuwH8KRHpMEh/gFHImi8R0Usx0sY+qJT6IRH9EUYRbR8jop/BKLT12y3P3Ajg/yKi5wH8fZGGMEAkaoYgCINAhpmCIAwCEWaCIAwCEWaCIAwCEWaCIAwCEWaCIAwCEWaCIAwCEWaCIAwCEWaCIAyC/x+aUbu0xKMmcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6052cda390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.imshow(Y)\n",
    "ax.set_xlabel('Users')\n",
    "ax.set_ylabel('Movies')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Collaborative filtering learning algorithm:\n",
    "### Collaborative filtering cost function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(params, Y, R, n):\n",
    "    Y = np.matrix(Y)\n",
    "    R = np.matrix(R)\n",
    "    nm = Y.shape[0]\n",
    "    nu = Y.shape[1]\n",
    "    \n",
    "    X = np.matrix(np.reshape(params[:nm * n], (nm, n)))\n",
    "    theta = np.matrix(np.reshape(params[nm * n:], (nu, n)))\n",
    "    \n",
    "    J = (1/2) * np.sum(np.power(np.multiply((X * theta.T) - Y, R), 2))\n",
    "    \n",
    "    return J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 10), (943, 10))"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "paramsData = loadmat('machine-learning-ex8/ex8/ex8_movieParams.mat')\n",
    "X = paramsData['X']\n",
    "theta = paramsData['Theta']\n",
    "X.shape, theta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(22.224603725685675,\n",
       " array([ -2.52899165,   7.57570308,  -1.89979026,  -0.56819597,\n",
       "          3.35265031,  -0.52339845,  -0.83240713,   4.91163297,\n",
       "         -0.76677878,  -0.38358278,   2.26333698,  -0.35334048,\n",
       "         -0.80378006,   4.74271842,  -0.74040871, -10.5680202 ,\n",
       "          4.62776019,  -7.16004443,  -3.05099006,   1.16441367,\n",
       "         -3.47410789,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ]))"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nu = 4\n",
    "nm = 5\n",
    "n = 3\n",
    "X_sub = X[:nm, :n]\n",
    "theta_sub = theta[:nu, :n]\n",
    "Y_sub = Y[:nm, :nu]\n",
    "R_sub = R[:nm, :nu]\n",
    "\n",
    "params = np.concatenate((X_sub.ravel(), theta_sub.ravel()))\n",
    "\n",
    "cost(params, Y_sub, R_sub, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Collaborative filtering gradient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(params, Y, R, n):\n",
    "    Y = np.matrix(Y)\n",
    "    R = np.matrix(R)\n",
    "    nm = Y.shape[0]\n",
    "    nu = Y.shape[1]\n",
    "    \n",
    "    # unroll parameters\n",
    "    X = np.matrix(np.reshape(params[:nm * n], (nm, n)))\n",
    "    theta = np.matrix(np.reshape(params[nm * n:], (nu, n)))\n",
    "    \n",
    "    # initializations\n",
    "    J = 0\n",
    "    X_grad = np.zeros(X.shape)  # (1682, 10)\n",
    "    theta_grad = np.zeros(theta.shape)\n",
    "    \n",
    "    error = np.multiply((X * theta.T) - Y, R)\n",
    "    \n",
    "    J = (1/2) * np.sum(np.power(error, 2))\n",
    "    \n",
    "    X_grad = error * theta\n",
    "    theta_grad = error.T * X\n",
    "    \n",
    "    grad = np.concatenate((np.ravel(X_grad), np.ravel(theta_grad)))\n",
    "    \n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(22.224603725685675,\n",
       " array([ -2.52899165,   7.57570308,  -1.89979026,  -0.56819597,\n",
       "          3.35265031,  -0.52339845,  -0.83240713,   4.91163297,\n",
       "         -0.76677878,  -0.38358278,   2.26333698,  -0.35334048,\n",
       "         -0.80378006,   4.74271842,  -0.74040871, -10.5680202 ,\n",
       "          4.62776019,  -7.16004443,  -3.05099006,   1.16441367,\n",
       "         -3.47410789,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ]))"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J, grad = cost(params, Y_sub, R_sub, n)\n",
    "J, grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularized cost function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(params, Y, R, n, lam):\n",
    "    Y = np.matrix(Y)\n",
    "    R = np.matrix(R)\n",
    "    nm = Y.shape[0]\n",
    "    nu = Y.shape[1]\n",
    "    \n",
    "    # unroll parameters\n",
    "    X = np.matrix(np.reshape(params[:nm * n], (nm, n)))\n",
    "    theta = np.matrix(np.reshape(params[nm * n:], (nu, n)))\n",
    "    \n",
    "    # initializations\n",
    "    J = 0\n",
    "    X_grad = np.zeros(X.shape)  # (1682, 10)\n",
    "    theta_grad = np.zeros(theta.shape)\n",
    "    \n",
    "    error = np.multiply((X * theta.T) - Y, R)\n",
    "    \n",
    "    J = (1/2) * np.sum(np.power(error, 2))\n",
    "    \n",
    "    # add regularization to the cost\n",
    "    J += (lam/2) * np.sum(np.power(theta, 2))\n",
    "    J += (lam/2) * np.sum(np.power(X, 2))\n",
    "    \n",
    "    # gradients with regularization\n",
    "    X_grad = error * theta + (lam * X)\n",
    "    theta_grad = error.T * X + (lam * theta)\n",
    "    \n",
    "    grad = np.concatenate((np.ravel(X_grad), np.ravel(theta_grad)))\n",
    "    \n",
    "    return J, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(31.34405624427422,\n",
       " array([ -0.95596339,   6.97535514,  -0.10861109,   0.60308088,\n",
       "          2.77421145,   0.25839822,   0.12985616,   4.0898522 ,\n",
       "         -0.89247334,   0.29684395,   1.06300933,   0.66738144,\n",
       "          0.60252677,   4.90185327,  -0.19747928, -10.13985478,\n",
       "          2.10136256,  -6.76563628,  -2.29347024,   0.48244098,\n",
       "         -2.99791422,  -0.64787484,  -0.71820673,   1.27006666,\n",
       "          1.09289758,  -0.40784086,   0.49026541]))"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J, grad = cost(params, Y_sub, R_sub, n, 1.5)\n",
    "J, grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning movie recommendations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie_idx = {}\n",
    "f = open('machine-learning-ex8/ex8/movie_ids.txt')\n",
    "for line in f:\n",
    "    s = line.split(' ')\n",
    "    s[-1] = s[-1][:-1]\n",
    "    movie_idx[int(s[0]) - 1] = ' '.join(s[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Toy Story (1995)'"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "movie_idx[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rated Toy Story (1995) with 5 stars.\n",
      "Rated Shawshank Redemption, The (1994) with 5 stars.\n",
      "Rated Forrest Gump (1994) with 5 stars.\n"
     ]
    }
   ],
   "source": [
    "# create our own ratings for some movies\n",
    "ratings = np.zeros((1682, 1))\n",
    "\n",
    "ratings[0] = 5\n",
    "ratings[63] = 5\n",
    "ratings[68] = 5\n",
    "\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[0], str(int(ratings[0]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[63], str(int(ratings[63]))))\n",
    "print('Rated {0} with {1} stars.'.format(movie_idx[68], str(int(ratings[63]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1682, 944), (1682, 944))"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# add our own ratings vector to the existing data\n",
    "R = np.append(data['R'], ratings != 0, axis=1)\n",
    "Y = np.append(data['Y'], ratings, axis=1)\n",
    "\n",
    "Y.shape, R.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26260,)"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nm = Y.shape[0]\n",
    "nu = Y.shape[1]\n",
    "n = 10\n",
    "lam = 10\n",
    "\n",
    "params = np.random.random(size=((nm * n) + (nu * n),))\n",
    "\n",
    "params.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalize the ratings (Y - mu)\n",
    "Ymean = np.zeros((nm, 1))\n",
    "Ynorm = np.zeros((nm, nu))\n",
    "\n",
    "for i in range(nm):\n",
    "    idx = np.where(R[i,:] == 1)[0]\n",
    "    Ymean[i] = Y[i,idx].mean()\n",
    "    Ynorm[i,idx] = Y[i,idx] - Ymean[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: 38970.80995695868\n",
       "     jac: array([-0.21624453, -0.39736421, -0.22766823, ...,  0.01078863,\n",
       "       -0.01019731,  0.00483273])\n",
       " message: 'Maximum number of iterations has been exceeded.'\n",
       "    nfev: 156\n",
       "     nit: 100\n",
       "    njev: 156\n",
       "  status: 1\n",
       " success: False\n",
       "       x: array([ 0.58802966, -0.19909935, -0.2288541 , ...,  0.06169603,\n",
       "       -0.10768292,  0.00061393])"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# minimize the cost function\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "fmin = minimize(fun=cost, x0=params, args=(Ynorm, R, n, lam), method='CG', jac=True, options={'maxiter': 100})\n",
    "\n",
    "fmin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "# unfold the parameters\n",
    "X = np.matrix(np.reshape(fmin.x[:nm * n], (nm, n)))\n",
    "theta = np.matrix(np.reshape(fmin.x[nm * n:], (nu, n)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1682, 1)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predictions = X * theta.T\n",
    "myPredictions = predictions[:, -1] + Ymean\n",
    "myPredictions.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[5.00000118],\n",
       "        [4.99999985],\n",
       "        [4.99999929],\n",
       "        [4.99999922],\n",
       "        [4.99999894],\n",
       "        [4.99999848],\n",
       "        [4.99999842],\n",
       "        [4.99999803],\n",
       "        [4.99999733],\n",
       "        [4.99999597]])"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get top 10 recommendations\n",
    "sortedPreds = np.sort(myPredictions, axis=0)[::-1]\n",
    "sortedPreds[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1200],\n",
       "        [1499],\n",
       "        [1292],\n",
       "        [1121],\n",
       "        [1652],\n",
       "        [1535],\n",
       "        [ 813],\n",
       "        [1188],\n",
       "        [1466],\n",
       "        [1598]])"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# need to get the indeces instead of ratings\n",
    "idx = np.argsort(myPredictions, axis=0)[::-1]\n",
    "idx[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1200, matrix([[1200]]))"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "int(idx[0]), idx[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.) Marlene Dietrich: Shadow and Light (1996)  with predicted rating of 5.00\n",
      "\n",
      "2.) Santa with Muscles (1996) with predicted rating of 5.00\n",
      "\n",
      "3.) Star Kid (1997) with predicted rating of 5.00\n",
      "\n",
      "4.) They Made Me a Criminal (1939) with predicted rating of 5.00\n",
      "\n",
      "5.) Entertaining Angels: The Dorothy Day Story (1996) with predicted rating of 5.00\n",
      "\n",
      "6.) Aiqing wansui (1994) with predicted rating of 5.00\n",
      "\n",
      "7.) Great Day in Harlem, A (1994) with predicted rating of 5.00\n",
      "\n",
      "8.) Prefontaine (1997) with predicted rating of 5.00\n",
      "\n",
      "9.) Saint of Fort Washington, The (1993) with predicted rating of 5.00\n",
      "\n",
      "10.) Someone Else's America (1995) with predicted rating of 5.00\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# get the movies names from the indeces\n",
    "for i in range(10):\n",
    "    j = int(idx[i])\n",
    "    print('%d.) %s with predicted rating of %0.2f\\n'% (i+1, movie_idx[j], myPredictions[j]))"
   ]
  },
  {
   "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
}