Assignment -  Week 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Expectation-maximization algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this assignment, we will derive and implement formulas for Gaussian Mixture Model — one of the most commonly used methods for performing soft clustering of the data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Installation\n",
    "\n",
    "We will need ```numpy```, ```scikit-learn```, ```matplotlib``` libraries for this assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import slogdet, det, solve\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from sklearn.datasets import load_digits\n",
    "from grader import Grader\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grading\n",
    "We will create a grader instance below and use it to collect your answers. Note that these outputs will be stored locally inside grader and will be uploaded to the platform only after running submitting function in the last part of this assignment. If you want to make a partial submission, you can run that cell anytime you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grader = Grader()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing EM for GMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For debugging we will use samples from gaussian mixture model with unknown mean, variance and priors. We also added inital values of parameters for grading purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX1sXOd15p8z/JIlK2O5Q9nRKFyrUCJFluikpLttCXRrpbvNsnG8BoIgEVtUaUBBxNZNuwWKfgDbvxYosEG3QVqIMNtEWZhqV3CzSLyYbZJWKhY7QANxXGlIWaKWtVFGpCTO2PI4kawhKZ79g7yTmTv3c+YO770zzw8wLM7cufflcOa5533ec84rqgpCCCHtQyLsARBCCAkWCjshhLQZFHZCCGkzKOyEENJmUNgJIaTNoLATQkibQWEnhJA2g8JOCCFtBoWdEELajO4wLppKpfSpp54K49KEEBJbcrlcUVX73Y4LRdifeuopzMzMhHFpQgiJLSLyL16OoxVDCCFtBoWdEELaDAo7IYS0GRR2QghpMyjshBDSZoSSFUPiSalUQjabxdLSEtLpNEZGRpBMJsMeFiHEBIWdeKJUKmFychKrq6vY2NjA7du3MTs7i9OnT1PcCYkYgVgxIvLbInJVROZE5K9EZEcQ5yXRIZvNVkQdADY2NrC6uopsNgtgU/gzmQympqaQyWRQKpXCHC4hHU3TEbuIpAH8JoAjqvq+iJwH8DkAZ5s9N4kOS0tLFVE32NjYwNLSEqN5QiJGUIun3QAeEZFuADsBLAd0XhIR0uk0Eonaj0sikUA6nXaN5g0Y1ROyPTQdsavqkoh8GcAigPcBfFdVv2s+TkROATgFAAMDA81elmwzIyMjmJ2drQh4IpFAb28vRkZGcP78edto3sAuqh8bG8OlS5ewsLAAEcHBgwfx3HPPMdInpAmajthFZA+AFwAcALAPwC4R+RXzcar6sqoOq+pwf79rDxsSMZLJJE6fPo2hoSHs27cPQ0NDFavFKZo3sIvqv/71ryOfz+P+/fu4d+8erly5gjNnzjCaJ6QJgsiK+UUAb6lqAQBE5JsAfg7AKwGcm0SIZDKJ0dHRusedonkDO4/einK5jGw2a3ktQog7QQj7IoCfEZGd2LRiPgGArRs7CCOad8pxT6fTuH37tq2Ym6m2cQgh/gjCY/++iLwK4HUA6wD+CcDLzZ6XxIvqaN4oZFpcXISqQkTw5JNPoqenB2tra5WoXkTw8OFDy/NV2ziEEH8EUqCkqn8E4I+COBeJN8YiablchqpWHi8UCujp6cHRo0dRLBaRTqcxODiIV155BeVyueYcfX19NTYOIcQfrDwlgWIsklaLOrDpp6+traGvrw/j4+OVxycmJnDhwgVmxRASIBR2EihWi6QGGxsbWFxcrHksmUzixRdf3I6hEdIxUNiJJ0qlkqfI2m2R1BzJt2qsbFZGOhnZji+ameHhYeWep9HFLIxOXvjExESNaBoe+4MHDyzPnUqlcODAAVfRbVSczYVQRuol2xuQdkBEcqo67HochZ1UYyWMTtkrzz77bF2+ealUwrlz57CysmL5GhGBqtqKbjPinMlkkMvlamYMiUQCQ0NDzIsnscersHOjDVKDVYWonagD1vnmyWQSJ06cwI4dOyAidc8bwcTGxkalGMltDFa9Z+zG49begJB2h8JOanBa/LTCLt/cKFpyax+hqjULqqVSCVevXm1YnL20NyCk3aGwkxqshNGJ1dVV274uyWQS3d3u6/NGBG9YMPfv3687xqs4j4yMoLe3t/I7WLU3IKTdYVYMqcGq74uq2mazzM7O4vr16zh06FCl8Kh6odNLKwHDrjEsGCu8irOX9gaEtDsUdlKDlTCWy2XMzc1ZirPhk+fzeQCo22TDfKMwk0gkKm2c7WygXbt2YXx83LM4W7U3oMiTToJWDHHl2WefrbE3nDAvdFa3+33iiSfQ1dVla5PY+eNHjhxpSIwNayeXy2F5eRm5XA6Tk5NsCUzaHkbspAanDTHy+Tzm5ubw/vvvO57DvNDpNYL20v7XD07ZNUx9JO0MhZ3UYCeG+XweIyMjFcvFCaeFTrue7sZzTv64X1uFqY+kU6GwkxrsxNBISVxbW6t7jWGfbGxsQEQgIlhcXEQmk6mIb6lUwsWLF7GwsABVxcGDB3H8+PE6YbYT/kY2zLZauGXqI+kEKOykhnQ6jVu3btVlwRSLRQDWux719/djYGAAi4uLKBaLUFXcuXMHhUKhYuOYWxLk83nMz8/XtCRwisgbsVWCtnYIiQtcPCU1jIyMWC6SGimPVoubAwMDGB0dxcDAAFS1Tnxfe+21uj4zAGqqTt0WOhuxVZz2aSWknWHETmpIJpNIpVK4c+dOzeOGzdLb22sbAduJ7927d22vZwizVUT+4MEDnDt3DidOnLDNh19fX0epVLIVaydPn5B2hRE7qWNgYMA2MneKgO3SFffs2WN7rXQ6jVKphDfeeMPS5llZWcHk5CQGBwctUy4LhQJTGAkxQWEndTiV5RsR8Pj4OEZHR2siZfPrgE0L5wMf+IClvdPX14fBwUFMTk7i3r17tuMxsnJOnz6NVCpV85yqem4QRuJLqVRCJpPB1NQUMpkMb+Qu0IohdfgpyzcveI6NjeHSpUuYnZ2t+PILCwt1r/voRz+KX/qlX3JsI2BgeOl2vWeYwtjeNJIR1elQ2CNK2KXwXrxpuy/coUOHKj3XrUgkEnj00UeRTCY9dZOsTlG089rduki6Efb7TexhoZl/KOwRJC4Rit0XbmFhwVGsqyNstyZh5gVao0jKnGVz/fp1x0VUJ+LyfrcbXm+mLDTzD4U9gkQlQrH64hnjW1pawrvvvmv5hfMTgVvlmvf09ODw4cMoFAp1X/hkMolDhw7VVcCura01/P5E5f3uJPzcTFlo5h8KewQJOkJxE2iraOnmzZs4e/ZsZfek27dvV8R0bW3NMcLeuXOnbT8ZcwTeSJtdo1iqmmbeH0aE24+fmykLzfxDYY8gQUYoVpGRWaDN0VKpVKoRdeDH7XmdML5wdptr9PT04GMf+1idcPvNNQ86gmNEuP34uZmyx75/KOwRJMgIxSoyMgu0OVrKZrOO+5ya2bVrF5LJZOULl81mUSgU6r64R44csezymEqlICKW1osVbu+P34VQRoTbj9+bKQvN/EFhjyBBRihe9zCtjpb8WBBGv3RD0M+fP4/+/n50d3fXpTEaC5wAamYRy8vLlWO8LFxavT+Dg4PIZrM1/WrsvFsr4WdEuL3wZtpaxC4lrZUMDw/rzMzMtl+3nbGLUjOZDHK5nKcFzaGhIYyOjiKTyWBmZqYuXbGrqwvd3d0VC8f4Mo6NjWF6errmSwpYNwwbHBxEX1+f45iqx+L1d6++UTidz3ys8TswA2b7YYqpf0Qkp6rDbscxYm8DnDIM3LamM+jp6alJKZydnUW5XK6Ie1dXF06ePAkAeO2113D37l3s2bMHzz//PPL5fJ3dY8fs7Cz27t3rOR3SC2a7yel8zICJDrRXWgdbCrQBTmKVTCYxNjZW2TDajkQigWw2W8kFP336NIaHh7Fv3z48++yzeOmll7B7925MT0+jWCxibW0NxWIR09PTeOuttzzZPYB9l0jzWPwsXLrZTYlEAv39/chkMrh8+TIzYLYZtgPYfhixxxhjKusmVvl83rYK1OD9999HLper8aPNC51vvPEGHjx4UHON1dVV/OhHP/I1bnOXyGqMPPZyuYypqSlPU3SnIifjfNevX7dN02QGTOtg8Vc4BBKxi8hjIvKqiFwXkWsi8rNBnJfYU92/3G5XI0Os/CygmhtqVV/HqlHXxsYG1tfXfY1dVTE2NlbpEjk4OIhnnnkG+/btw7FjxwAAc3NznjegNjcfExF0dXXhiSeewLFjx/Doo4+iXC7bijoX7VqH02yStI6gIvavAPhbVf2MiPQC2BnQedueRheQnHxls1i5le1XY7Yl3PzrRCKBxx9/HIVCwXVWYFAoFDA9PV0XtZVKJZw7d64mHdOLB26XRQRsZt9UzzKqscurJ8HB4q9waFrYRSQJ4OcBnAQAVV0F4NyujwBobppqF4VbiZXXBVSD6s0rnKJ94wby/PPPY3p6umax1YnqVrvVdo+dCHsRAquFuEwmY9s5MpFIVN4nc9pkPp9npkZAsPgrHIKI2A8AKAD4uog8AyAH4EuqWjNvF5FTAE4Bmxs5kOYyNOy+MB/72MfqXmuOaN99913cv3/f9twrKys4c+YMDh8+jHfffdfymJ07d+Lpp5+uCN/p06dx7tw5rKys1ByXSCSwY8eOuuvZzQysaFQInG5Kqoof/ehHOHPmTMV7v3XrFi5duoREIkE/OCCYrx4OQXjs3QB+CsAZVf04gHsAfs98kKq+rKrDqjrcbIvVdqGZaarTZhhWVG+Q8fTTTztmpQCb+5FeuXLF8gbQ1dWFz3/+85UbSCaTwfnz5y0XUY0t7syYxdpJhKtTMf1gtaOTgari2rVrNd67MdugHxwc3Hc2HIKI2G8CuKmq39/6+VVYCDupp5lpajPVqX6tGTOqiunp6brCJDusMl/MN6F0Ol1TgVrNT/7kTzYkBM3+ngD94CBgvvr207Swq+ptEfmBiBxS1XkAnwDwRvNDa3+anaY2+oWpvilcunTJ9+uNKPwb3/iG74yYHTt24NixY3U3ocHBQduxNNpr3cjhNwqqVNX3eOkHNw4rS8MjqKyYlwBMb2XEvAngCwGdt62JQte67u5u32Jn0MjrHjx4YNkl0txfvRpVbagytFQqeZpRGBi7PhkeO/3gxmH+ergEIuyqehmAa/8CUk8Y01S33irApq/98OHDhi0MJ/L5PObn5zExMVH5krvZHY3YIV5SNc2benR6VkxQUTZbN4QLK087EDvBExH09/fj+eefx+7du2va6s7PzztusOGXcrmMCxcu4MUXXwTgXj3qZodYCZLdgqy5zbBZuPbv39/EbxZfgoyy/SQG0LIJHgp7B2IneB/84AcxPj5e+bk6xxwAFhYWsLq62rB1Y2ZhYaFyfrsceKP9gJMd4rSpttXidHVfePJjgoyyvSYG0LJpDRT2DsRPNo5Vm1snRASPPfYY7t696zqO1dVV3Lx5s+KDG8IuInj88cfR3d2NgYEB1wjOTpBUtaYnDT1zZ4KsEvWaGEDLpjVQ2NsYuymun2wcqy+ecbxVfvrg4CAKhYInYV9fX6/bgg/YXCxNp9MVm8bt97ETpGKxGPridJwIskrUa2IAWw60Bgp7m+I2xfUqeIuLi5ZfPLut84xFSK+9aezOMzs7i+PHj1fGZLW5tvH7OAkSc6i9Y5X3b7R/qE439eqJe3nv2XKgNbAfe4xx6nPt1lWvuhJ1dHTU8otZKpVQLBbrHk8kEtizZ0+dLWN8IUdGRtDT01PznIi49oSvxkhxNMZhtbm28fv4rcIl1hg3/KNHj1b+VqqK2dnZSofN6m6fXrtvOsG/XWtgxB5T3CLyIKa42WzWMuoWkUrjL6/+dSNbMFbvemQV2Ru/TxTqAdqFZDKJvr6+Sk4/UB8UBOmJ82/XGijsMcVt0clPVoLdl2ppaclSkFOpFPbv32/7hcxkMpY94v1g7idvh3EMLZfgcAsKgvbE+bcLHgp7THH78nlZIHWL+u1uDkZ3TrsvpNeNPewwpzim02ncunXL8iZTLpcbajdA7HELCuiJRx8Ke0wwR9b9/f2OXzDzFLe/vx+qivPnz1eia7eo30/2TPX4rFoGeCWVSuHAgQOW/eStct3n5uZw48YN5j0HiNvf3eq5wcFBZDIZX3YKC5NahzTifTbL8PCwzszMbPt144pVLrmxOGlUgxpfMCuBs3p9b28vkskk7ty5U3e9nTt34tSpU0gmk56+fMb5vW60Ycfhw4fx3nvvWV6net9V8xZ9iUQCQ0NDnM4HiPF+Ly4uQlUhIpWaAgA1AcPq6iquX79e+ds7fRarz2/1meQN2hkRyamqa/sWRuwxwCqyXltbw9GjR9HX1+ca8TgV8BgNr6q5f/8+Jicn6za1dhqfF1Hv6enBRz/6UczNzVlaNdevXwcALC8v4/XXX8fJkycr5f3GOJaWluqEnXnPwVF9I+/v78e7775bCR4KhULFqhsdHbXd9crLgioLk1oL0x1jgFMBjlvKotPrNzY2bFMQ/WwwYbfIWo2xu9Px48fR19fnes6HDx/i7NmzdWl0Vptn0OMNBnMq45UrV2o2IjFnxzjteuV2s2VhUmthxB4Dmi3isHv93bt3bRc5zV8yJ0vGaZMMAxHBW2+9hddff922KMnMw4cPceHChZpZyeDgILdaaxFu3TCB2s+F2364Tp9PFia1Fgp7DGh2Qw6r14sINjY2bCPt6i+ZW/aM0yYZBqpqWezkxuzsbGWsxnXHxsY6urVuq/CSzVT9uXCqMHb7fHIv1NZCYY8BzRZxWL1+cXHRcuEUqK/+c/ND8/l8TUGLFU6C8cgjj9RM+atR1bpCmXw+Tx+2BfT39zvOvMyfC6sWBCKCY8eO1bSDsIKFSa2Fwh4Tmi3iML8+k8mgUCjUiemuXbtw5MiRukIlJz/Ui8fuxO7du/HpT38ar776asWmqS5pt7ouU+WCx+5v+BM/8RPo6+ure5+DCDh4g24NFPYOoloMU6kUenp66tIlx8fH676YdlPu9fV1lEolxyl5V1eXo+UDbDYO+9a3voWTJ0/WWCyrq6uYnZ2t82GTySS++tWvVm4Ct27dYg/vALCzyvr6+mr69FdDcY4mFPYOwcon7+npwbFjx1AoFGoKmFKpFESkZrs4qwKhYrGIyclJjI2NWU7JBwcHMTw8XNNTxsqyMToImi2WUqmE+fn5uvz9Gzdu1CzAqirK5TJT5ZqEC5rtA4W9Q7DLhe/t7cVnP/vZGtGv9lmrFyxfe+01rKysVJ6r9rydpuRe/H2rVDerqf7q6iquXLlS9/upKlPlmoQLmu0Dhb1DcPLJndLcqsW7u7v+41LdYdEuWvbi79tFhubXTk1N2f6OjCybgwua7QOFvUNIpVKWGQ+pVMo1zc0Q76Cm6s1EhnZ+fldXFyPLAKBn3h5Q2DsEuza6a2trrjseVW+g4VeQ7bJXGo0MrVLsurq6cPLkSUaWhGzBJmAdwpe//OW6HivAZnrj+Ph4jcdejbk5k9emYEYDqWKxWMmKCarRE1MdSafCJmARJExBsruBq2pdBG3OijHnLjtN1c3ZN9UE1eiJdgEhzlDYtwm3svxWc/DgQeTzecvHgeDE0q3fSJCNnhi5Nwbft/aHwr5NhN2m9Pjx45ifn6/ZBKOvrw/Hjx8P9DpuC7GJRAKpVKqyKYORP18sFn2JTNg3yrjC960zoLBvE2G3KU0mk5iYmGh5pOa0EGsUGM3Pz1cqXq1y5r2ITNg3yjhgFZnzfesMKOzbRBSq+vzYLW7Tdbvn7TpJplIpDAwMWLYJMPAjMmHfKKOOXWSeTCb5vnUAFPZtIk5VfU7TdQC4ePEi8vl8ZUHWHGk7pTJOTU157vfthN3MoFQqIZPJdLxv7GfXLLYNaD8o7NtEnKr67ETh4sWLmJ+fd90KzWlm4KU1rBeRscpnB4B79+4hl8t1vG9sN6MREfT29sYiwCCNE5iwi0gXgBkAS6r6qaDO207EJU3PThQWFhYa3goN2IymjX1NrfAjMtU3yqtXr+L+/fs1Y+k039hsjaVSKUvrz9iQOg4BBmmcICP2LwG4BuADAZ6ThIDdeoCqOtooqVTK8bzZbNayAtau37cb1RtcVws7EB3fOIjUQi/rHVadO63aMhuv7ZQbXqcSiLCLyH4AvwzgvwD4T0Gck4SH3XrAwYMHMTc3Zyvu8/PzKJVKvjbVBpz7fXshCgvTVgSRWujlHNlstqalsjFjGRwcRG9vLyPzDiThfogn/hTA7wKwDedE5JSIzIjITKFQCOiypBUYNsfQ0BD27duHoaEhnD59GsePH0dvby8SCeuPzdraWmUHeyvS6XTda4MQ4JGRkZpxRcU3dkotDPIci4uLlj3ub9++jdHRUYyPj1fWPkhn0HTELiKfArCiqjkR+QW741T1ZQAvA5u9Ypq9LmktdtN1w9e+fPlyna3iZn+0KjMoqgvTQaRkejmHXbuIt99+mxlCHUoQVswIgE+LyCiAHQA+ICKvqOqvBHBuEjGqBT+Xy/myPxoRYK8edRR94yAsIi/nMPaHNbO+vs4MoQ6laStGVX9fVfer6lMAPgfgAkW9/WnU/jAE2Is9YPjLuVwOy8vLyOVymJycRKlUCvR3aRVBWERezjEwMGBrjzVi/5D4wzx20hDbYX/Evfw9iPfI7hwAavrtVGfAmIlKhhDZPgIVdlX9BwD/EOQ5iTthdetrtf3RDm0DgniPzOewS288evQoFhYW6lI/o5AhRLYXRuwxp9GUuji0bo1qGmPYWM1kyuUyRASnTp2q+TxEJUOIbC9BpTuSkGgkpS4u3nVU0xjDwuiDc/nyZUvLxei3b5WqGrWbNmktjNhjTiN2RVy8a8NfvnjxIhYWFqCqlY1BOg2nnakMVLXyN4zS35FsPxT2mNOIXRGmd92IBTQ/P18RtLm5Ody4caPjolC3nakM4rT+QFoHrZiY04hd0aoKUDcasYCCqN5sB9x2pgK4/kB+DIU95tiV/ztFs2F5142IdDtkxgSB1c24mk5ffyC10IppA/ym1IVVgt+ISKdSKcv+7W6dJNsNq3YMPT09OHz4MAqFQmQzm0g4UNg7lDBK8BtZD7Arl7d7vF2Jaj8cEk0o7GTbaKQJmF0n0E7oEGq10MxsF+IFCjvZNhqJOju1SMlL4VkcisxIOFDYSUM0Kip+LaA4bQIeJG61BkFs4kHaFwo78U1QouLl5tCp3rLbQnNcisxIOFDYiW+CEBU/N4co9lpvNW4WFNNAiRPMYye+CUJUOq3wyOjzMjU1hUwm49qXx63WIKwiMxIPGLETT1TbJuvr6xCRmi3Z/IpK3CLOZhYqvc5OzNcYGxtDPp+3vGanrj0Qb1DYiStmYUokElDVirg3IipxynZpdk3Bi3Xl9xqduvZAvEFhJ65YCVMikUAqlUJ3d7cvUTGi0sXFRYhIUzeH7aLZNQUvs5NGrtGJaw/EGxT2DsXOWrB63E6Yuru7MT4+7uua5sjfuEEMDAxENuJs1jbyMjuJmzVFog2FPWaUSqW6/uTHjx/3HS0Xi0Woas20f2xsDNPT0xXhXV5exszMDB5//HEkEommbRO7yH9gYCDSkWeztpEXPzxO1hSJPhT2GFEqlXDmzBmUy+XKY/l8HvPz85iYmHDdCs9uowZj2v/aa6/VPa+qePvttwGgIu6N2iZxiUrNs5bBwcGmFiq9+OFcDCVBQmGPEdlstkbUDcrlsqvf67ZRw8bGBu7evWv7vIg05KlXE4eo1G4R0ylDxQtufjgXQ0mQUNhjhFNk6xb1um3UkEgksGfPHhSLRcvjVNW3p24mDlGp3SJmPp9vuV3ExVASFCxQihFOka1b1Ou0UYMhsM8//zx6e3ttj2k2sm5kU5DtJi52ESFOMGKPESMjI8jn83V2TF9fn2vUa46WRcQyI8XYPDqfz1cKkIKMrKMelcbBLiLEDamuHtwuhoeHdWZmZtuv2w4EkRXjxcPt1JawVimZvb29kZtZkM5ERHKqOux6HIW9c4mSeHMshLhDYSeORCkyjdJYCIkyXoWdHnsHYW7kFZV+3uwtTkiwUNg7BKcCJQOr7I/tsCWYiUJIsFDYOwS3AiWgPvujVcU6ZpiJQkiwUNhDZDsX6bwUKJlTGq0skgcPHuBrX/taJRUyiL0241C4REicaFrYReRDAP47gCcAKICXVfUrzZ633XHqvw0gcMG3iopFBP39/bZtAuxuBtUL7o364X42lSCE+COIiH0dwO+o6usishtATkS+p6pvBHDutsVuwfDixYuYn58PfPd5u6j4xIkTtue1uhlY4dcPD2ozbEKINU0Lu6reAnBr698/FJFrANIAKOwO2C0YLiwstCRDpJEmU+abgRN+/HA7i+fcuXOONxo7rCwt4zqcAZBOJFCPXUSeAvBxAN+3eO4UgFMAMDAwEORlY4ndgqHRI72aoDJE/JbzGzeDCxcuYHZ2FnY1D11dXb78cDuLZ2VlBZOTk74id6voP5/PAwDW1tY4IyAdSWBNwETkUQB/A+C3VPU98/Oq+rKqDqvqcH9/f1CXjS12u9AfPHiw5bvPl0olZDIZTE1NIZPJoFQq2R6bTCbR19dnK+oA8JnPfMaXYDo1JDNmJ16xiv7L5TLK5bLlrIeQTiCQiF1EerAp6tOq+s0gztnu2FkjAHDjxo2WZYg04m+7zRbefPNNHD582PMYDIvnwYMHdc/5nZ24Zfs0el5C4kwQWTEC4C8BXFPVP2l+SJ2DnTXSqg0XSqUSzp07VyOoXjz8dDqN5eVl2/P6FUzjpnbu3DmsrKzUPOd3duJ1gZd58aSTCMKKGQHwqwCOi8jlrf9YB94EhuCPj49jdHQ0MFGfnJysE1LAPZodGRlBX1+f5XONCmYymcSJEyewY8eOOjvKz+zEytLq6+tDX19fU+clJM4EkRXzfwFIAGMhVQRdvGR40Va4iXMymcTExAS+853v4Nq1azWva0Ywg9gOzsnSYlYM6VTY3TGCtKLb4dTUlKOd8sUvfhH79+/3NDYKJiHhwO6OMaYV3Q6dvGgRQT6f9yTsdusCrRR83kwI8QeFPYK0otuhUyaKqjZ17lZWkt68eRNnz57Fw4cPAQTTm4aQdoebWUcQqzzvZrM6DC967969dc81e26nGUYzlEqlGlEP8tyEtDMU9ghiV7zUbFZHUJkoZlrVTz2bzdaIepDnJqSdoRUTQYLIFqmm1Z0UW9VP3Um8mZNOiD0U9ojit6+LHdvRSTGofurmG1AqlcKtW7fq2hl47U3DRVfSqTDdsc3JZDLI5XJ10fTQ0BBGR0cDE79mz2OV4tnT0wNgs39M9ec0lUrhwIEDjtfgBtmkHWG6IwHg7H8HGc03O8OwWoBdW1vDsWPH0Nvbi8XFRRSLRagqisUi3nnnHcexcoNs0slw8bTNseqkafjfrcpmaQS7G1ChUMDo6CgGBgZqWhq7jZUbZJNOhsIk3NprAAAOR0lEQVQeE/y02q1+zfXr1+se7+npwcjISKTEzy3F0+9YW5EySkhcoBUTAxrdHzWbzWJtba3ufIcOHUIymWxZNosx5mZ2azIvwPodKzfIJp0MhT3CGOJ49epVy1a7bvuj2vUqLxaLAFonfo14924pnn7HGnTKKCFxgsIeUcziaMbL/qheotyPfOQjWFhYgKpi586d6O7uRjabbUoEG124dFqAbUSog0oZJSRuUNgjilkczXjZH9UpyrW6cbz//vsAgEKh0FSue6u8e7NQG+sOjMgJqYXCHlGctnyr3h91bm7ONiJ3inIzmYzjbMBvamC1p76+vg4Rqck9b8W+ra0uvCIkrlDYI4pdm91du3bhyJEjnvdHtbMj3PYK3djYwOLioqeI2KoYyFz4ZmTiBEWjdg+rUUknQGGPKHY2yvj4eI0QNbpA6LZXqIigWCyiUCi4RsRWImtmbW0NP/zhD+te26jQNmL3MMonnQLz2COKYaMMDQ1h3759GBoashSgRvdHNXeQrCaRSNR5+E4FQW7Rv/H6s2fP1uTfG0I7MzOD5eVlXLp0CV/96ldx8+ZN1/E3kqcepYIsQloJI/YI08qsDrP/nkqlICIoFApIp9NYXFzEnTt3al6zsbGBy5cvAwAGBwcrHSLX19eRSCRcxf3hw4c4d+4cTpw4gWQyiWw2i3K5XGPbPHz4EGfPnsVLL70EALbjGxwc9J2qGaWCLEJaCYU9JKLg9TrdODKZTMWGqWZtbQ25XA6XLl2qLJAa//ci7isrKzhz5gwmJiawtLRU58UDm+J+4cKFmvWD6v1aDQvFb/vhVhZkERIlKOwhEAev1+zxV2P8bIiyIe6pVArd3d1IpVJ12TrVlMtlXLhwAel02naDbXOOvvn6q6uryOfzvmY0rEYlnQI99hCIg9db7fEb7XOdUFV0d3djfHwcL774Ir7whS+gq6vL9viFhQWMjIxYHpNIJCAirlk7fi0Ur+sWhMQdCnsIxM3r7e52n9iJCNbX1ytNynbv3o2XXnrJ9rWrq6sAgJMnT9aIe3WOvtXCbvVxjVgojS42ExInKOwhEIfOg4ZdlMvlKhWpZkSk8n+jT/ry8jJyuRwmJycBAE8//bTla9fX1zE5OVm5ATz77LM1UfRzzz3nmLVDC4UQe+ixh0AcvF63lgYHDx7Ee++9h7t376KrqwvlctnSWnruuedw/fp1lMvlunNUFxRZeeVOWTssLCLEHgp7CMSh86Bbbvo///M/A9j01q1aAxvWUjKZxMTEBKampnDv3j3LY+xgEy9CGoPCHhJRFy23ylS3vXLNPWuOHDliufdqs/ZTFNJGCYkaFHZiiV26o+GrOwm7sRH16uoqpqamGi4ociMOaaOEhIG4RV6tYHh4WGdmZrb9unEhKlGoMY7FxUWsr6/j/v37EBE88sgjePvtty1f093djWQyiXfeeaci/oaI+y0ociOTyVjOAoaGhiI9GyKkUUQkp6rDbscxYo8YUYpCDbvIPCa7LBmjv4xZ9Dc2NlAul30XFLkRt7RRQraLQIRdRD4J4CsAugD8har+cRDn7UQabUfrhJcZgNMx5jEZkbi557rTYquqYnFxsaGx2cEWAYRY07Swi0gXgD8H8G8B3ARwSUS+rapvNHvuTiTIKLRUKuHChQuYnZ2tCLDVDMBtlmCXIbN37148+eSTyOfzroupQL0v3+zsJOi00ahYYIQ0SxAFSj8NYEFV31TVVQB/DeCFAM7bkQRVvGSIpll0rdoXuLU4sBvTwMAAent7KwuqbpiPa7a1QpAtAqoLsqqLrKrbDBMSF4KwYtIAflD1800A/9p8kIicAnAKAAYGBgK4bHsSVBRqiKYV5hmA2yzBaUznz5937egI/PhGUI2f2YldNB1U2mgrLDBCwmLbFk9V9WUALwObWTHbdd24EVTxktueqdUzADev2mlM/f39th0aq88lIpWt9ozXevXIb968ibNnz+Lhw4cAgFu3bjW1oGx1k+BCLGknghD2JQAfqvp5/9ZjpEGCiEKdCozMMwAvswS7MTl563v37q1ssbexsYE7d+6gUChURNnLdUulUo2oG9csl8sNRdN2vv5HPvIRLsSStiEIYb8E4MMicgCbgv45ACcCOC9B4wt6VgVGIoLBwUE899xzNedoZpZQLBZtn+vu7kY6nUahUKjcAMwWh9t1s9lsjagbqGpD0bSd5SIi6O3tjXT/HkK80rSwq+q6iPwGgO9gM93xa6p6temRdSjVQp5KpTA/P4+1tTXfWSN+xbrRWYLTZhnpdNrV4nC7rpN4NxJN242nUChEvn8PIV4JxGNX1QyATBDn6mTMNoFZMP0u6G1HP5qRkRHk8/m67o19fX0YGRlBNpttyuKws5S6uroaiqadfP2o9+8hxCvsxx4h3FrlAtFb0DO6Nw4ODmLnzp3YtWsXnnnmGUxMTCCZTGJkZKSmr7pfi8P8emBT1E+ePNlQNN3oeEqlEjKZTGUjEaZBkijDXjERYmpqylOGSdx6oTRb+BN04ZDf85lnUsbNgM3GyHbDXjExxK1VblwX9Jq1OIK2SPyejznuJG5Q2COEVfpfT08PDh8+zJ2DQoQ57iRuUNgjRBx2VupE2GyMxA0Ke8RgZkb0iMMetYRUQ2EnFdjd0BrOpEjcYFYMAcDMD0LiALNiiC+inPnBmQQh/qCwEwDRzfyI0laBhMQFVp4SAI1v8NHqisxmN+MgpBNhxE4ANJb5sR3RdFRnEoREGQo7AdBY5kcjvrxfv5w55IT4h8JOKvjNofcbTTcS4TOHnBD/0GMnDWPlywM/9t3NfnsjfnmQG1YT0ikwYicNY7VLEwDcu3cPuVyuLhpv1C9nNS4h/mDEThqmOpreuXNnzXNW0XijmTdhwj7sJI5Q2ElTGNH0Y489VvecORpvdtON7cZYE8jlclheXkYul8Pk5CTFnUQeCjsJBC/ReNz8cubQk7hCj50EgtfslTj55cyhJ3GFwk4CoR07IDKHnsQVCjsJjDhF415gDj2JKxR2Qmxox1kI6Qwo7CQSRLU1b7vNQkhnQGEnocPWvIQEC9MdSegwrZCQYKGwk9BhWiEhwUIrhgQOW/MSEi4UdhIobM1LSPjQiiGBwta8hIQPI3YSKGzNS0j4NBWxi8h/FZHrIpIXkf8pIvUt/khHEcfWvIS0G81aMd8DcFRVBwHcAPD7zQ+JxJm4teYlpB1pyopR1e9W/fiPAD7T3HBI3GEZPiHhE6TH/usA/ofdkyJyCsApABgYGAjwsiRq0C8nJFxchV1E/g7AkxZP/aGqfmvrmD8EsA5g2u48qvoygJcBYHh4WBsaLSGEEFdchV1Vf9HpeRE5CeBTAD6hqhRsQggJmaasGBH5JIDfBfBvVPV+MEMihBDSDM1mxfwZgN0Avicil0VkMoAxEUIIaYJms2IOBjUQQgghwSBh2OIiUgDwLwGcKgWgGMB5WgXH1xwcX3NwfM0RxfH9K1XtdzsoFGEPChGZUdXhsMdhB8fXHBxfc3B8zRH18TnBJmCEENJmUNgJIaTNiLuwvxz2AFzg+JqD42sOjq85oj4+W2LtsRNCCKkn7hE7IYQQE20j7CLyOyKiIpIKeyzVRLVnvYh8UkTmRWRBRH4v7PFUIyIfEpGLIvKGiFwVkS+FPSYrRKRLRP5JRP5X2GMxIyKPicirW5+9ayLys2GPqRoR+e2tv+2ciPyViOwIeTxfE5EVEZmreuxxEfmeiPy/rf/vCXOMfmgLYReRDwH4dwAWwx6LBZHrWS8iXQD+HMC/B3AEwOdF5Ei4o6phHcDvqOoRAD8D4D9GbHwGXwJwLexB2PAVAH+rqocBPIMIjVNE0gB+E8Cwqh4F0AXgc+GOCmcBfNL02O8B+HtV/TCAv9/6ORa0hbAD+G/Y7FkTuQUDVf2uqq5v/fiPAPaHOZ4tfhrAgqq+qaqrAP4awAshj6mCqt5S1de3/v1DbIpSpLZgEpH9AH4ZwF+EPRYzIpIE8PMA/hIAVHVVVd8Nd1R1dAN4RES6AewEsBzmYFT1/wB4x/TwCwC+sfXvbwD4D9s6qCaIvbCLyAsAllT1Sthj8cCvA/jfYQ8CmyL5g6qfbyJiwmkgIk8B+DiA74c7kjr+FJvBxIbbgSFwAEABwNe3rKK/EJFdYQ/KQFWXAHwZmzPsWwBKpk17osITqnpr69+3ATwR5mD8EAthF5G/2/LizP+9AOAPAPznCI/POMa1Zz2pRUQeBfA3AH5LVd8LezwGIvIpACuqmgt7LDZ0A/gpAGdU9eMA7iFCNsKWV/0CNm9A+wDsEpFfCXdUzmy1JI+cI2BHkDsotQy7nvAicgybH44rIgJs2hyvi8hPq+rtsMdnEMGe9UsAPlT18/6txyKDiPRgU9SnVfWbYY/HxAiAT4vIKIAdAD4gIq+oalTE6SaAm6pqzHJeRYSEHcAvAnhLVQsAICLfBPBzAF4JdVT13BGRD6rqLRH5IICVsAfklVhE7Hao6qyq7lXVp1T1KWx+oH9qO0Xdjaqe9Z+OUM/6SwA+LCIHRKQXmwtX3w55TBVk8y79lwCuqeqfhD0eM6r6+6q6f+sz9zkAFyIk6tj6/P9ARA5tPfQJAG+EOCQziwB+RkR2bv2tP4EILe5W8W0Av7b1718D8K0Qx+KLWETsMefPAPRhs2c9APyjqp4Oc0Cqui4ivwHgO9jMSPiaql4Nc0wmRgD8KoBZEbm89dgfqGomxDHFjZcATG/duN8E8IWQx1NBVb8vIq8CeB2b9uQ/IeQqTxH5KwC/ACAlIjcB/BGAPwZwXkS+iM1utJ8Nb4T+YOUpIYS0GbG2YgghhNRDYSeEkDaDwk4IIW0GhZ0QQtoMCjshhLQZFHZCCGkzKOyEENJmUNgJIaTN+P+glUVW3DxoyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = np.load('samples.npz')\n",
    "X = samples['data']\n",
    "pi0 = samples['pi0']\n",
    "mu0 = samples['mu0']\n",
    "sigma0 = samples['sigma0']\n",
    "plt.scatter(X[:, 0], X[:, 1], c='grey', s=30)\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reminder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, that EM algorithm is a coordinate descent optimization of variational lower bound $\\mathcal{L}(\\theta, q) = \\int q(T) \\log\\frac{P(X, T|\\theta)}{q(T)}dT\\to \\max$.\n",
    "\n",
    "<b>E-step</b>:<br>\n",
    "$\\mathcal{L}(\\theta, q) \\to \\max\\limits_{q} \\Leftrightarrow \\mathcal{KL} [q(T) \\,\\|\\, p(T|X, \\theta)] \\to \\min \\limits_{q\\in Q} \\Rightarrow q(T) = p(T|X, \\theta)$<br>\n",
    "<b>M-step</b>:<br> \n",
    "$\\mathcal{L}(\\theta, q) \\to \\max\\limits_{\\theta} \\Leftrightarrow \\mathbb{E}_{q(T)}\\log p(X,T | \\theta) \\to \\max\\limits_{\\theta}$\n",
    "\n",
    "For GMM, $\\theta$ is a set of parameters that consists of mean vectors $\\mu_c$, covariance matrices $\\Sigma_c$ and priors $\\pi_c$ for each component.\n",
    "\n",
    "Latent variables $T$ are indices of components to which each data point is assigned. $T_i$ (cluster index for object $i$) is a binary vector with only one active bit in position corresponding to the true component. For example, if we have $C=3$ components and object $i$ lies in first component, $T_i = [1, 0, 0]$.\n",
    "\n",
    "The joint distribution can be written as follows: $p(T, X \\mid \\theta) =  \\prod\\limits_{i=1}^N p(T_i, X_i \\mid \\theta) = \\prod\\limits_{i=1}^N \\prod\\limits_{c=1}^C [\\pi_c \\mathcal{N}(X_i \\mid \\mu_c, \\Sigma_c)]^{T_{ic}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### E-step\n",
    "In this step we need to estimate the posterior distribution over the latent variables with fixed values of parameters: $q(T) = p(T|X, \\theta)$. We will assume that $T_i$ (cluster index for object $i$) is a binary vector with only one '1' in position corresponding to the true component. To do so we need to compute $\\gamma_{ic} = P(T_{ic} = 1 \\mid X, \\theta)$. Note that $\\sum\\limits_{c=1}^C\\gamma_{ic}=1$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Important trick 1:</b> It is important to avoid numerical errors. At some point you will have to compute the formula of the following form: $\\frac{e^{x_i}}{\\sum_j e^{x_j}}$. When you compute exponents of large numbers, you get huge numerical errors (some numbers will simply become infinity). You can avoid this by dividing numerator and denominator by $e^{\\max(x)}$: $\\frac{e^{x_i-\\max(x)}}{\\sum_j e^{x_j - \\max(x)}}$. After this transformation maximum value in the denominator will be equal to one. All other terms will contribute smaller values. This trick is called log-sum-exp. So, to compute desired formula you first subtract maximum value from each component in vector $X$ and then compute everything else as before.\n",
    "\n",
    "<b>Important trick 2:</b> You will probably need to compute formula of the form $A^{-1}x$ at some point. You would normally inverse $A$ and then multiply it by $x$. A bit faster and more numerically accurate way to do this is to solve the equation $Ay = x$. Its solution is $y=A^{-1}x$, but the equation $Ay = x$ can be solved by Gaussian elimination procedure. You can use ```np.linalg.solve``` for this.\n",
    "\n",
    "<b>Other usefull functions: </b> <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.slogdet.html\">```slogdet```</a> and <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.det.html#numpy.linalg.det\">```det```</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Task 1:</b> Implement E-step for GMM using template below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def E_step(X, pi, mu, sigma):\n",
    "    \"\"\"\n",
    "    Performs E-step on GMM model\n",
    "    Each input is numpy array:\n",
    "    X: (N x d), data points\n",
    "    pi: (C), mixture component weights \n",
    "    mu: (C x d), mixture component means\n",
    "    sigma: (C x d x d), mixture component covariance matrices\n",
    "    \n",
    "    Returns:\n",
    "    gamma: (N x C), probabilities of clusters for objects\n",
    "    \"\"\"\n",
    "    N = X.shape[0] # number of objects\n",
    "    C = pi.shape[0] # number of clusters\n",
    "    d = mu.shape[1] # dimension of each object\n",
    "    gamma = np.zeros((N, C)) # distribution q(T)\n",
    "#     print(N,C,d,gamma)\n",
    "    \n",
    "    \n",
    "    for n in range(N):\n",
    "        x = X[n]\n",
    "        for c in range(C):\n",
    "            diff = np.reshape(x - mu[c],(d,))\n",
    "#             print(x,mu[c],diff,sigma[c])\n",
    "#             print(np.linalg.slogdet(sigma[c]))\n",
    "#             print(np.dot(np.dot(np.transpose(diff),np.linalg.inv(sigma[c])),diff),\"\\n\")\n",
    "            \n",
    "            gamma[n][c] = np.log(pi[c]) - (0.5)*np.linalg.slogdet(sigma[c])[1] - (0.5)*np.dot(np.transpose(diff),\n",
    "                                                                        np.linalg.solve(sigma[c],diff))\n",
    "#         print(gamma[n])\n",
    "        gamma[n] = gamma[n] - np.max(gamma[n])\n",
    "#         print(gamma[n])\n",
    "        gamma[n] = np.exp(gamma[n])\n",
    "#         print(gamma[n])\n",
    "        gamma[n] = gamma[n]/np.sum(gamma[n])\n",
    "#         print(gamma[n])\n",
    "        \n",
    "        \n",
    "#         break\n",
    "    \n",
    "    return gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = E_step(X, pi0, mu0, sigma0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current answer for task Task 1 (E-step) is: 0.5337178741081263\n",
      "Submitted to Coursera platform. See results on assignment page!\n"
     ]
    }
   ],
   "source": [
    "grader.submit_e_step(gamma)\n",
    "grader.submit(\"goutham7r@gmail.com\", \"oLKeOUdLHBd2oQ6x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### M-step\n",
    "\n",
    "In M-step we need to maximize $\\mathbb{E}_{q(T)}\\log p(X,T | \\theta)$ with respect to $\\theta$. In our model this means that we need to find optimal values of $\\pi$, $\\mu$, $\\Sigma$. To do so, you need to compute the derivatives and \n",
    "set them to zero. You should start by deriving formulas for $\\mu$ as it is the easiest part. Then move on to $\\Sigma$. Here it is crucial to optimize function w.r.t. to $\\Lambda = \\Sigma^{-1}$ and then inverse obtained result. Finaly, to compute $\\pi$, you will need <a href=\"https://www3.nd.edu/~jstiver/FIN360/Constrained%20Optimization.pdf\">Lagrange Multipliers technique</a> to satisfy constraint $\\sum\\limits_{i=1}^{n}\\pi_i = 1$.\n",
    "\n",
    "<br>\n",
    "<b>Important note:</b> You will need to compute derivatives of scalars with respect to matrices. To refresh this technique from previous courses, see <a href=\"https://en.wikipedia.org/wiki/Matrix_calculus\"> wiki article</a> about it . Main formulas of matrix derivatives can be found in <a href=\"http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf\">Chapter 2 of The Matrix Cookbook</a>. For example, there you may find that $\\frac{\\partial}{\\partial A}\\log |A| = A^{-T}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Task 2:</b> Implement M-step for GMM using template below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def M_step(X, gamma):\n",
    "    \"\"\"\n",
    "    Performs M-step on GMM model\n",
    "    Each input is numpy array:\n",
    "    X: (N x d), data points\n",
    "    gamma: (N x C), distribution q(T)  \n",
    "    \n",
    "    Returns:\n",
    "    pi: (C)\n",
    "    mu: (C x d)\n",
    "    sigma: (C x d x d)\n",
    "    \"\"\"\n",
    "    N = X.shape[0] # number of objects\n",
    "    C = gamma.shape[1] # number of clusters\n",
    "    d = X.shape[1] # dimension of each object\n",
    "    \n",
    "    prob_sums = np.reshape(np.sum(gamma,axis=0),(C,1))\n",
    "    mu = np.divide(np.dot(np.transpose(gamma),X),prob_sums)\n",
    "    \n",
    "    sigma = np.zeros((C,d,d))\n",
    "    pi = np.zeros(C)\n",
    "        \n",
    "    for c in range(C):\n",
    "        tmp_gamma = np.array([gamma[:,c],gamma[:,c]])\n",
    "        sigma[c] = np.transpose(np.matmul(tmp_gamma*np.transpose(X-mu[c]),(X-mu[c])))/prob_sums[c]\n",
    "        pass\n",
    "    \n",
    "    pi = np.reshape(prob_sums,(C,))\n",
    "#     print(pi)\n",
    "    pi = pi/np.sum(pi)\n",
    "#     print(pi)\n",
    "    \n",
    "    \n",
    "    return pi, mu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3,) (3, 2) (3, 2, 2)\n"
     ]
    }
   ],
   "source": [
    "pi, mu, sigma = M_step(X, gamma)\n",
    "print(pi.shape,mu.shape,sigma.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current answer for task Task 2 (M-step: mu) is: 2.8993918820503812\n",
      "Current answer for task Task 2 (M-step: sigma) is: 5.977105216897527\n",
      "Current answer for task Task 2 (M-step: pi) is: 0.5507624459218776\n",
      "Submitted to Coursera platform. See results on assignment page!\n"
     ]
    }
   ],
   "source": [
    "gamma = E_step(X, pi0, mu0, sigma0)\n",
    "pi, mu, sigma = M_step(X, gamma)\n",
    "grader.submit_m_step(pi, mu, sigma)\n",
    "grader.submit(\"goutham7r@gmail.com\", \"oLKeOUdLHBd2oQ6x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loss function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we need some function to track convergence. We will use variational lower bound $\\mathcal{L}$ for this purpose. We will stop our EM iterations when $\\mathcal{L}$ will saturate. Usually, you will need only about 10-20 iterations to converge. It is also useful to check that this function never decreases during training. If it does, you have a bug in your code.\n",
    "\n",
    "<b>Task 3:</b> Implement a function that will compute $\\mathcal{L}$ using template below.\n",
    "\n",
    "$$\\mathcal{L} = \\sum_{n=1}^{N} \\sum_{k=1}^{K} \\mathbb{E}[z_{n, k}] (\\log \\pi_k + \\log \\mathcal{N}(x_n | \\mu_k, \\sigma_k)) - \\sum_{n=1}^{N} \\sum_{k=1}^{K} \\mathbb{E}[z_{n, k}] \\log \\mathbb{E}[z_{n, k}]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_vlb(X, pi, mu, sigma, gamma):\n",
    "    \"\"\"\n",
    "    Each input is numpy array:\n",
    "    X: (N x d), data points\n",
    "    gamma: (N x C), distribution q(T)  \n",
    "    pi: (C)\n",
    "    mu: (C x d)\n",
    "    sigma: (C x d x d)\n",
    "    \n",
    "    Returns value of variational lower bound\n",
    "    \"\"\"\n",
    "    N = X.shape[0] # number of objects\n",
    "    C = gamma.shape[1] # number of clusters\n",
    "    d = X.shape[1] # dimension of each object\n",
    "    \n",
    "    loss = 0\n",
    "    \n",
    "    for n in range(N):\n",
    "        for c in range(C):\n",
    "            diff = X[n] - mu[c]\n",
    "            if pi[c]>0 and gamma[n][c]>0:\n",
    "                loss += gamma[n][c]*(np.log(pi[c]) - 0.5*d*np.log(2*np.pi) \n",
    "                                         - 0.5*np.linalg.slogdet(sigma[c])[1] \n",
    "                                         - 0.5*np.dot(np.transpose(diff),np.linalg.solve(sigma[c],diff))\n",
    "                                         - np.log(gamma[n][c]))\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current answer for task Task 3 (VLB) is: -1213.973464306019\n",
      "Submitted to Coursera platform. See results on assignment page!\n"
     ]
    }
   ],
   "source": [
    "pi, mu, sigma = pi0, mu0, sigma0\n",
    "gamma = E_step(X, pi, mu, sigma)\n",
    "pi, mu, sigma = M_step(X, gamma)\n",
    "loss = compute_vlb(X, pi, mu, sigma, gamma)\n",
    "grader.submit_VLB(loss)\n",
    "grader.submit(\"goutham7r@gmail.com\", \"oLKeOUdLHBd2oQ6x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bringing it all together"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have E step, M step and VLB, we can implement training loop. We will start at random values of $\\pi$, $\\mu$ and $\\Sigma$, train until $\\mathcal{L}$ stops changing and return the resulting points. We also know that EM algorithm sometimes stops at local optima. To avoid this we should restart algorithm multiple times from different starting positions. Each training trial should stop either when maximum number of iterations is reached or when relative improvement is smaller than given tolerance ($|\\frac{\\mathcal{L}_i-\\mathcal{L}_{i-1}}{\\mathcal{L}_{i-1}}| \\le \\text{rtol}$).\n",
    "\n",
    "Remember, that values of $\\pi$ that you generate must be non-negative and sum up to 1. Also, $\\Sigma$ matrices must be symmetric and positive semi-definite. If you don't know how to generate those matrices, you can use $\\Sigma=I$ as initialization.\n",
    "\n",
    "You will also sometimes get numerical errors because of component collapsing. The easiest way to deal with this problems is to simply restart the procedure.\n",
    "\n",
    "<b>Task 4:</b> Implement training procedure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def train_EM(X, C, rtol=1e-3, max_iter=100, restarts=10):\n",
    "    '''\n",
    "    Starts with random initialization *restarts* times\n",
    "    Runs optimization until saturation with *rtol* reached\n",
    "    or *max_iter* iterations were made.\n",
    "    \n",
    "    X: (N, d), data points\n",
    "    C: int, number of clusters\n",
    "    '''\n",
    "    N = X.shape[0] # number of objects\n",
    "    d = X.shape[1] # dimension of each object\n",
    "    best_loss = -float(\"inf\")\n",
    "    best_pi = None\n",
    "    best_mu = None\n",
    "    best_sigma = None\n",
    "    \n",
    "    avg = np.mean(X,axis=0)\n",
    "    \n",
    "    for _ in range(restarts):\n",
    "        print(\"\\n\\nStart\",_)\n",
    "        \n",
    "        pi = np.random.random(C,)\n",
    "        pi = pi/np.sum(pi)\n",
    "        \n",
    "        mu = np.random.random((C,d))*avg\n",
    "        \n",
    "        sigma = np.zeros((C,d,d))\n",
    "        for c in range(C):\n",
    "            sig = np.random.randn(d,d)\n",
    "            sigma[c] = np.dot(sig,sig.transpose())\n",
    "#             sigma[c] = np.eye(d)\n",
    "        \n",
    "        try:\n",
    "            for i in range(max_iter):\n",
    "                print(\"iter\",i)\n",
    "                gamma = E_step(X, pi, mu, sigma)\n",
    "                loss1 = compute_vlb(X, pi, mu, sigma, gamma)\n",
    "                pi, mu, sigma = M_step(X, gamma)\n",
    "                loss2 = compute_vlb(X, pi, mu, sigma, gamma)\n",
    "                print(loss1,loss2)\n",
    "                \n",
    "                perc = (1.0*abs((loss1-loss2)/loss1))\n",
    "                if perc<rtol:\n",
    "                    break\n",
    "            \n",
    "            if loss2>best_loss:\n",
    "                best_loss, best_pi, best_mu, best_sigma = loss2, pi, mu, sigma\n",
    "            print(best_loss)\n",
    "\n",
    "        except np.linalg.LinAlgError:\n",
    "            print(\"Singular matrix: components collapsed\")\n",
    "            pass\n",
    "\n",
    "    return best_loss, best_pi, best_mu, best_sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Start 0\n",
      "iter 0\n",
      "-1431.1225359817538 -1193.7263305309043\n",
      "iter 1\n",
      "-1179.1313632035553 -1172.290389876311\n",
      "iter 2\n",
      "-1167.5972493402003 -1163.946729608073\n",
      "iter 3\n",
      "-1160.7005459383038 -1157.6951415318158\n",
      "iter 4\n",
      "-1154.6299580364628 -1151.6326966637882\n",
      "iter 5\n",
      "-1148.4942713345943 -1145.3227061045125\n",
      "iter 6\n",
      "-1141.9198945294859 -1138.2718086118316\n",
      "iter 7\n",
      "-1134.4364613760697 -1130.3230564660614\n",
      "iter 8\n",
      "-1126.4258588359655 -1122.3059522382716\n",
      "iter 9\n",
      "-1118.4203410444632 -1114.0063701520874\n",
      "iter 10\n",
      "-1109.594980388274 -1104.0024967630259\n",
      "iter 11\n",
      "-1098.5469211543048 -1091.2684930387768\n",
      "iter 12\n",
      "-1085.6905047148075 -1079.1282857223264\n",
      "iter 13\n",
      "-1074.437767987172 -1069.601478663527\n",
      "iter 14\n",
      "-1067.505200058602 -1066.2598620866797\n",
      "iter 15\n",
      "-1065.5127985337772 -1064.9758118895377\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 1\n",
      "iter 0\n",
      "-2245.6632687575266 -1244.7212638380363\n",
      "iter 1\n",
      "Singular matrix: components collapsed\n",
      "\n",
      "\n",
      "Start 2\n",
      "iter 0\n",
      "-1443.1913462370283 -1248.9426387576634\n",
      "iter 1\n",
      "-1217.420304929371 -1192.9713507530225\n",
      "iter 2\n",
      "-1178.0501863703923 -1168.2102050749859\n",
      "iter 3\n",
      "-1164.5440181395074 -1162.7617676834818\n",
      "iter 4\n",
      "-1162.1191354477194 -1161.8193862564744\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 3\n",
      "iter 0\n",
      "-13902.759126400799 -1226.4482682932025\n",
      "iter 1\n",
      "-1172.0762807118049 -1162.2420102619897\n",
      "iter 2\n",
      "-1155.276961064641 -1148.8403527762373\n",
      "iter 3\n",
      "-1144.2885743571956 -1140.6261875486891\n",
      "iter 4\n",
      "-1138.5216014214564 -1137.12268397898\n",
      "iter 5\n",
      "-1136.290494194074 -1135.6159539684736\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 4\n",
      "iter 0\n",
      "-4836.318827012205 -1260.903561149659\n",
      "iter 1\n",
      "-1204.5466032238758 -1188.7303902304454\n",
      "iter 2\n",
      "-1178.8111139630703 -1169.1935835818608\n",
      "iter 3\n",
      "-1163.4998177104205 -1159.982849502748\n",
      "iter 4\n",
      "-1158.4344815257891 -1157.505018646558\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 5\n",
      "iter 0\n",
      "-1822.9891809604649 -1243.3843399309958\n",
      "iter 1\n",
      "-1205.3224431091187 -1190.0627102410838\n",
      "iter 2\n",
      "-1177.7277647199323 -1165.8163659291392\n",
      "iter 3\n",
      "-1154.0575373858178 -1142.7745922274826\n",
      "iter 4\n",
      "-1132.752076541957 -1124.4406323128721\n",
      "iter 5\n",
      "-1118.4058389306522 -1113.9711128437657\n",
      "iter 6\n",
      "-1110.8991466019263 -1108.5534522699047\n",
      "iter 7\n",
      "-1106.7475826767889 -1105.2680096252607\n",
      "iter 8\n",
      "-1104.0498154291786 -1103.0169592159016\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 6\n",
      "iter 0\n",
      "-1575.8285693103505 -1202.2402082573908\n",
      "iter 1\n",
      "-1182.8597122196752 -1175.8191997201916\n",
      "iter 2\n",
      "-1170.7245524043865 -1167.0772448075813\n",
      "iter 3\n",
      "-1164.1610129481207 -1161.7274147567284\n",
      "iter 4\n",
      "-1159.5757228036175 -1157.520691424848\n",
      "iter 5\n",
      "-1155.5980534512464 -1153.6526155128686\n",
      "iter 6\n",
      "-1151.8592089664799 -1150.0526858780752\n",
      "iter 7\n",
      "-1148.4224028347105 -1146.8211457205591\n",
      "iter 8\n",
      "-1145.46363810166 -1144.2518507937205\n",
      "iter 9\n",
      "-1143.3458897916155 -1142.6381702858623\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 7\n",
      "iter 0\n",
      "-1692.3652073574485 -1235.0214787330156\n",
      "iter 1\n",
      "-1201.8371681013716 -1182.3403162332813\n",
      "iter 2\n",
      "-1167.1515958068032 -1155.3202528766708\n",
      "iter 3\n",
      "-1145.8543617923763 -1138.5799161055481\n",
      "iter 4\n",
      "-1133.3522291863974 -1129.0081436661633\n",
      "iter 5\n",
      "-1125.316302322203 -1121.782913600516\n",
      "iter 6\n",
      "-1118.4794998410794 -1115.1175687932166\n",
      "iter 7\n",
      "-1111.9687619205845 -1108.748933059172\n",
      "iter 8\n",
      "-1105.767949418826 -1102.7003284350856\n",
      "iter 9\n",
      "-1099.8365143748854 -1096.7811730981275\n",
      "iter 10\n",
      "-1094.083021654991 -1091.1729848019477\n",
      "iter 11\n",
      "-1088.8673518455548 -1086.5282892675548\n",
      "iter 12\n",
      "-1084.8072595758133 -1083.1263168864791\n",
      "iter 13\n",
      "-1081.7802827048927 -1080.314923346102\n",
      "iter 14\n",
      "-1078.966017671016 -1077.2495133289535\n",
      "iter 15\n",
      "-1075.572200059535 -1073.2776130634677\n",
      "iter 16\n",
      "-1071.3741944859123 -1069.172834484769\n",
      "iter 17\n",
      "-1067.7737819912272 -1066.5684002775165\n",
      "iter 18\n",
      "-1065.9182765319388 -1065.3963138660329\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 8\n",
      "iter 0\n",
      "-3104.7492820069574 -1297.4239376234798\n",
      "iter 1\n",
      "-1205.8479505258406 -1176.941484757821\n",
      "iter 2\n",
      "-1162.969968193041 -1152.9512093843912\n",
      "iter 3\n",
      "-1146.6361735960681 -1141.6008795805049\n",
      "iter 4\n",
      "-1138.418590609205 -1136.1397295973345\n",
      "iter 5\n",
      "-1134.613711104865 -1133.4298831985047\n",
      "iter 6\n",
      "-1132.4544127892843 -1131.5453213106707\n",
      "-1064.9758118895377\n",
      "\n",
      "\n",
      "Start 9\n",
      "iter 0\n",
      "-4141.742455248762 -1288.0180108908735\n",
      "iter 1\n",
      "-1193.4891602747662 -1169.332906756382\n",
      "iter 2\n",
      "-1157.2553118823407 -1149.4119082764364\n",
      "iter 3\n",
      "-1145.7232470502104 -1142.4849101891684\n",
      "iter 4\n",
      "-1139.9545743471929 -1137.2070987078762\n",
      "iter 5\n",
      "-1134.8480027565013 -1132.3326886969426\n",
      "iter 6\n",
      "-1129.761583583516 -1126.8912843889045\n",
      "iter 7\n",
      "-1123.4809126725684 -1119.5572988839237\n",
      "iter 8\n",
      "-1115.1961705128306 -1111.1995617807975\n",
      "iter 9\n",
      "-1107.4545830819457 -1103.5315515410505\n",
      "iter 10\n",
      "-1099.1693119184476 -1094.5321853618705\n",
      "iter 11\n",
      "-1090.890134536773 -1087.8525716905224\n",
      "iter 12\n",
      "-1085.4705081263342 -1083.2429463631702\n",
      "iter 13\n",
      "-1081.1096378634463 -1078.9668822812232\n",
      "iter 14\n",
      "-1076.8589372807342 -1074.768925322694\n",
      "iter 15\n",
      "-1072.819874662163 -1070.9833791293418\n",
      "iter 16\n",
      "-1069.4023560879095 -1068.0197970582274\n",
      "iter 17\n",
      "-1066.9245652719085 -1066.0407192356577\n",
      "-1064.9758118895377\n",
      "Current answer for task Task 4 (EM) is: -1064.9758118895377\n"
     ]
    }
   ],
   "source": [
    "best_loss, best_pi, best_mu, best_sigma = train_EM(X, 3,restarts=10)\n",
    "grader.submit_EM(best_loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you implemented all the steps correctly, your algorithm should converge in about 20 iterations. Let's plot the clusters to see it. We will assign a cluster label as the most probable cluster index. This can be found using matrix $\\gamma$ computed on last E-step. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXm4TdUbxz9rn/ncyTxnnud5SDIlUSoVhUq/SKMGKkqljEVSkUSa0EAqIaJIMmUOmef5Gu58z7zX74+jw3HOuQMX17U+z9PD3WdN+3S9e+13ve/3FVJKFAqFQpF30K72AhQKhUKRsyjDrlAoFHkMZdgVCoUij6EMu0KhUOQxlGFXKBSKPIYy7AqFQpHHUIZdoVAo8hjKsCsUCkUeQxl2hUKhyGMYr8akhQoVkmXLlr0aUysUCsU1y7p1605JKQtn1u6qGPayZcuydu3aqzG1QqFQXLMIIQ5kpZ1yxSgUCkUeQxl2hUKhyGMow65QKBR5DGXYFQqFIo+hDLtCoVDkMa5KVIzi2iM9xcGPH85j5Zx1FC1diK4v302VhhWu9rIUCkUYlGFXZIrH7eHZZq9ybO8J3E4PO9fsYfUv6xkyeyD129YKtDtzPAGj2UhsgZiruFqFQpEjrhghxAtCiK1CiC1CiG+EENacGFeRO1j+49/EHzyF2+kBQEqJK93NxH5fAHB451H61OnPg+We5oESfRhw61CSz6RcxRUrFNc3l2zYhRAlgWeBhlLKmoABeOBSx1XkHvZtPogj1Rly/ciuY/h8Pl5s8yb7txzC4/LgcXv5Z+lWhnZ57yqsVKFQQM65YoyATQjhAezA0RwaV5ELKF+7DLZoa4hxL1W5BFuWbSc9xcn5RdG9Hh9bV2wnIT6J/EXiANi/9RCTB0xlx9+7KVGxGL1G9KBOqxoc2HaYWWPncnxfPM06NaDjY7dgsVmu6P0pFHmNSzbsUsojQoh3gYOAA1gopVx4YTshRB+gD0Dp0qUvdVrFFeTGuxsxbdj3HN19HLfTgxACs83Ek2MfwZnuQojQPkLT8DjdAMQfOsVzNw7CkepASkg6lcKg20fw+JieTOz/BW6H38WzcfEWfhq3gCn/jsVoUsc/CsXFkhOumPzAXUA5oAQQJYR48MJ2UspJUsqGUsqGhQtnqmGjyEWYzCY+XDGch97oQvVmlWl5/428v2wYdVvXpHbL6ug+PaRP4VIFKXxDIQBmj1+A2+XmvE09LoebSS99FTDq4PfdH91znG9G/nDZ70mhyMvkxOHpLcA+KeVJKaUH+AG4MQfGVeQibNE2HhjYmQ+WD2fQ189TsV45//UoK6/P6I/VbsEea8MeayOucCyDZ72IOLuVP7zzKF63L2RMZ5or7FwLPlty+W5EobgOyIn33YNAUyGEHb8rpi2gpBuvIxp3qMeM45PZ8PsWTFYT9drUxGgycmjHEWaNncuh7UcwmAz4POeMu8Ho31P4vKG7fZfDfcXWrlDkRXLCx75aCPE9sB7wAhuASZc6ruLawhZt48a7GgGwe+M+PnzqU7at3gln3S9CgBACKSUWmxlbjI3Gt9dj4ed/hIzVrFPDK7hyhSLvkSMnVFLKwcDgnBhLcW2zf+shXmjxeoibRUowWQxUb1aFZp0a0P5/bTBbTfy7fAdH95xA9+kYjAZiCkbzyFAVLatQXAoq9ECRo3wz8kdc6eF95x6Xl4r1y3HvC50C1yb9M4Zls1azZfl2ylQvxS0P3kxUrP1KLVehyJMow67IEl6Plx8+mMeSb5ZToHg+7uvXiXptaoW0O7jtcFD0y/loBo2yNYJDXU1mE2263USbbjfl2FqTTiUzb/Jv7N98kFotqtGuZyusdhUbr7h+UIZdEUJ6ioP1v/2DwWigQbvaaAaNB8s9xemjCf4GG2D9b5t57uPHuO1/bYL61mldgz2b9iP1UOsudUm9NjX49YslmCwmmt7RAHuMLUfXfuLASZ5q+DLONBdup4cVP6/lx3Hz+WjN29iilNKF4vpAyEjbq8tIw4YNpap5mjtZ/9s/DO48Cs3gj1oRmqBu65os//HvkLZRcXZmnfoMg8EQuJYQn8QjlfuSnuwIaZ+vSCxpyQ4MBg0hBEITjFzwGtWbVg5pu/qX9Xw78kcSTybR/O7GdHulM1FxUZmuf9Qj4/l9+rKg2HqL3Uzvtx/k7mc6ZOk7UChyK0KIdVLKTKMLlB67IoDb6eat+97FmeYiPdlBerKDtMR0VvwUatQBHKlO0hLTg67lLxLHpE1jsEVb4YKM1MT4ZDxOD840F45UJ+nJDoZ2HcOFm4tFU5cytOt7bFm+ncM7j/HDB/N4rvlr+LyhsfAXsvnPbSEJU650Nxt+35yFb0ChyBsow64IsG3VrrDXI73UGYwGovKFHnQWLVOYz7a9z60PtyImfzRCC6M5cJbk06kc3XMcgJSEVKYNncl7j00MOoD1uLzEHzrF3/M3ZHoPJSsXD7lmMhspU6NUpn0ViryCMuyKANYoS1jfeCTDXKl+eRJOJIX9rFDJgrz0+dOUrFw87Jj/4XF5sMfaSUtK44l6L/HNyB/xur2h7ZweDu88luk9PPxmVyx2c9DazTYzdz51W6Z9FYq8gjLsigCVG1Ygf9G4gH8dwGgyYLaZw7bf/vcuHqn8LBsWR3ZzhNOROR+z1Uz+InH88unvJJ1MDmi+X4jRbKRyw/KZ3kP1ppV5Z+Eb1GlVg4IlCnBT58aMXz2SQiUKZNpXocgrqKgYRQAhBKN+G8zQ+99jz8b9ANS8qSpHdx8nPoyui+7TcaW7eOu+d7HYzCScSKJ87TI8O+GxwIFox95tObjtcERdmLI1bwDg35U7I0oJWKMs1Gheldo3V8/SfdS4sQrvLn4TAEeak2N7TpCamEZ0vswPXxWKvIDasSuCKFqmMG/9+DJPju3JS58/zYhfXuWme5pgNEfeA6QlpnPmWCJSl+zZuJ8B7YZw4sBJADr0bkvLrjeGleG12C10f+UeACrVL4fZagppk69IHH3H92b43FcComJZ5YcP5tKlaG9eaPE695d4jIn9vgg5qFUo8iLKsCuCmPvJQh6u8DSfDpjOB09OokeZp2j7YAtKVCga0SVzIV6Pl3mTFgGgaRovTnmKqXvH02tEdwqVKoimCWILxfDk2J4BfZnb+7TDHmPDYPKHTgpNYLVbGLlgELf2bIXBaIg4Xzg2/bGVzwZ9iyvdRXqKA7fTw7xJv7Hgs8XZGkehuBZRceyKACcPn+aRyn2D/NxCE1RtVJGxfw3l+/fmMuWV6Rkehv5Hu54tefnzZ0KuSynxuDyYLKaQHfipo2f4ZuSP/LN0K2Wql6LHoHspV6tM4PM1v25kzse/4nV7af+/Ntx8X9OIu/gRPT5gyTd/hVyvWK8cH68blen6FYrcSFbj2JWPXRFg7a8bgw5OwZ8tun3NbpJPpzJ92PchRt0eY8PpcKGfJ79rNBsxGAwc23eC4uWKApCe6mDFT2s4svsYJSoUo8W9TUPS/AuVKEDfcb3Cru270bOZNmRmwFe/edk2/vlzK33H9Q7b3usOfwiblVh4heJaRxl2RQB7jC1saKPUJbPHzQ/bx+VwUahEAVLOpOJMdyF1ie7TWfz1MhZ/vYzHx/TE6/byyYtfBSJkDCYDn7z4FeNXj6RY2SKBsdJTHBzfF0+xckWCpAYcaU6mvjkj6HDVmeZi/qeL6TawM4VKFgxZV/tHWrNmwcagQ1uL3UyHXm1C2ioUeQ3lY1cEaHJHg4i1RmeO+TmCC0YwccNo+k1+IhB1ovt03E4PbqeHj1/4gskDpwWFPfo8PpJOJvPRs58BfvfM1CEzAwedXYr2ZuqQmYGDzuP74tHC+NjNVhP7thwKu97GHetz51PtMVlMRMXZMVtNNOvUUMWzK64L1I5dEcBqtzBi/iD6Nnkl5DOD0YDHFeze0AwatVpUIyZ/NHVa1Qgb0ih1GdH9sXHxFgBWzF7DjNGzcTvduJ3+z2aMnk2xckVo26MFRUoXwucJTVpypDg4sPUQdVvXwGQOjqgRQvDYOw9x7wt3sG/zQUpULBZwCykUeR21Y1cEUaZ6qfARKALufrYj1igL9lg71mgrJSsVZ8DUvoA/1jy7ROf37/B/nrAg5KHgTHMx+pHxdCnWm1Vz13Ffv04hc+i65MvB3/Fi6zfxhjH8AAWK5adBuzrKqCuuK5RhVwRhi7JS86aqgZqk/6H7JD0G3cvME1MYPOtF3l38JlO2jg1kdNqibbS8/8agkEihCewxNkzW0BdDg8nAAwPvJv7QKTYv2xZ2LVJC8qkUxvb5hMYd69Fv8pNEXZBk5ExzsW/zQZb/tOZSb12Ri5F6CtK1FOnZrHIRsoAy7IoQXpn+HKWqlDi7O7dhj7Ux+Pv+ROeLwmq3UL9tLao0rIAQAkeak41LtrD3nwM8P7EP7R9pFXgoCCFwO92YTMFuEqPZSK8R3bnzqdsY0+tjvO6MI1XcDhc/jV9A6weaEyiieh6OVCeb/9yaY/evyF3o6T8i45sjE19AnnkIefpOpH7mai8rV6N87LmU3Rv3sfDLP/B5fbTtcXNYzfLLRcHi+Zn8z3vs/ecAaUnpVG1SCbMlNCv0z+9XMvqRjzAYNXw+nZKVivPiZ0/x6xd/4PO6/ZIDF8gEmG1mur96D13634mu62xcsiXTHZiUkHImFYAipQuzL/FA0OcWu4Ubqpa8xLtW5Eak7ygkvwG4zj3TvXuQSa8j8n90NZeWq1E79lzI/M9+5/nmrzF73HzmTFjIy7e8xXejZ1/RNQghqFCnLLVvrh7WqJ8+lsA7PcfjTHeRluzAmeZi/9ZDjOo5PrLOL+B2uP11UR0uhBAZShX8hzXKQtvuLQDoNaJ7iPSA2+mmfO0y4bpmCV3X+WbkD3Qp1ptOMQ8ypOsYzhxPuOjxFJkj9UT0lA/RT/dAT3oL6Q0f3YTztzAXveBaolwyGaAMey7D5XAx4bnPcTnc6LpESokr3c1Xg78j+UzK1V5egBWz13Bh0qfP42P/1gj/QM9DCDh9NAEhBB0ebRMiVaAZNCx2C9YoCxa7hQa31qFND39N1CYd61PwAqVGqUveum/MRScffTboG6YP/4HE+CScaS5W/PQ3z92YtcIeiuwj9RTkqTshbRJ41oDjO797xbs7tLEwE1KxBYDsSUxcbyhXTC7j0I6jaFro89ZkMbFn4/6wBaQz48juY8wcM4eD2w5Tv20t7u7bEVu0laN7jhNXOJbYAjEhfXw+H5+/9i2Lpy9DSknrbs2586n2zP1kERt+34LBZAi7Y9I0Da8nY4MohKBQKX9SUZ93H8aZ7mLx18uQEoqWKcTAac8hdcnBbYepULcsFeuWC/SNP3iS00dD/asel4ctf22nTqsa2fpufF4fP437BVe6+7xrOkmnU1i36B8ad6iXrfEUmSPTvwM9AfjvO/eC9CFTPkDkHxfc2HorpIy84GjFArY7sy0Kdz2hDHsuo1DJAnjCFJrwur0ULVM42+Pt23yA55q/htvpxufV2fH3bmZ/9Ctejxev24vPq9OyazP6TX4iEAuecCKR3jX7kXz63BvCrLHz+OGDXxCA1+NDM2ghWusmi5EW9zZjxc9/40wNL9NrsVvoNbJ7wL1jtph4ccpTPDOuF840J3GFYgP/YKs1qRTSX0oIeVU4i65nrP0eDpfDHbawh9R1Th1RB3SXBc8G4MLfDwmeUF1/oRWAfBORSf1AOkB6wXIzIva1K7LUaxVl2HMZ+QrH0fqB5iyduSKwizTbzNRtU5MSFYple7xPX/kaZ5oz4Pb2Z4QGVz368/tVFC5VkEeHdwfgw6c/DTLqcLZghu+CnwUYjQZMFhM+r49aLarx3MePsW3VTo6lngjqL4SgdPVSPP3B/6jXphYnDpxk5Zy1mK1mmnduRNLJFCw2c6a7sKJlClO8fFEObjsclAlrMBqo1aJadr8ebNFWipUrypFdwdWZdJ+kTqus6b8rsompJrj+JNi4CzBWCdtcWJpB4eXg2w9aPr+xV2SIUnfMhfi8Pma++zNzJy1C9+rc2rMl3V+7L+whZmbcX7IPZ45lfhAYVyiG7+P9Kf6327tHrGR0IQVL5Gfk/EHEFIgOaLbMn/I7Hz79achOuPuge/jf0G7Mn/I74/tOCVx3uzz+exOC8rXLMOSnl8lfNF/EOY/uOc7A24aReCIJIQQmq4lnxj3K+kWbObzrKPYYG/s2H8TtcNPi3qb0GtmdqLiooPlWz11HYnwS9drWIvl0CgPbD8Pn9eFxezFbTXTu25FeI3tk6TtQZA+pJyBPdQQ9CfACGggLosC3CFP2H87XE1lVd1SGPY8gpWTbqp3s23yQcrVKU61pZYQQDGg/lPWL/sm0f1ScnZ8SvgTg7gI9SUtMz7SPENCoQ32Gz32FP79fybRh35MUn0zdtrVY8u1fSF/w75bZauLTrWPpXeOFiA8Og9FAtWaVGbt0SKb3u2fjfjxuL1JKBrQbgtvpCXEPGc1GylQvxcfrRiGE4Ni+Ezx342s4013oXh8SuK9/J+5++jaWfLOclIRUmnZqSJWGFTK9f8XFI33xyLRPwL0GjBURUU8iTKGuN0UwSrb3OsLj9jDo9pFsW7UTKSVCCKo1rczwea/w6PDubF2+A7fTjdQlBpMB3wWHm0azkZu7NAv8fOdT7Zk5Zg5eV/CO22wzI3W/nrpm0DBbzTw6vBu/frGED5/+FPfZmPXFXy8Ll0eErussmLIYo9kY0bD7vD52rN5F0qlk4grFRrxnIQQV6/kPVfu3Ghyx9J7X7eXo7uNsXbGDms2rMrbPJySdTEI/z40za8wcWnVpxj3P3x5xPkXOIgxFELGvX+1l5FlUuGMeYN6k3/h35Q6caS5c6W6caS7+XbmDXyb/RpWGFRjzx1sUK1cEBCFGHfzhhb1Gdg/83PPN++nQqy0ms9FfySjKQq+RPRi3ajhla96AxW6haNnCvD7zBSrUKcukl6YGjDoQ1qgDeN0+5kxcmPkhpxDZCjXcu/lAxg0ExJ8t1ffP0q1BRh38D5M1CzZmeT5F9pHSh9RTVez5FUIZ9jzAkm+XB4XrAbjS3Sz+ZjkAcycu5PTRhAwMrocXbnqduZ8swuf1YTAaeHZ8b35M+ILvjk7m5+SptOvZkhdbv8Xefw7gSndxYv9JhnUdy76tB0MOWjPC7XBjsZpDtGj+Q2iCAsXz8fOEX1k6YwWeCAUzzueGKhlnnXrdXv78fiXdyzwRNnfKaDYSUyA6S+tXZA8pJXrqJGR8Q2R8Y+TJ1kjX8qu9rDyP8rFfo3jcHr55+ycWfr6ExJNJIYYdoFmnhrw2ox9353sYjyu8+uH5WOwWGrWvy+BZLwauHd8fz5heE9j4x9aQB4PQBE061mPVvPURHxrhaNzRHxu+9tdNCE0ghMBsMaFLidflwWAy4Ep3Y4u2UqRMYT5cMTyo8MaFbPlrGwNvG4bb4QnaEQrhdx/5vDo+jy/ibjEqzs70/ROCDlgVOYOe/j0kDwUc5121IgrNQRgvPlv4euWK+tiFEPmAT4Ga+P+JPyqlXJkTY18PJJ1KxuvxUbB4/iz3Gdr1PdYv+idEi+U/LHYL975wBx6nG92XNavrSnexZsEG9m89RNkaN+D1eHn+ptdIOJ4Y1nBLXXJw+xGMJkOmQl7ns+H3zXy4cgTD5vh13z1uL9tW7mTpjBUs+Hxx4CHlSHVybM9xfhr3C91fvTfieDVvqsbYP4fy9YgfOLbnBFUbV8QaY8WR4uTPWatwndWZOR+hCQRQrnYZXv7iGWXULxdpnxJs1AG8SMcMRMxLV2NF1wU55Yr5AFggpawK1AHC67Aqgkg6lcxLbd+i2w2P81D5p+lTpz+HL4inDsexvSdYt3BTWKNuMBqIKxTL8xMfo06rGkTFRVG6WtYFsiSwZ9N+wL+jdqQ6Q3zS/6EZNKo2qcztfdqF6LdkhMfl5f3HJyHEud16lcYVWTV3XcibhdvpYeWcdZmOWal+eQZ//yITN4zm+U8e54l3e6LrOmmJaWHbR+WzMydtOhPWvgPA4Z1HAfB6vMQfPInblbVwT0UmyNCHKnjPhjoqLheXvGMXQsQBNwOPAEgp3ZzLFVZkwJAuY9i6YkfgQHP/lkO82OZNpu+fgMEQWQvj+P54TBZTSGSJLcbK8LmvUqN5lSBZglemPUv/VoPxen34PL4MY9TdDjcT+39JtcYVWTF7TcQ3As2oYbVbeHhwF4qVK0K+onF8NXhGSPk8k9kYNpN259o9gQgeKSUvtX2Lk0dOh7QTQlDkhtCapllhxU9rIpTzg/K1y7Bv8wEG3z2atOR0pK4TnT+K9CQHui4RAnq8di8PDOh8UXMrzmJpB44ZwPm/czaE9dartaLrgpxwxZQDTgKfCyHqAOuA56SUQVslIUQfoA9A6dKlc2Daa5szxxPYtmpXUJSKlBJHsoOty3dQ++bIWY8V6pYNKVMHftdIpQblQ7RmytUqw9eHPmHlz2vZtW4Pcz5ZGDHlHyAxPonH673kL2sXJorGaDbQofctdH3xToxmI+t/20z5mmWwRllxpAS/dkeSH4gpEB3IMv3nz385sPVQWHePwaTR9eW7I641I6zRlogHu9tW7qRvk1eD/O4XnlNMHzaLkhWL0+Lephc1vwJEzAtIz9/gOwpSByTYbgdzi6u9tDxNTrhijEB94GMpZT0gDRh4YSMp5SQpZUMpZcPChbOveZLXcDs9aFqY9Hnh93VnRGyBGB4a3AWL3YIQoGkCi93ME+/1xGoPX6LOarfQ+oHmPDKsW4ZvAwBIf2WiiLt1TSMtMY0fP5jHwxWfYfgDYxn2wHth1y2lDDHqVruFbq+e2wkf2n4Urzf84W5MgRgqNyif8XojcM9zt2OJ8H14XN5MQ++caS5mvT/vouZW+BFaLKLgHES+8YjYVxEFv0eLG6EEvC4zObFjPwwcllKuPvvz94Qx7IpgipYpTMGSBTi6+3jQdd0nqZXBbv0/HhjQmZo3VeO3qX9iMGrc+kjrLGVLmi0m3vj+RQbf/Q4uhzuiqyIj3E4PS2esQAgNr8eLJ4vyA5pBo2CJ/Dww4G46Pdk+cL1Sg/L4POFj2xNPJHLqyBkKl8q+O+b2Prew/rfNrF24ESQhD5is4ExzZruPIhghNLA0v9rLuK64ZMMupTwuhDgkhKgipdwBtAX+vfSl5W2EELz5w0u81PatgFtF9+kMnvVixF33hdRsXpWazatme+76bWvx3dHJPFarH/EHT0VYIBmGMPq8OpB9QzlszishRTGsdnPEpCUpYXDnUUxY80625nE5XDzbbBDH9p5A9+qYrKZsG3aL3Uy7h27OVh+FH+k7jkweBu7lIKIhqjfC/rDaqV8hckpSoC8wXQhhBvYC/8uhcfM05WqW5tvDn7BxyVY8Lg9129TEFmW97PO6nW7ee+zjiLK0+YvlIz05PWxs/KWg+3SeqP8SL3/5DLf0OGcwl/+0JnCIGo5D24+wb/MBytXKetzzwi/+4OieEwH3UOCtIswDyxptRUrp15J3ezGZjfh8Puq0qsldz3TI1j0qQEo38nQX0E8COsg0SHkPKZ2I6Mev9vKuC3LEsEspNwKZBs0rQjGajDS8tc4VnXPSS1NZ+fPasDtYi93C2D+HcPpoAp+8+BWHdhwlJr+dxPjkLCs+ZoTUJaMeHk9aUhp3PeU3mmarCaPJGPZAGPwunIQTSZTLpMZI8ukUzDYzVruFdb/9E9bnbzAYMFtNCE3g9fjoO74XpauWxGw1U6FuWY7sOsaeTQcoXbVEth4keQEpdWT6NEj/EqQTrLchol9AaNnMynUtPhvmeP7vlwPSJiGj+oTdtUvpAYxqR59DqMzT65BOMQ+GiGYJIah+Y2UGfNWX4uWKhvRZOWctM0bPZuvyHTmi92G2mZiTMg1N09i35SBPNXg5YuUli83MjOOfRsw+3bflIMMfGMuR3ccRAlrc25RCJQvw4we/hIRa1m1Tk/8N60byqRRq3lSV6HwqMek/9OTRkD6NcwlFZjBWRhSclS2DK9O+RKaMJjTqWSCK/osQ5w7vpWcLMulV8O4AYYeoXoiop5WBj4BSd1REJJzAlsFkoPndTUKMusftYfb4BSyauhSr3UKJisX8fuuLOIg8H7fDgyPVyemjCfS7+Y2QqkiaJtCMGprBwLMfPxbRqLscLvq3GkzKedmlf/2wmsa3N8BkNeH1eAP6MBabmZ5vdqV608qXtPa8iJQuSJ8KnH9Y7AbfXvBsAnPdrA9mbkTYgDtjlWCjrp9BnnnI76oB/5+pk5FYEdG9L+Y2FGdRhj0Pc+rIaX794g8SjifQ5I6GNLy1DkIImt3ZiBU//R20Q9Y0QfO7G4WM8eY977Lpjy0Bf7vFbqZImUIc3xcf4qs2WU3kLxpH/IEIB7IXsPqX9Sz8fAmpF2SHmm1mbn24JaUql6B558YUK1sEAEeqgxWz1+JMc9K4Y30KlyrI6nnr8XlCs1VXz1nLuL9H8sXr37F99U5KVCjGoyO6U/MmVcghLHoS4XUjBPgOA1k37MJUHWnrDI6f8O/+LSCMiLgRwQ0dc0FeuMlwQPrnoAz7JaEMex5l2+pdDGg3xB+O6PLy6xd/0LxzEwZ+1ZdnP+rNoR1HOLY3Hk0TeN1envrgfyGl9/ZtOcjGJVuCJHld6W7OHE3EZDaF+MQFUL9tbRZNXRo2selCxvT6GLcz9IDW7XBjjbZy7wt3BK7tWr+XF9u8ia5LpK4z4fnPeWZcL4SmhZU80HWdEhWKMXT2gEzXoQC0QqBFg37huYQXTMFGXepnwDnfX4PU0hZhLMeFiNg3wdYJ6VoKogDCdifCUPCCcRIJrX0K6FlXC1WERxn2PMrYxz/BkXrutdqZ5mL5j6vZsaYDVRpV5JMN77Jz7R4STiRRo3kVYvKHHpCtnrsuWGf9P6TEaDaEGHaj2UjjjvVYOnMlDs+5DFRhEKDLEMncsGOfZd8/5zTWU5PS6NfqjZBs2fF9pzD+75EhbiGhCSrVL39FIozyCkJoEDscmfg8/vR/H6CBVtDv/zaWAkC61yATHjubReqDlA+QMS+iRfW8YDwB5oYIc2R3sLC0RKZNIVgkTMW85wRKj/0a5cC2w7x+1zt0Ld6bF27wX3hqAAAgAElEQVR+nS1/ndNd03U9yDD+h9fjY/MyfzshBFUaVaTpHQ3CGnWAhV/+Efa6MIiwO3Kv20f9W2rTa0Q3DEaDfwsv8Ee8ZfO8dee6vYG/D7nv3bASCAaTkQNbD/PM+F6YrSZsMVbssTbyF4lj4LRnszehAmFtgyg0G4zlAQOgg34UmdQPPXUiUkpkYn+Q6fh98R7ABSnvIn1Zc78FzWeuA/augBW/uyYKtMKI2Ddy8rauS9SO/Rok/tApnm32Ko4UB1JCwokkBrYfxqjf3qB6M78AWGzBaJJPByvrmSxGipY5J+eg6zq/T1vGr18uwWq3cNfTt9HotnqBzw7tOBp2/qZ3NKRKowp8OXgGbocLEJhtZh4Zcj+OVCeTX54WdEArsyPWfpb/+iedSg48jC5E9/mIKRBNq/ub0/SOhmz4fTNRcXYatKuN0aR+tS8OI3gP4t+xn0U6IHUC0toe9DCF0YURPGvBcFu2Z9NiByFtXcC9CgxFwdIafzqM4lJQv/3XILPHL/DXMD3PXrocbr56ayZvL3gNgB6v38dnr34TiOU2GA3E5I+maacGgT6jH/mIv35cHQh93PTHVv437AHuee4OhBBE54sKe7B5831Nufm+ZlRrUolfz+7q2z/SmprNq/LT+PkXYcaDMZoM3HxWeMuZ5sJgMOAl9A1B12XgbSN/kTjadLvpEmdW4N0NwgTygjckYQA9lYjpyNrFKXACCFNlMKlIpZxEuWKuEXRdZ9PSrSyduZI9m/aHLWxxbM853ZnOfTvyzLhHKVmxGLEFY2j1QHOGzRnI0T0n8Pl8HNl9jD+/XxkUz+5Mc/H569/hdroRQtDt1XuC5A0MRo3YgjE0u9PvN615UzX6T36S/pOfpGbzqv5XdV3P0sFpJAwmAxXrleOJ9/w+2yKlC5G/aL6wbT1OD/1avsHSGSsuej7FBRgrgQyTKCZ9CGMFsN6B33US6ABaETA1QEoH0nda1TXNBagd+zXA6WMJ9Gv5BgknEhEInA6Xv2rReQbUYNSo3fKceJgQgtv+14bb/tcGR5qTEd0/4OnGr2Aw+TMv73rqNoxmY0g2qcflYe8/B6jauBJd+nfCFm3hu1GzSU1Io2H7ujzxXk9M5uCiGlJKZoyezbdv/xSyw88yAirVr8DzEx+jcoNzYmZCCF6f2Y+XbxmCz+sLSaxyOdx8+PRkbrq3SeaqlYpMEcZSSNvt56JeAIQNop5GaHaIG4LU4vwa69INllYQ0x+Z9BI4f/W3NxSDuDF+H3oGSOev/sNTPRVsdyCiHkUIdeCdE6jM02uA1+58m7ULNpwV3jqLAIvVjMvhxmIzY4u2MmHdqLAqiGN6f8zv05cFRbGYrSYkhFVmjM4XxYS171C8fGgGajhmfzSfyQOmZaotk69ILCkJaeF39Gfzkyw2M8XKFaVj7zZ06H1LILLFme7is1e/Zu6kRSFrttjMfLHzQwqVvHh3gMKPdG9Cps8E376zsgCa/0+ZDuamiJgXEYbi/rZ6mj9r1PUrIYJwIgpR+A+EFhd2Hj3tc0h5n3MRMRYw1UQU+FplnWZAVjNPlSsmlyOlDDXq+DVmmnVqyC0Pt+Tht+5nyrb3wxp1KWWIUQd/SGCBYvkCBvV80pPTmTxwWpbX+N2o2Zka9TLVSzHz+BSen9gnfAPp/8+V7ubA1kNMemkaTzcaiPPsGYHVbqHtgzdjMIT+ygpNEFsoNsvrVYRHT5vuzwR1zvQfhnr3gHc7+A74Bb2c85CnOiPPxpnLpIHg+p3wKp8SnAvDziOlF1I/JDjM0QXef/1ZropLRrlicjlCCIxmEz5vsAvCYDRQt01Nbu/TLtMxwqX/e1weTh85E/YsTNclWy6IRDm4/Qg/vD+XEwdO0axTQ257tDVmqz96ISUhY/eLxW7mpnuasODzJSyYsjjT9YI/Kib+4CkWfbWUW3u2JC0pnUr1y1Gpfnl2rNkdcCFZ7Ra6vnwXZkvWa64qQpF6OqSMIkRSIAgdpAPpmA22zn6xLyIIw0nv2bDIcJ8lh/fjI/wSBtnIclWERxn2a4DbHm3N/CmLgxJ6hCZocV/mJduEEDTr1JBVc9cG+eR1n0T3RT7kPD8scsvy7bxy2zDcTg+6T2fzsm0s+uoP3v9rmP8B06o6q+auj7AAiC0Yw4xRs8PWPs0IV7qLH96fx8T+XyKlJK5gDM9P7MP2NbtZPP0vbGezU29RmumXjm+/P/IlU8+sA7z7CZsxGoQAS+sIH+Xzx6zLCx4cUgdjjaysVpEJyhVzDdBn9MPc3KUZJosJs9VEkdKFGDl/ELEFYrLU//lP+lCm+g1YoyzYY20YzUZM1sg7XIvNzP+GdQv8/NFzn+FMcwV2/q50Fwf+PcyquesAaNi+XsSx7NE2kk6lRDbqwv+QCodm1Di27wRuhxuP08OpI2cYdv9YOvRqy5e7xjFxw2jaPdxS+WRzAKkVjbzDPh9hR1gaIbQCYChLWF8eZoh5CWEMX9vYn+X6Bv7omrP9hQ2styBMVS7uBhRBqB37NYDZYmLAF8/Qd1wv0pPTKViiQLaMWVyhWD5eP4o9G/eTeDKZfEVieO7G18K2rdSgPL1GdKf+LbUD1w5sPRzSzpHqZNf6vTS/uzGLpi6NOLfb5cGbwU69XI3SVGpQnj++Wx4UoSM0gfRJdBnsRvJ6fSyaupQy1Uqxeu46CpUqSIfebSlyQ6GIcygyR3h3IsMaaQAzfreMDYzVwNLW3yff+8gzDwKesxIDHrB1RUQ/hTBkXNdYs92ONN6ATJ8OejLCejtYO+bgHV3fKMN+DWGPsUWUr80MIQQV650Ta2reuQkrf14TCB+0Rlno8dq9PDCgc0jfYmULh2ShGs3+whi6ruNMDV8XVDNoxBSIIuF4UsR1Hdh2mAbt6/DOb2/w+7RlbPlrG45UF3Vb12DZ96tIT3EEtfd5fcz75DcS4hPxurwYTQZmjZ3L6MVvZqnmqyI80r0CwiSBgQWi+oDvIMLSAqwdEMJvNoSpEhRZBq5l/t2+pbl/J59FhKk2Iq525g0V2UaFO14neD1eZr77M/OnLAYBtz3ahrhCMSz6cilmu5m7nr6N/EXz4Xa4qX5jFTRNkJaUTkyBaNbM38DQru/hukC0yxptpVK9ctRvV5tvRv4YIup126NtuKFqSb5687sMo2YsdguPDLmf+/p1Cro+6pHxLPnmr6CzAaGJsAW4qzSqwPjVb1/MV6MAZNpXyJR3CT48BbQSaEX+uBpLUoQhq+GOyrBfJwzv/j4rZ68JGGeLzUzzzo15ZdpzHN1znJdvGULyGX8Ym9ftRdM0dF0nOn80z37Um5gC0Yx6ZLy/+PV5vzJWu4XHx/Tkrx9WsWX5DgwGDY/bQ5vuLej/6ZPoPp0xvT9m6YwVGEwGnKnOsIJgRUoXYvr+j4OuJZ9JoV/LwcQfOAnCH8nj8+phDbvBaGCB+9uc+8KuM6SegDx5y9nY9bPfr7BB9CtoUQ9c1bUpzqEMuyLAycOn6VnpGTyuYF+3yWLiqz3jGdRxBPu2HAxrMMEfrvjhihEM7TqGwzuPhXze9I4GDP15ILs37uPIzmNUqFeOUpWKB7U5czzBH744dSlzJiwMSTuPLRTDrPjPQsaWUrLlr+2cPHyaf5ZuZd6k38KusXCpgnx9cGKG34MiY6RnFzL5TfBsBC0fRD2FsHdXh9O5CFUaTxHgyO5jIQlOAEazgV3r9nJo59GIRh3A4/Iy95OFFCieP8SwawaNIqX9B5cV65ajYt3QogsABYrlp0Cx/HjdXn79/I+gQtMms5EW94QP3RRCUKuFv+rRiQPxmCzGkAcUAnoOvT/i+hVZQ5gqIQpOv9rLUOQAKtzxOmDnmt1hk5TcDjdlatyQqVi67tNJSUijx6B7sdiDJVXNVhN3PdMhYt+kU8n8Nu1P/vpxNW6nm5o3VaNz3w6Yzuqn26KtlK5eit5v98j0Pto91NKv834eQghufbgV7XtGiJlWKK5D1I79OmDlz+vCXrfYLZQoX5QqjSqybdWusEWuwR8x0/r+5tS/pTYDvuzL5AFTiT94mjLVS/HMuF6UrloybL/F3/7FmF4fYzBoCCEwmAy8u/hNeo3swe2Pt+PflTspUroQNW6skqXX/UIlCzJy/iDefXQCx/fHY7KYuK9/Jx4e3DXrX4ZCcR2gfOxXECklB7cdxmQxhdQXvZwMuHUI63/bHHJdaIKfEr7Ele7ijbtHsXfTfhAC3af7jbEm0HVJ6/tvpP+Up7Lsa92xdg9Th8zg73kbQnzpN1QtwWf/fnDJ9+RIdWC2mZWio+K6QvnYcxl7/znA63e+TfLpFKQuKVm5OMPmvBJWuCunqX9L7bCG3WQxcerIGUpXLcm4lSOIP3gSl8NNyUrF+XflTo7uPk6VRhUoU/2GLM+1aelWBt0+wh99E2bPcHzfSU4fS6Bg8fyXckvYom2kJKSyZv4GjGYjjTrUUzVOs4jUz4B7jb84hqmBOhzNg6gd+xXA5/XxQMk+JJ5MDlzTDBoV65Xlo7/fuezznzpymh5lnwrxs1ujrcyKnxIQ88oJnm48kJ1r90T83Gg28PzHfVg6cyVmq4lOT7anQbuMdbvDsXreOoZ2fQ/NqAECTQhG/voa1ZpUuoTV5330tOmQ8ra/ShI6aMUQBaYjDEry+FpA7dhzEVv+2o77Atlc3aezf8shTh4+fdl37YVKFqTzcx2ZO3FhIFHIYrfw+OiHctSoAxzeGb5OKvjfEAoUz8f4Zz8LZLyu/mU9JosJn1enXpuaPDOuF8XKFslwDme6i+Hd3g9JmBpy37tMP/AxmqZiAsIhvfv9Rh3XudJ3voPI5DcQ+T+6mktT5DDqX8AVINKhJISX1L0cPD76YV6f0Z+WXZrR7uGWvLPwde54/NZLGjPpVDI71+3BkXou7b9szfDCT5pBo17bWpw5nhhUBcnr9uFIceJ2uFkzfwN9m76Ky5GxcuC/K3aEFQ5LTUzjyK7QOPvrEelain6yHfrxqugn2yGdSyJop3vB9cdVWKHicqJ27FeAWjdXC9lFCiEoXr5okDzu5UQIQZOO9WnSsX6mbQ9sO8ynA6ezc81uSlUuwaMjulPjRr/qntvl4cDWQ8waO5c/Z63CZDbi8/h4dEQ37nnuDp4Y05OX2r7l15Hx6ZhtZoqULsSHK4ax5a8dbPlrW9iqTeDXgXc5XCz/aU2Ghantsbawcfc+r47tIrV08hLSvQmZ0JeAPIDvADLxObA/CBgI0VAXSss+r6EM+xXAZDYxfN4rvNbpbXxeH1KXxBaK4a2fXr7aSwsh/uBJnm36Ko5UB1LCmeOJDLh1CO8ufpPj++IZ2+cTPG5voCLTf0b6s0HfUrlBBWreVI2P/h7J92PmcGTPcZp0rE+nJ9tjj7FRulrJDJUewS8bcOZYQoZtqjSqSIHi+Tm290TgjUdoAqEJ3n10Ag+90SXwILoekWlTCNVLd4F3NwhxwaG2FWz3XbnFKa4IyrBfIao3q8KMY5PZvnoXJouJyg0r5MpohB/Hzcft9ATlLLnS3UweMI0df+8O8Wv/h9vhYv6UxdS8qRplqt9A/ylPhbQ5deQMegYZruDXfKnTKuNiC0IIRv8+mBE93ufflTvRvTpSStwON+sWbmLLsm0MnTOQem1qZX7DeRH9BKEhSRL0BET+T5GJ/UBPAnSwtkfE5L4NhuLSyDEfuxDCIITYIISYm1Nj5jWMJiM1b6pGlUYVc6VRBzi47TBeT+iuet/mg3gzOCuQEnQ98nlBWlIar3d6O2THrmkaZqsJIQQWu4V2D7eiUv3yma6zcKmCjF06lMfeedCfDXv+g8jh5tOBuSc1XkqJlJG/u6z1d4XkBAQ+92xHP/0g+vE66Cfbg1YGfxGL87GCtR3C3AhR+E9E4QWIIqvQ8o1GiJw9QFdcfXJyx/4csA1QVYWvYerfUptNS7YG7cwNRgNFyxZm76YDGfZNSUhDShn2obV6XoTSeRr0GtmDlIRUGrWvS/Vm2XOhHNl1PKwk8NE9x7M1zuVASolM+xjSPgWZhjRWR8QNR5iqZ30M52Jk8lv+XbjIh4x5Cc1+77nPfSeQZ7qBPFt31rcPfMdAKwLyFEgnCCsYyiLsDwL+Nx4MJXL0XhW5ixzZsQshSgG3A5/mxHiKq0fH3m0pUroQFrsFALPNTEz+KPqMegijOeN9wKYlW9j0x9awn/kiRP8IBLf3uYWeb96fbaMOUOPGKlijLCHXzy8qcrWQaVMg9ZNzUrjercgzDyL1jM8QAv09W5GJz4N+DH8h6TOQPATpOlexSqZ/F1o7FCdo0Yi4sYjo5xFx7yEKzkJo9hy7N0XuJqdcMe8DLxMaSxVACNFHCLFWCLH25MmTOTStIqexRduYsG4UT4zpSdsHW9DzrfuZsu196rWpxZPv9cRsMxMVF95AuNJdrPh5TdjPmtxeP8S//p8/3WILNcxZ5eYuzShRoVhAnMxgMmCLtvL4uw9f9Jg5RvpnQHAFKKQXHPOy1F2mfYW/JN35OJCpk8/96D1ESJQLgPcgwtoaEf0EwtoGIZT0wvXEJbtihBB3APFSynVCiFaR2kkpJwGTwJ95eqnzKi4fVruFOx5vxx2Ptwu6fsfjt9L6gebsWr+Pd3qO49ThM0GfG81G4gqGL7AdWyCG12f0Y3i3sYizejTFyxdl4NS+EdeRlpzOgim/s2PNHqo0qsBtj7YhKi4qqI3ZYuKDFcP59fPFrJm/kZKVi3N33w4UL1f0Iu8+B9FTwlx0I/WEiNVFg/ufIuxeSU8893djmdAAGABS0U92APv9CPtDyrBfZ1yypIAQYiTwEODFf2ITC/wgpXwwUp/rTVIgL/LzhAVMenlakK66NcrClH/fz7CwtNvpZvvfu4nOF0W5WqUjHiInn07hiXovkXw6BZfDjcVmJqZgNBPXjyau0LVxjKOf6QPuPwk2zjZEga8Q5sxlFGT6LGTyEIJ3/f4apFqM/4Gop38HyYOJ/LJsA2t7tHyjLu4mFLmKrEoKXLIrRkr5ipSylJSyLPAAsDgjo67IG3R6sj3dXumMPcaGZtAoWak4I34ZlKFRBzBbzdS+uTrla5fJMDLo+/fmkBifFDjEdTncJJ1M5vv3rp2gKxH3pl9oS0QBZsAK9i5ZMuoA2O4CcyN/iTosIOxgrISI6oX0HUc6F577LCIOcM5H+k5c6u0oriFUHHse4czxBIxmI7EFwrtCchohBD0G3Uu3VzrjcXkuyU8ejg2Lt+C5IDTS4/Ky4fdQlcrcijCUgMJnU/l98WBugjBVzXp/YYT8k/2l6jxbwVgezE2RKWMh/XMQZpC+s39KQgpRBwYyg+8wGHKBe0pxRchRwy6l/AP4IyfHVGTM4Z1HGdJljF98S0LtltV59Zvnr5iB1zQtx406QOmqJdm5dk+Qlo5m0ChdLXxRj9yKEGawRq4wlXl/AeZ6/v8A6VoF6WcPVf+LhpE+MDcB3QHe9YS4ZaQbjEr18npCiYBdw/i8Pvq3fpP9Ww7hcXnxuL1s+mMrw+4fm2nfZT+s5uFKz9DB0o0nG7zEluXbr8CKs879A+4OUZ40W03cP+Duq7Si3IF0zCJ0Z+4Gz0ZEwakQ3R+wwX/Hs8IG0c8gtGvjXEKRMyjDfg2zedk2nGnOoIxEr8fHlr+2kXgyKWK/dYs28c7DH3Jszwm8Hi+7N+xnYPthGUruXmlKVy3J+38NpXGHehQskZ9Gt9Vl7LKhlKlW6mov7aog9RT0M73AOYewFUxkKjJpEFr0Y4gCk8F6F1jvQuSbhBb9+BVfr+Lqogz7Ncz58rfnI4TAHUFBEeDrET+EZGt63R5+Gj8/R9d3IY40Jz9+OI/X73ybyQOmEX/oVIbtK9Qpy/B5r/LI0AfYs+kATzUcQJ+6L7J1xY7Lus7ciEx6FdyriRz9IsE5D+k7iTA3Rss3yi8XYGlyJZepyCUow34NU6dV9bB67kXKFM6weMepI2dCrvm8OvEHMja0l4LL4eKZxgOZ8urXrJq7jh8+mMtjtfpxYNvhDPstm7WK8X2ncOZYAlKX7PvnAAPbD+XYvusnykNK11kt9fACbAGECXy5561LcfVQhv0axhZt47Xv+mG1W7DH2rDH2ogrHMvg71/MMJSw6R0NQuQBrFEWbryr0WVb6+/TlhF/8FTgTcFfYMORqVjX9OGzwrxd+Jg3adFlW2uuI6sCYtILxgqXdy2KawIV7niN06Rjfb47NpkNv2/GbDNTr01NjKaM/7f2eO1eVv68lsT4JBypTmzRVirULUvbB1tctnVuXr49xHUkJez4e1eG/RLjQ88KvB4vp49mTW8lLyA0O9JUFzxhIl4COaxWiHkeoUVf4dUpciPKsOcB7DE2mt/dOMvtYwvE8OnWsSz/8W8O7zxK5YYVaNi+DgbD5Us7L1ezNGabGfcFeu4lK2esMligWL6wRrxem5o5ur7cjsg3BnnmQdDPAMK/O496AvSDgEDY7kOYG1ztZSpyCcqwX6eYLSZaP9D8is1326OtmTl6Nj6PF5/Xv+u02M30Gt4tw35HdofK72oGLcPD4byElG7AhzAUh0KL/MlKMgVMDRFaVKb9FdcnyrArrgixBWKYsG4UX701g3+W/kuJCkV5aHBXqjetnGG/cFrrQhOkJTnCtM47SD0NmfwaOBcCOtJUF5HvXYT5XM1aKd3+jFQtFqF864rzUIZdccUoXKog/Sc/ma0+dVvXYMPvm4Mkf41GA41uq5vTy8tVyKSXwPUnAUlezwbk6R5QeDFCaEjXn36tdgDpRRorIQpMRmgFrtqaFbkHFRWjyDb7ExPoPmsGFce9R52J43l/1Qp8GZTFi4TP5+PY3hOkJadHbPP8J4+Tr0gc9hgrZpsZs9XEvf07Ub52mUu5hVyN1JPAtZTg8EYdZBJ41iL1BGTCM/4CHjIVcIJ3GzLplau0YkVuQ+3YFdnC6fVw38xvSHQ60aUkxe1i8vo1+KRO/2Y3ZXmcFT+vYUzvj3E73Ph8Ou17tuKZcb0wGIMPcIuVLcLUfRNYPW89iScSqde2FqUyOXC95pHphN9zCdCTwbUkzOdecP2JlG5Vw1ShDLsieyzauweX14t+noyBw+vli43r6de0eZaKdB/dc5wR3d4Pqqu6aOpSCpYsyIOv3RvS3mwx0eKeazuDUuqJ/sxREQfmRhkXvtCKgaGwX5ExaBCPX+zLdR3F8CsuCuWKUWSL0+npeMO4XdI9niBjnxG/Tf0Trzc46caV7mbuxIU5ssbLgZQS6dmOdP+NlBHkcSOgp/+AjG+BTHoFmfgk8mRbZJgMUendh578NjLpZbD3AKL9Wu7CDlggbiRCiwFLG0L1YkxgaaV26wpA7dgVWeBYSgoT1q5m/bGjlIiJDTEpAqhdtBgGLWv7BJfDhe4NfTh4XLkzhFHqZ5Bn/ge+/YAB0JGxo9Fs7TLpib/ARfJgwAXybIKWdCATX0YUnHaunWsFMuEJ/IXIvP5dueVmhLW9v5+lJcLgl4kQWj7IPwGZ+Bx+37sXTNUQcSNz9sYV1yzKsCsy5GRaGrd/8xWpbjdeXWfn6VNoQmA2GJBSYjIYMGoa79zSPstj3nxfM2Z/tCAolNFkMdKya7PLcQuXjEwcCN5d+I3uWZL6I81LAsY2Iq6lhL4Y6+BZF+QPl8mvEyTHK9P9faMfR5hqhAwrLM2hyErwbgMRhzDm3cNkRfZRhv065mR6GuluD6Xj4hBCIKVk3bGjbDpxnLL58tGqTDm+2LSedLcn4H7xSYkEbipVmualS5PPaqNDxcpEmzN3AcSnpfLO8mX8eWA/2tAmmL7bQcHtKeg+nbI1S9NrZI/LfMfZR0oPuJcBF+i1CAGuxWDvkvEAwgZCC6O0q/GfwZfSBb4jYTrr4N4AYQy7fwkmMNXOwl0orjeUYb+GSHW7+XbLPyw/dIBKBQrRs249SsZkrYCClJINx4+x7tgR8ltt/LT9X9YeO4pBCPJbbYzrcAcfrVnNqiOH8Oo6Jk2jWHQMhaOicOvBRk2XksMpyTxWP+uiYS6vl7u/m87JtDR8UoIBrA9VpkNsObpUr0Gl+uWzdPB65RGc02M5Dynwu2UywdIGGHzhRbB29Je+A8AMIhpk8gVTG8FYOtsrViiUYb9GSPd4uPPbqRxPScXp87L0wH4+37iOwS3b8GDtjJN1dCl5fsE8ft+3F6/uQ5fSb1zP4khNocePMxH4I1wA3D4fexLOsCfhDJoQQQejGoI6RYtla/0L9uwixeUKmtfp8zLXfZRBDe7M1liXG+negHTOB2FF2DqD5ZazsrnnnQEICda2mY4ltCgo8KU/mcgXD0iwtPYXuv6vjRDI6L6Q8h7wX0atyR8dY75ysg+KvIMy7NcIP2zbyolUv1H/D5+UvLl0MUCGxn3pgX0s3r8Xhzfy4aTH5wsyuuejSxkw7kZNw2Y08mzj7PnDjyQnBx4a53MqPXJy0tVATx0HqZMBF2BApn0BccNBJvrdIsIAmBH5PkBocVkaU5hqQaHfQD8BIsof2XIBWlRPpFYQmfYJ6ElgbYeIfi7jsEiFIgLKsF8jrDl6JKxh1KVk9Iq/6FazdsSolMV795LuyTjiRAiBOOs/D4uUNChegpqFi9K7QcMsu4D+o2GJkliNxpB1VC9cJFvjXE6k7wSkfsK5jM+zESrJQxFFVoIe7ze6xornuVGyhhACDBm/5QjbHQjbHRe1doXifFQc+zVCtUKFMUYw3E6vlxR3+DJ5AIWjojBrGe/8DEJgyUC216BpTLmzM4Nbtcm2UQdoVKIkN91QBrvJBIBZM2A3mRjS+pZsj3XZ8GyCcHHg0gm+YwhDcYSparaNukJxpVG/odcI99eoxcR1f5PsCjXgMWYzsRZrxL5dqtdk0vo1QYegAjAbDJgMBny6zvvtb8duNvHq74s4lBxa3KJsvvwZzpEZQggm3H4ni/ftYcn+vRSJiqZr9VoUjwl1SzU6YgIAACAASURBVFw1DCUiVCuSoOW/4stRKC4WIbOYLZiTNGzYUK5du/aKz3utczw1hS4zv+VoSnLAZWIzGhnZ9lburFItw76bjh9j0OJF/HvqJAVsNvo2akqtosVIcbloVLJUYCcN8PnGdYxavgwhBAYhMGga39zTlWq5yG1yOZBSIs90Bc82zrljrGC/Dy32jau5NIUCACHEOillw0zbKcN+beHTdWbv2MbPO7aTz2rl4Tr1qF8866JYUspAWKHL6+VYagpFo6KxnWfYAQ4nJ/HH/n1Emy20K1+BqCzEqV8Kp9LTmbdrO2luD23LV6BKwUKXdb5ISD0VmTIWXL8AZrA/iIh6VB1iKnIFyrArMmTa/9u77/goq6yB4787z8ykF3pJpfceuoAIiLACFkTFgrqrrrurq9td9911i7vqFnHXfdfFgvvaCygoSBGkQ+g91EAagYSSPpmZ55n7/jEhm8lMwiQZmUm438+Hj5mSOycQz9y5z73n7N/LC5s3IHFfgH106HCeGDk6KHvJd5zJ5cElizFcEt1lYNE0nhgxmu+m+d/uT1GuBf4mdrXGfg05euE8R88XYjcM/rRpvccumwW7d9C1dWtm9Ox9VWOSUvLjlV967JYxdJ2X07cwq1ef0FqDV5RmQiX2a4DucvH4l5+zPus0mhDYapXdhf+W3r3aif2CzUZBRbnX/WaTie1ncpl1hWsHiqJ4U4n9GvDxoQNsyDpNpY998DX5ejw9N4c/b91EVlERgzt24hdjx9Gt9RUKXzVAtNXi68A+IGgXqZo1K0pjqMQeRLklxXxx7AgOw8VN3XvQ8xu6YPjZ0Qyfh5tqCjebmd3Xs9jUjjO5PLh0cXXCX3vqJNtys5nbfyCLjhzG5nRyfWoXnp0wiXZRjUvC4WYLd/TtzycZh6pfxywE7aIiGZWY1KgxFeVapxJ7kKzJPMnjK77AcLlwScmru7bzi7HjuX/QEIoqbWzJySHaamVMUnKdB5P8FWXxvaPFJASRFgu6y8W45BTuHeBZluDlbVs8ZvESKHc6eW3Prur7Vp88waHCAtbc95Df9djBXeN9VeZxTMLEY2kjiQ0L490D+6jUDSamduG3EydhCsmiYIoS+pqc2IUQScD/AR1w/7+/QEr5clPHbcl0l4ufrF7hkTQNXedPm9ZjNWn8dsPX1ck80mLhg9l30iW+8QdkHhw8lPS8HI9Ze4TZzO8mTibWGkbXVq18Lq9kF3sfVPL6WaTkQoWNzTnZjE9J9SueZceO8tPVK5C4t17+adN6Xpk2g5+MGef3z6QoSt0CUVJAB34spewLjAK+L4ToG4BxW6ysoks4Xd4nHC0mjV+vX4Pd0Cl3Oih3OjhfUc4Pv/yiSa83LiWVp6+bQLTVSrjZTLTVyi/Gjuf2Pv2Y0q17nWvmIxMT0fyYNUskBeVlfsViczr5+VcrqTR07IZBpa5Tqes8uWIZjissFzVEqd3Oewf28cKmDaw9lel32z5FaQmaPGOXUuYD+VVflwohMoAE4HBTx26pWkdE+uwb6jAMLCaTx2MS9zbFEntlk4703ztwMHf2G8AFWwVtIiKx1FMX5rKnRo3l69OnKHM4cBi+jtq7GS4XIxP8Ww8/UHAOk8n7zaLM6eC6t17jpanTGZvkfzcgm9PJm3t3seLEcVpHRPLIsDRS41txywfvUOF0YtN13j6wl2GdOvPmzNsatFykKM1VQH/LhRCpwBAg3cdjjwghdgohdhYWFgbyZZudVhERTOvWgzDtv++rVk0jITYWTfj+J7FcoYiXPyyaRoeoaA4VFvDhwf3sPZtPfQfUOsfE8tV9DxJ/hTeUh4cOJynOvxK2rSMifL6pgfv06SOff0ZeSYnPx2tzScndiz/ile3pHCosYGP2aR75/DMeX/45l2yV1UtPFU4nu/LPsCrzhF/jKkpzF7DELoSIBhYBT0pZuxUMSCkXSCnTpJRp7dq1C9TLNlsvTLmJBwYPIT48nGirlVt69eHD2+/ErHn+k1g1jSldu3sd+W8Mh2Ew77NF3LP4Y3634WvuWfwx8z5bVO9s3KqZuWCru2Z6q7BwfjTa/2YQ3Vu3oXfbdljqmDnrLheLjxzya6wtOdmcvHgBe40a9TZdZ9+5s7hqFSCucDrZlH3a7zgVpTkLSGIXQlhwJ/V3pZSLAzFmS2fVNH4+djy7H/k++7/7OM9PnkrbqCjev/1Ourdu4668aNKYmNqlQY2i6/PBwf3sys/DpruXKGy6k135eXx4cH+9cda3K6fEYae4srLOx31ZOPM2pnbr4fMx3eXyWcHSlxMXL+D0Mfv39RkkTNNIiPHvU4WiNHeB2BUjgDeADCnl35oe0rWtV5u2rLr3AQorygnXzMSEhQVk3GMXzvPSti1e+9ltus7SY0e4b9AQn99nNpm4s98A3ju43+cSSpjZzJnSEuLC/V//jwsP5+/TbibCbObToxke44abzXUm/dou16iv/YmjVXg4tqqLspdZNI07+vX3O0ZFac4CMWMfC9wH3CCE2Fv1Z3oAxr2mtYuMClhSzy4u4raP3qPY7ntmHR8eUe/3PzPueu7o6zspSilJaeRWzKfHTSA1vhVRFguRFgthmsa9AweT1jnBr+8fkZDIwPYdiTC75ycC9xvD85Om8uLkqaTGxRNptjA6MYlP7rhbnWRVrhmB2BWzCZ9t3JXGOldWxqu7trM7/wz92rXnu2kjSI6Lb/R4r+/eWedWQgGMS06u9/stmsZzN0whzhrGW/v2UGno1Un0J2PGedRyb4j48AhW3DOP9Nwc8stKGdYpgZR4/39OIQRv3XI7Hx8+yPLjR2kbGcVDg4cyqGMnAG6+ynVvFCVUqLK9IaawvJyb3n2LUocD3eVCE4IIi4Wld91LaiNnxnMXf8S23Jw6Hw83m1k8Zy6929Z/UVtKyYas03yccRCLycTd/QcxIiERgBJ7JacuXSI5Lp5WEfV/Amio4spKMs4X0ik6pkGJX1FaGlW2t5lauHcXZVVJHcCQkgqnk1e2b+MvN05r1JhjEpPYk38Gex27XxyGwb937eClqfWvoAkhmJDahQmpXarvk1IyP30LC3btwKJpOA2DewYM4plx1wektvs/d6TzcvqW6o+EY5NS+Ne3ZhJmVr+6ilIXdVojxOw7d9Zrp4dLSg4WnGv0mPcPGkK7qCjC60iGLil99jn1x9pTmby+eyd2w6DM4cBuGLx/cD9Ljh5pdLyXrck8yV+3bkJ3uXBW/VmfdYr56VuaPLaitGQqsYeYwR07Ya11GMkkBAM6dGj0mLFh4SyfO4+fjRlHj9atvcoEhGka1/tZ56W29w7u87nT5t0DexsbbrU/bVrvdZ/EvW1TUZS6qcQeYh4YPJToMGv1AR5zVQXGHwwf3aRxo61WHhg8lI/vuJvE2LjqC56RFgudY2KZN2hoo8at6xJNIGqznC3zXX+mvgNViqKoNfaQ0y4yiuVz7+e13TvZVbUr5tFhw0mMbdzhmkrdyR82rOPTI4fRXS6mdO3Oe7fNYUtuNocLCujXvj3f6tGr0WvWd/YfQHpeLjb9v63tIsxm5g4Y1KjxakqKjePoxfNe9w/q0KnJYytKS6Z2xbRwjy1bwrrTp6ovnFpMJnq0acvnd91LqcPOhwcPsOdsPoM6dOSu/gMbdNAI3BdPn9+8gf/btwerpmE3DO7o25/fXt/weupLjhzmpfQtnK+oIK1TAtN79ORXa1ej1/odTevUmYldunLPgMHE1rPXf8WJ4/xl60bOlpUxuEMnfj1h4jfWzERRrgZ/d8WoxN6CXaioYOzCBV5LF5FmC6/NuIUfr/6SIlsllYZefcp12dz7aRsZ2ajXOnnpIl3iWzWqm9KSI4f55drV1ev1AogNC+N310/mlR3byCkuwuFyIXDvFArTzHSIjmLZ3fcTZfVuJPJV5gmeWLHM4/RptMXKmvsfanS3J0UJNn8Tu1pjb8GKKm0+67xoJsH7B/dz0WajsqqAVqWhU1Rp49+7tjfqtdpERjIiIbHRSXN++laPi7ASsBsGF20VrLz3Aa5LTkVKiVE1EbEbOoXl5XyS4btg2Pxa3Z8AnC6DDw+pC69Ky6cSezNRYrfzwcH9vLJ9G3vyz9RbbveyTtExGD7quzgNFzklxV4zeafLRXo9B5m+SecrvCtIVuo6eaXuQqEHCs55Ffey6Tq78/N8jldQXu51n90wyPGzJLCiNGcqsTcDJy5eYMJbr/H7DeuYv20L9376Mb/6+isAcoqL+TTjMFtzsr12ojy7fq3P3Sm/mTCRPm3beW171ISgR4DWoO26zvLjx3hzzy6/9uCPSEjwWpOPtFgYU9V0o3vr1l7fE66Z6du2vc/xRiUm+RxvfHKqnz+BojRfaldMM/CrtV9RYrdXz1htus5nRw5j13WWHT+KZjIhgE4xsXx4+520ioigwulk6bEjXoedIsxmkuLiGZmYxNJjR7A5nUjca9phZjOPpY1ocrznysq49aN3KbXbcRouNJNgZs8+/HHSlDpPo/7P+Ins+fA97IaOTdeJtLiLd13uo/qTMePYvehDKnUdibvqZJTVyp39B/gc7+fXjWdrbnZ1F6VIi4WB7Tsytbt/lSMVpTlTF09DVJnDwYubN7Ds+FEu+ah3bhICkxAeJW8tJhMze/Xmz1OmUVhRzriFr3ktt0RZrPxx0hRm9OzNkcICXtiykWMXLtA2MpIB7TswKjGJG7v1wOpH67y6/ODLz1l54nj1ejhAhNnCGzNvZVRi3S30SuyVfHYkg7zSEsYkJjMuJdVj1n2o4Byv7NjGyUsXGZWQxA9GjKJ9VHSd45Xa7Sw9doSsokuMSEhkYmpX1RpPadbUrphmbvbH73Pw3DkcPppeg3vG6qs+emxYGHsf/QFSSq7/zxtepQLCNI0NDzzM+qxTPLdxHRVOJ4bLhdmk4XAZRFkstI+K5tM75za6x+qgV/9BqcPhcZ8AvjMkjafHTWjUmL7sO5vPipPHibJYubV3XxJiYwM2tqKEIrUrphk7XFhARmFhnUld4E7sZh/LGhFm94lSIQTzp04nymIhwmzBqmmEVXVtyrx0kd+sW0Ox3Y7T5cIF1a9V7nSSV1rCgl07Gh2/rzcEq6bRphHbKOvy0rYtzF38EQt27eAf27cy5Z2FbMw6HbDxFaU5U2vsISi/rBTN5Hst2mwyMaJzIk+MHMWDSxaj19jSF2E2M69GJ6QhnTqz6cFH+PLEMSqcTm7o0pXU+FY8tmyJV32XmhyGwdpTmczu25+dZ/LoGB3DmKTkeg8c7TiTy6s7t3OuvJzurVtzrqzU42CR2WTi1t59G/LXUCf3G8/26kNXlwuE/fSrFWx56NErHoyyV/3sqkKk0lKp3+wQNKhDJ5w+6qHEhYWz4+HHqvemvzbjVn66egUXbBWYhODeAYN5ZNhwz+8JD+eu/gM97vOnp2iZw8G0d/9TfWG2Y3QMH82+y2et9TWZJ3l8xRfV+8YF3n1Hnxo11uce91K7nY3ZpxFCMD451edho9p255/BbDJ5lSEurqzkfEV5nevu5ysq+MnqL9mcnYUQgnHJKfxlyrSA149XlGBTiT0EtY2M5MmRY3h5+1Ycuo7JZMJsMvHC5Bs9DhyNSUpm04MPc9FmI9pq9XsGOqtXH/aeza9z1m7VNAoryt2Jsyp5ZhcX8cdN6/jzFO+a8H/YuM7jMJCvqzbPbVxHSlw8k7p2q75vc04Wjyz9DBdguFwIAa9+axYTu3StN/6O0XVdMBX1lhh4cMkijl44776oKyWbsrP4zuefsmjO3HpfT1GaG7XGHqIeTRvBh7Pv4uFhw/le2kiWz72fG300eRZC0CYyskHLCrf16cd1yamEm81EWdzr71EWC1aTRtdWrZiY2sXn4aVVJ09QqTuRUrItN4c/blzHv3dt96uWuwS+u2wJS49mIKXEYRg89sVSbIaO3dDRpXs55Tuff0rmxYv1jpXWKYHE2LjqCpjgXoa6u/8Aws2+2/Qdv3CBzEsXPS44O10uMs4XklVUdMX4FaU5UTP2IDpyvpBzZWUM6tjRZ0PpAe07MKB94+uw18VsMvHvm2eRUVhAxvlCerRp6/E6/9y+jbWnMr32wJc6HAz+9z/pHBPLudJSKg0dq6b5dQoW3DVefrp6BV+fPsU9Awb53NUjgWe+Xs2fp9zE02tWsS03h0izhVGJSZwrL0N3uZjdtz/v3HoHL27ZyMqTxwk3m5k3cIjXMlRNJY5Kn1sdNWGixHHlpSlFaU7UdscgKLXbeXDpYjIKCzCbTDgMg5+PHc8DgxtXEz3Q8kpLuPHthfVeYK1N4P704E8d9jBN44XJU/nxqi899rpfFm2xEmGxcMFW4XO8CLOZG7p05R/TZvgdn13XGf76vyirtQ2z9nULRQllartjCPvTpvUcPHcWm65TWtVO7sUtGzl6wbv2eDAkxMTy+oxbSYiJ9So7UBcJ3NyjFyMTEvnZmHGkxtXddNpRVdyrlY9PKeDei1/hdNb5JmHTdb7KPNmgJZQws5lXps0gwmwm0mKp3gb6yvSbVVJXWhy1FBMEy44fw1FrGcJpGHx5/Ci9QqRe+OikZFbeM4+X0rfw9v69V+xaFGsN46Wp06tLBswdMJDbPnyPzKJLXs+VuD+1vH3LbG7+4G2PWXu4Zub61C4sPnK43tezahqnii6REl/3G0ht41NS2fLQo6w9lYkQcEOXro0+hKUooUwl9iDwNUM0CYFFC51/jtd27WB++pbqRtL1Cdc0bujalfs++wTdMJjTbwC39O7Lj0Zfx1Mrl/n8/ld37WDvubOsuGcef9i4jr1n8+kYHcOPR4+la6vWLMqoP7E7DIPebRv+JhgXHs6tfQKzn15RQlXoZJJryJ39BvDWvt0eWwQ1k4mZPXsHMar/2pB1mvnpW+pcYw/XzIxJSmbHmTxiw8Lo0boNK08cr37+gYJzbM/L5dnrbyDKaqXIR60bm66TnpfLocICFs663evxh4em8fqenRguF1JKj8NOEWYLd/YbQMfomAD9xIrSsqjEHgRPjhpDfmkJK04eR6vao/7i5KkkxTWur2mgvb1/b70XTgUws1dvShx27LrOppwsjx0uNl1nydEMnhg5mg9uv5MnVy7nyPlCr3EqnE5WnDjOzF59vB770eixTO/Rk/VZp4gPC6dtZBSfHT2M3TCY3acfU7p2D8jPqigtkUrsQWDVNF666Vv82mbjgq2ClLh4LE2ophhoDqP+3TAGkqfXrPZoYF2bVdM4efEi41JSWT73fmZ98A4HatVl14Sot+NS77bt6N22XfXtmoebLisoL2PNqUzMJhM3du3e4J6titISqe0AQdQqIoLurduEVFIHmN23f3UxsdrCNA3dMOpN6uDuVtS9dZvq24+PGEVErUNUVk3jvoGDGx3n58eOMOGt1/nDhq/57bq1jF24gO15uY0eT1FaCpXYFS839+jFLb37EKZp1clYEwKLyURa54Q6k/5lJty7ZGZ//B7j33qN+du2MC45lWcn3EDbCHeFx9T4eBbMuMUj+TdEmcPBz79aid0wsOk6FbqTCqeTx7/8wq+99IrSkqmlmBBzrqyM+elb2JqbQ0pcHD8cOYahnTpf1RiEEDx3wxS+N3wkxy6cJyU2no8OH+D9gwfYnJPNlXa2u4Dztv/2MH1153YOFJzljZm3cUe/ARguV5MbXlwuBFZbmcPB6aJLdG3l3UpPUa4VKrGHkBJ7JTe//zbFlTZ0KckuLmLHmTzemnU7IxISGz1uYUU5Xx4/hsMwuLFbd5J9HB7SXS42ZJ3mXHkZIzon0K11GxJiYkmIieX9g/s9Lqg2dD7scBlsyckm89JFurZq7ZHUXVKy80weBeVlDO+cSIc6C3x5igsP9zkzN6SLmHoKgSnKtSAgiV0IcRPwMqABr0spnw/EuNeaTw4fotzp8NjaV6nrvLh5I5/MubtBY0kp2XfuLMuOHeGdA/vcx/1dkr9u3cyzE27w6BVaWFHO7I/e52LVEX4J3N1/IP8zfiLgnnHX3iWjCYGseh1/Er1LSnKKiz1m0hcqKrhr0YecLStFIHC4DJ4aOYZH/ei7OrB9BzrHxHD60qXqvy+rpjE2KZl2kXVfkK1Pbkkxe8/mkxQbx8AOHevsz6oooa7JiV0IoQH/BKYAucAOIcRSKWX9J0wUL0cvnPfY237Z6WLv05v1cRoG3/n8U3bk5VHpY4fLs+vXMq1Hj+pTl79f/zX5pSUebygfHNzP9B49GdYpgWK79z50Q0ru6jeA9LwcThcVXTG5O10u+tUqaPbs+jVkFRd5bJV8eftWJqR28dgN44sQgndvm8PPVq9kU/ZpTEIwrXtP/jBx8hUi8Sal5E+b1vP2/r2YTSYk0LN1G96+9Q6/6sMrSqgJxMXTEcAJKWWmlNIBfADMCsC415y0zglE1rowKYCB7Ts2aJwPDx1g5xnfSR3AopnYlX+m+va6rFMeSR3cnxTWZJ4EYEJKF6+aMZFmC1O79SC7uNivGXtcWBhta7XGW3sq06vCo9MwWHPqpB8jQrvIKBbOuo2M7z/J4e/9kJemTm9UIt6am8O7B/ZjNwzKne6LsIfPFzI/fUuDx1KUUBCIxJ4A5NS4nVt1nwchxCNCiJ1CiJ2Fhd6HVRSY2bM3iXFxhFftRLGYTERarPziuvENGmfp0Yx6Dxi5pKRtjeWKKIt3MrRqWnUp4WfGTaBDVDRRFgsWk4kIs5lpPXowNinZr+UKi8lE91ZteGnbZvJKSmq8hvcHRrPJ5LXrJq+khLs++YDe/5zPyNf+xRu7d3qUCjabTE26GLvs+FGv7ZsOw+CLY0caPaaiBNNVu3gqpVwALAB32d6r9brNSZjZzOI5c/no0AE252TRNb419w8eQkJMbIPG8VXb/TKLyURKXDz927Wvvu/bQ4bx0rbNHm8GmsnELb3dJ0LbR0Wzdt63WXPqJGdKSxneOYGBHdyfIjpERZFXWur1OpffnJyGgZSSXWfPsL/gLK/v3snCqovBcwcMZOFez9IKJiGYUaO0QkZhATM+eKf6QmmhrYK/bN1EpaHz/eGjGvT3Au7Trq/v3sny40eJj4jg0WHDibJY0ITwKiFcV9MORQl1gUjseUBSjduJVfcpjRBpsfDA4KEetdkzL13kgq2Cge07+tUp6aEhw9ick+U1azebTEzp2p3fTZzkMdN+aMgwSh123tizC7uukxwXzwuTp3r0DrVqGtO69/QY75LNRr6PpA5gNWnMv+lbPLZsiVfT6V+sWcXa+x/iyZFjKKqsZFHGIQTQOiKCv944vfo0qpSSeUsWee1+sRsGr+7czvfSRjboAqdLSuYu/oij5wurYzpw7iyPpY3AomkYtRqDPxgi9fEVpaECkdh3AD2EEF1wJ/S7ANVEMgCKKyt5aOlijpwvRDOZkFLylynTmNrdu0VeTaMSk/j9xMk8t3E9ZQ470VYrv7xuArf37e/z+SYheGrUWJ4YMRq7YRBp8W+meq68DM2k4XJ5l/QtdzqqWuZ5J97s4iLsuk6Y2cxzN0zhmXHXU2q30z4qyiNRnyq6RLGPAmLu8Z0YUmJuQGLfmpvNiYsXPJpg23SdhXv3MP/G6Tzz9WrKne4lmfsHDWnSqVhFCaYmJ3YppS6E+AGwEvd2xzellIeaHNk1yCUl/7sjnYV7d1HudBIfFs7FSpvHBcanVi5nXadve8ymfbmtTz9u6d23KrGHYfIjAWomE5ENWKvuEt8KTYCv4gLDOnUmISYWs0lgr5X3L/dZvSzSYvH5ZqIJU50z8tT4+AY3yMi8dAnD5b0KWGyvZHxqKtu+/V0KK8qJDw9XyzBKsxaQkgJSyuVSyp5Sym5SyucCMea16G9bN/GvnelcqqzEYRgUVJR79wUVsOrkCb/GMwlBbFi4X0m9McLMZn4/cbLXjpkIs5nnJ09lfEoqbSOjvJpOf2+4f0soKfHxpMa38przC+ClG6c3ON5+7dr7/LtoHxVFmGZGM5noGB1TZ1LPLy1l79l8Kpz118lRlGBTJ09DhEtKFu7dc8U+owK+sUTdGLf37U/fdu15ddcO8kpKGJOUzHfTRlTPwBfNuZsXN2/kq1MniQ0L49FhI5hTx5KQL2/MuJVvL13M6eIiXFISbbWy4OZZDOrYqcGxDunYiZGJiaTn5mDTdUxCYNU0fnf9pHrfaOy6zlMrl/P16Uz3WrzLxf+Mn8hd/Qc2OAZFuRpUM+sQYdd1+v3r71csYBVuNrPhgYe99oS3dKeKLmG4XHRr1bpJJ0J1l4vPjhxm2fGjtImIZN7goQyodXCqtr9u3cTru3dhr3EuINxsZtGcufS5wkEqRQkkf5tZqxl7iAgzm+nRuo1XQ2sBWDQNq0nDJAQv3TT9mkvq4F7PDwSzycTsvv2Z3YBPDR8fPuiR1MG9z33JkcP0uW5CQOJSlEBSiT2EvDjlJu5Z9BG6dOEwDKyaxgODhjB3wCAuVVbSu03bkKvdfi1QVYCV5kYl9hAyoH0H1j/wHb44fpQSeyUTUrrQv2qZIDE2NNrmXYtm9+3Hm3t2eWyTtGoas3y09FOUUKASe4hpFRER1P3TRZU2yp1OOkfHqOqGVZ4YMZoTFy+wIes0Vk3D6XLxzLjr6VPj9K6ihBKV2BUAyh0OfrRqOeuz3JUS20RE8sr0GQzq0LACZC1RmNnMv2++hbzSEvJLS+ndth3RquqjEsJUazwFgF+tXc36rNM4DINKXSevtIT7Pv2Ycocj2KFRWF7Ov3dt548b17E1J5tg7OQCSIiJJa1zgkrqSshTM3YFp2Gw/MQxnLUOQ0kpWXPqJDODuJa879xZ7l38EbrLhd0wePfAfmb07MXzk6cGLSZFCXVqxq7gktLn/nkp8dn4w/M5ksOFBWzMPv2NzO5/uWYV5U5n9YVLm+7k82NHOFhwLuCvpSgthZqxK4SZzaR1TmDnmTyP0rWGdDExtWud33fRVsF9n35CXQ0TqAAACaRJREFUVnERmhDoLhfPT5rKjF696/yehnBJScZ579r9hkuyPS+3eseQoiie1IxdAeAvN06jU0wMURYL0VYr4Zq73svlErq+PL1mFccvXqDC6aTU4cCm6/zsqxWcLfNdyvcyf9fIBRDrozG1RXPXdFEUxTc1Y1cA94XBdfO+w84zeRRXVjIyMbG6J6ovLilZeyrTqzmFEILVmSd9btk8duE8T69Zxd6z+cSEhfHw0LR6a6oLIfhe2kheTt9SXUNHE4KYsDAmd+3WhJ9WUVo2ldiVaiYhGJGQ6NdzLxcjq53YkdLnUc0Su505n3xAqd2OrLr9vzvSMQsTj6aNqPN1Hh6aRpTFwqu7dlBitzM+JYVfjZvoUfZXURRPKrErjSKEYGbP3nx+/CiOGicyKw2D5zat5/jFC/x6wg3VNdO/PH4U3XB5NL626ToLdu+sN7ELIbhn4GDuUU0vFMVvao1dabRnr5/E2KRkrLUaXjgMg08yDvG3rZur7ztvs3kV0gIoddi/8TgV5VqjErvSaFFWK2/MvI0nRo72Whqp1HXeO7iv+va45BSv55iEYKSfSz/BsCk7i8n/9ybd/v5XrntzAcuPHw12SIriF5XYlSYTQmDU7vSE5x74gR06MqffAMLNZiwmjSiLhfjwcP4wccrVDNVvR84X8ugXn5FZdAkJnCkr5SerV7A5JyvYoSnKFak1dqXJbujSjX9s34ZRI5FrQjAhpYvH834z4Qbu6NufLTnZtI2MYmq37kT42Tj7antr726PawfgfqP61450xialBCkqRfGPSuxKk/Vq05bHh4/i79u3Vl0sFbSJjOD3Eyd5Pbdvu/b0bQZVEc+Wl3nv+AEKKsqDEI2iNIxK7EpAPDZ8JLf07kt6Xi5tIyMZnZiEZmq+K31Tu3ZnR16uRw9aq6YxtWuPIEalKP5pvv/nKSGnU0wMt/Tuw3XJKc06qYO7SfeADh2JtFjQhCDKYiElLp5Hhg0PdmiKckVqxq6EDMPl3uduDoE3Baum8d5tc9ick8WhggK6tWrNxC5dQyI2RbkSldiVoCt3OPjV11+x/PhRXFIyNjmFP0++qd46NVeDSQjGJacyLjk1qHEoSkOp6YcSdE+s+IIvq+rBG1KyOTuLexZ/FLSGGorS3KnErgTV+YoKNudke2wtNKTkTFkp+1XNdUVpFLUUowTcmdISFmUc4qLNxqQu3RiblFxnBccSeyWaj8dMQlBks33ToSpKi6QSuxJQO8/kMe+zReguA6fLxceHDjKjZ2/+NPlGn89PjW9FjDXMY1shgO5yMaxzwtUIWVFaHLUUowTU02tWYdOd1f1TK3QnS45lkFFY4PP5JiH4+7SbibRYiLRYiDCbCdPMvDBpqmoarSiNpGbsSsDYdZ1TRZe87pcSdpzJo08dJ05HJCSy5aFHWJ15EodhMLlLt6DviFGU5kwldiVgrJq7uFdprabWZpMgISa23u+NDQvn9j79vsnwFOWa0aSlGCHEn4UQR4QQ+4UQnwoh4gMVmNL8XG5lF2H+73zBbDLROiKSCald6vlORVECqalr7KuB/lLKgcAx4Ommh6Q0Z48MG86vxk8kOS6O+PBwZvXqzaI5c9WJTUW5ipq0FCOlXFXj5jZgdtPCUZo7IQR39x/I3f0HBjsURblmBXIa9RDwZV0PCiEeEULsFELsLCwsDODLKoqiKDVdccYuhPgK6OjjoWeklEuqnvMMoAPv1jWOlHIBsAAgLS1NnRVXFEX5hlwxsUspJ9f3uBDiAeBmYJJUxT0URVGCrklr7EKIm4CfAROklBWBCUlRFEVpiqausb8CxACrhRB7hRCvBiAmRVEUpQmauiume6ACURRFUQJDBGNZXAhRCGQFYKi2wPkAjPNNUfE1jYqvaVR8TROK8aVIKdtd6UlBSeyBIoTYKaVMC3YcdVHxNY2Kr2lUfE0T6vHVRx0HVBRFaWFUYlcURWlhmntiXxDsAK5Axdc0Kr6mUfE1TajHV6dmvcauKIqieGvuM3ZFURSllhaT2IUQPxZCSCFE22DHUlOo1qwXQtwkhDgqhDghhPhFsOOpSQiRJIT4WghxWAhxSAjxw2DH5IsQQhNC7BFCfBHsWGoTQsQLIT6p+t3LEEKMDnZMNQkhnqr6tz0ohHhfCBEe5HjeFEIUCCEO1rivtRBitRDieNV/WwUzxoZoEYldCJEE3AhkBzsWH0KuZr0QQgP+CUwD+gJ3CyH6BjcqDzrwYyllX2AU8P0Qi++yHwIZwQ6iDi8DK6SUvYFBhFCcQogE4AkgTUrZH9CAu4IbFW8BN9W67xfAGillD2BN1e1moUUkduAl3DVrQu6CgZRylZRSr7q5DUgMZjxVRgAnpJSZUkoH8AEwK8gxVZNS5kspd1d9XYo7KSUENypPQohE4FvA68GOpTYhRBwwHngDQErpkFIWBTcqL2YgQghhBiKBM8EMRkq5AbhY6+5ZwH+qvv4PcMtVDaoJmn1iF0LMAvKklPuCHYsf6q1ZfxUlADk1bucSYonzMiFEKjAESA9uJF7m455MuIIdiA9dgEJgYdVS0etCiJDpDi6lzAP+gvsTdj5QXKtpT6joIKXMr/r6LNAhmME0RLNI7EKIr6rW4mr/mQX8Evh1CMd3+TlXrFmveBJCRAOLgCellCXBjucyIcTNQIGUclewY6mDGRgK/EtKOQQoJ4SWEarWqmfhfgPqDEQJIe4NblT1qypJHnIrAnVpUhGwq6WumvBCiAG4fzn2CSHAvcyxWwgxQkp5NtjxXRaCNevzgKQatxOr7gsZQggL7qT+rpRycbDjqWUsMFMIMR0IB2KFEO9IKUMlOeUCuVLKy59yPiGEEjswGTglpSwEEEIsBsYA7wQ1Km/nhBCdpJT5QohOQEGwA/JXs5ix10VKeUBK2V5KmSqlTMX9Cz30aib1K6lRs35mCNWs3wH0EEJ0EUJYcV+4WhrkmKoJ97v0G0CGlPJvwY6nNinl01LKxKrfubuAtSGU1Kn6/c8RQvSqumsScDiIIdWWDYwSQkRW/VtPIoQu7tawFJhX9fU8YEkQY2mQZjFjb+ZeAcJw16wH2Cal/G4wA5JS6kKIHwArce9IeFNKeSiYMdUyFrgPOCCE2Ft13y+llMuDGFNz8zjwbtUbdybwYJDjqSalTBdCfALsxr08uYcgn/IUQrwPXA+0FULkAr8Bngc+EkJ8G3c12jnBi7Bh1MlTRVGUFqZZL8UoiqIo3lRiVxRFaWFUYlcURWlhVGJXFEVpYVRiVxRFaWFUYlcURWlhVGJXFEVpYVRiVxRFaWH+H91kMyIqq3bDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gamma = E_step(X, best_pi, best_mu, best_sigma)\n",
    "labels = gamma.argmax(1)\n",
    "plt.scatter(X[:, 0], X[:, 1], c=labels, s=30)\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Authorization & Submission\n",
    "To submit assignment parts to Cousera platform, please, enter your e-mail and your token into variables below. You can generate the token on this programming assignment page. <b>Note:</b> Token expires 30 minutes after generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "You want to submit these numbers:\n",
      "Task Task 1 (E-step): 0.5337178741081263\n",
      "Task Task 2 (M-step: mu): 2.8993918820503812\n",
      "Task Task 2 (M-step: sigma): 5.977105216897527\n",
      "Task Task 2 (M-step: pi): 0.5507624459218776\n",
      "Task Task 3 (VLB): -1213.973464306019\n",
      "Task Task 4 (EM): -1064.9758118895377\n"
     ]
    }
   ],
   "source": [
    "STUDENT_EMAIL = \"goutham7r@gmail.com\"\n",
    "STUDENT_TOKEN = \"yo9ffuvfA1lOkk2N\"\n",
    "grader.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to submit these answers, run cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitted to Coursera platform. See results on assignment page!\n"
     ]
    }
   ],
   "source": [
    "grader.submit(\"goutham7r@gmail.com\", \"oLKeOUdLHBd2oQ6x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}