IRIS data set.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\"> IRIS data set </h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the IRIS data set from: <p></p>\n",
    " https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data <p></p>\n",
    " \n",
    " This is a data set of 150 points in R<sup>4</sup>, with three classes; refer to the website for more details of the features and classes. <p></p>\n",
    " \n",
    " <b>(a)</b> Use a PCA projection to 2d to visualize the entire data set. You should plot different classes using different colors/shapes. Do the classes seem well-separated from each other? <p></p>\n",
    " \n",
    " <b>(b)</b> Now build a classifier for this data set, based on a generative model. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "from sklearn.decomposition import PCA\n",
    "from scipy.stats import multivariate_normal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing iris dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2016-03-13 22:30:19--  https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data\n",
      "Resolving archive.ics.uci.edu (archive.ics.uci.edu)... 128.195.10.249\n",
      "Connecting to archive.ics.uci.edu (archive.ics.uci.edu)|128.195.10.249|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 4551 (4.4K) [text/plain]\n",
      "Saving to: ‘iris.txt’\n",
      "\n",
      "iris.txt            100%[=====================>]   4.44K  --.-KB/s   in 0s     \n",
      "\n",
      "2016-03-13 22:30:19 (72.3 MB/s) - ‘iris.txt’ saved [4551/4551]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data -O iris.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (a) Read Data and Perform PCA  to 2 components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal_len</th>\n",
       "      <th>sepal_wid</th>\n",
       "      <th>petal_len</th>\n",
       "      <th>petal wid</th>\n",
       "      <th>class</th>\n",
       "      <th>P1</th>\n",
       "      <th>P2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>Iris-setosa</td>\n",
       "      <td>-2.684207</td>\n",
       "      <td>-0.326607</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>Iris-setosa</td>\n",
       "      <td>-2.715391</td>\n",
       "      <td>0.169557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>Iris-setosa</td>\n",
       "      <td>-2.889820</td>\n",
       "      <td>0.137346</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>Iris-setosa</td>\n",
       "      <td>-2.746437</td>\n",
       "      <td>0.311124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>Iris-setosa</td>\n",
       "      <td>-2.728593</td>\n",
       "      <td>-0.333925</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal_len  sepal_wid  petal_len  petal wid        class        P1        P2\n",
       "0        5.1        3.5        1.4        0.2  Iris-setosa -2.684207 -0.326607\n",
       "1        4.9        3.0        1.4        0.2  Iris-setosa -2.715391  0.169557\n",
       "2        4.7        3.2        1.3        0.2  Iris-setosa -2.889820  0.137346\n",
       "3        4.6        3.1        1.5        0.2  Iris-setosa -2.746437  0.311124\n",
       "4        5.0        3.6        1.4        0.2  Iris-setosa -2.728593 -0.333925"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iris_df = pd.read_csv('iris.txt', names=['sepal_len','sepal_wid','petal_len','petal wid','class'])\n",
    "pca = PCA(n_components = 2)\n",
    "pca.fit(iris_df.ix[:,0:4].values)\n",
    "X = pca.transform(iris_df.ix[:,0:4].values)\n",
    "iris_df['P1'] = X[:,0]\n",
    "iris_df['P2'] = X[:,1]\n",
    "iris_df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAICCAYAAADI/mSeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VGX6+P/3k4QkJCQkoYQqQcACKyWiolISQIqCbf0o\nICi6FlaxoBRdkWpDUNddXXZdGyBE2f19V0GlCUxEUdqAoigqAlJFCCUEkpDJ8/vjTGJIZpKZyZQz\nc+7Xdc0VZs6Zc557MuQ+56lKa40QQgghIldUqAsghBBCiMCSZC+EEEJEOEn2QgghRISTZC+EEEJE\nOEn2QgghRISTZC+EEEJEOEn2QgghRISTZC8ijlIqTSl1l1Lq/ymlflRKnVJKHVNKrVFK3amUUqEu\no/AvpVQvpVSpUmqSD+/d5Xxv2cOhlDqqlPpcKXWfUirazftaKKWeU0ptVErlKaWKlVK/KqVWKKUe\nVEolV3POWyucr6+3ZRbCWzGhLoAQAfB/wGxgP7Aa+AVIB24EXgcGADeHrHTCbLTz8TJwDIgGWgN/\nBC4HegM3VXyDUuou4O9ALPAVsAA4CqQBVwB/BSYCjd2c826gFFDAPcAn/gxIiMok2YtItB0YrLX+\nqOKLSqm/ABuAPyqlbtBa/y8kpROB4I/amr9qrX8pP6BSzwIbgRuUUj201mucr98KvAbkATdorZdW\nKYxSlwL/cFlQpc4DegIrMC4OrlVKNdJa/+aHGIRwSarxRcTRWtsqJ3rn64eAf2IkhixvjqmUOl8p\n9aZSaqdSqtBZXfupUmqUi337KKWWKqWOOPfdrpR61lW1rlLK5qw2jlFKTVJK/aSUOq2U+t5591i2\n3yil1NfOJok9SqkpLo7Vylkt/KazvO87y3DS2YRxlZvYYpVSjzmPX6CUOu6M7f9qOEcrpdS7Sqnf\nnGXeoJS6pprPcKhSarWzivy0UmqbUuoJpVSsi31LlVKrlFINlFKvKaX2Oz/Lb5RSIyvt+xawCuPu\nfEql6vie7spTE631d4DN+fRS57nqAX9znusWV4ne+d71GLUCrtzjfP9bwNsYtQMj3ewrhF/Inb2w\nmjPOnyWevsGZwBZi/FFeilFlmwJ0AsZhXECU7Xsvxh3dSeA/wCGMC4sJwCCl1JVa6xMVDl+2OMW7\nGAnlY2cZbwL+pZQ64zzPCOBDjOrea4FJSqkCrfVMF0U+F/gC+NpZtqbALcASpdRQrfV/KpS3DrAc\n407zO+AVIMF5/veUUp201hNdnCMDWA/sAOZi3KHeAryvlOqrtc6t9Bm+iZHQ9gD/xagu7wZMB3or\npa7SWpdWOkcK8DlQ5Pws4zCaaN5USjm01vOc+/3P+TmOxEjOtgrH2OWi7N4oqzEo+z39H5AKfKG1\nXlndG7XWZyq/5vy8bwOOY5Q7AXgBuAtw9bsUwj+01vKQhyUeGG2xWwEHcJWH72mA8Ye5EOjuYnuz\nCv8+x7nfMaBdpf1exWij/Wel11c7X18HJFV4vTVGksvDSKhNKmyrD/wG/ApEVXi9lfNYDuC5SufJ\nBIqBI0C9Cq8/7nzP4krHagjsdB6rm5tzTKx0jn7ObR9Wen2k8/X/ALGVtk1yHuuBSq+XneNfgKrw\n+oUYF0PfVNq/l/M9k3z4XpTFeU6l1zsABc5tVzpfe935fJqP38EhznL+o8Jr/3EeMzvU/0fkEbmP\nkBdAHvII1gOY5fxDu8iL9zzqfM+LHuz7hPOP9nQX21KcFw0FQJ0Kr692vifLxXtWOrfd7mLbm85t\nLSu8VpaI84BEF+95y/meERVe+xGjlqOdi/3vdB7vdRfn+LliEq6wfRdwqNJrmzEuXJJd7B+FceHy\nZaXXS4F8KlyYVNhmc8aRUOE1fyT7l4DJwDTgnQqJ/j8V9v3I+do9Pn4Hy36nl1V47Rpn2XNC/X9E\nHpH7kGp8YQlKqQeBR4BtGNWonroMowrXZdtsJV2cP1dX3qC1PqaU2gz0AC7AqGGoaJOL4+13/rS7\n2LbP+bMFRtV4RXatdYGL99iA253lnOdsf24D7NVa/+hi/1XOn11cbNuitXa1PvYejOp5AJRSdYGO\nGAl9jItRjwrjQuBCF8f6UWt90s05wKhOP+Viu68edP7UGM0wXwHzMGoXak0p1QajSed7rfW6CpuW\nAgeB65VSaVrrPH+cT4iKJNmLiKeUGo0xFOoboK/W+pgXb09x/txX7V6G+s6fB9xsL3s9pfIGrXW+\ni/3L+hUcr2ZbHRfbfnVz/oPOn/Ur/fS6vBhNFa6UcHbH31SMhN4Io8reHVcXDtWdA4xmGX/K0FpX\nvnCqrOwzae7D8e9x/pxT8UWttUMpNR/jYnQk8KIPxxaiWtIbX0Q0pdTDGL2nvwZ6a6NHvjfKEo4n\nf9zLknITN9ubVtovUNLdvF5WruOVfgayvGXv3ay1jq7mYYYbD0+G733m3K+PVwdWKgajVgXguUqT\n+JRiNBcpjPH3QvidJHsRsZRSEzDukuwYnZ8O+3CYLzH+CA/0YN/NuBnWp5SqD3TG6MD3nQ/l8Eam\nUirRxevZGHfQmwGcVeQ7gObOKubKejt/umpi8IizOeFboINSylUNgb84nD/9fbdf2X8x+kRcrpTq\nXd2OlYYUXo8xwc52jE5+rh4/A+cppXoEoNzC4iTZi4iklHoSeBZjEp2+WuujPh5qDnAC+LOrP8JK\nqYp3/O9g9BR/wEXyfApIBuZpF0Oy/Kw+RkezckqprsAwjJqKipMJvYnxd2CmUiqqwv4NgSf5fTx4\nbbyIMWzuLedFz1mUUilKKVf9ArxxxPnznFoep1rOC6QHMS7qFiql+rnaTyl1OcaFYpmysfVPaq3v\ncfUAnuH3GfWE8CszVJ0J4VdKqduBqRhtu58DD7noGLZLaz2n8ouVaa2PKKWGYQyPWq2UWoLRJJCM\n0fGsBUYnN7TWu53NBq8AdqXUQoyOab0wJljZBjxW+whr9CnwJ6XUZRjxN8OYHlgB91bq9DYLo9bi\nOuArpdTHGGO//w+jnX2G1nptbQqjtX5LKZUJ3AfsUEotw5jCOA1jiGFPjIuO+2pxmu0Y/SqGKKVK\ngN0YyXWuB+3wXtFaL1BKxWP8npcqpbYAazGmy22A8bvuhDHHAkqp1hjV/oeBD6o59HsYfUv+qJR6\nwMu+JUJUS5K9iEQZGH/oo4GH3OyTS6WOUu5orT923hlPwPijfRXGH/bvMe7GKu47Wyn1IzAWYy7+\nBIze4zOAZ/XZE+qUv62603tSxkp2AqOA54B7Me6qN2KMDT9rDnat9RllLMTyCMad/2iMi6QtwINa\n64VuyuRVmbXWDzgvlEZhfIYpGNXhv2B8NvNrcw6tdalS6nqMmG8CkjAubtZQdbSCR2Wudmet33Re\ntIzG+D4MAxIxak6+wfjeldWI/Mn5c67W2u1kTlrrAqVUDsYEO7djzNUvhF8o16NnQksptQujY08p\ncEZrfWloSySE+SmlWmEk+re11neGujxCCPMw6519KcYkI762swohhBDCyawd9BTmLZsQQggRVsya\nUDWwwrmKlow7FcJzNbV1CyEsyKxt9k211geUUo0w1nwerbX+rNI+5iu4EEIIEUBaa08mf6rClHf2\nWusDzp+/YYwJdtlBL9QLC4TyMXny5JCXQWKX+CV+iV/iD96jNkyX7JVSCc4FOnDOAtYPYyiLqGDX\nrl2hLkLIWDl2kPgl/l2hLkJIWT1+X5mxN3468D9nNX0MMF9rvTzEZRJCCCHClumSvdZ6J8Yc4qIa\nI0eODHURQsbKsYPEL/GPDHURQsrq8fvKlB30PKGU0uFadiGEEMJbSil0JHXQEzWz2WyhLkLIWDl2\nkPitGH9GRgZKKXlY5JGRkeH375DpqvGFEEKcbffu3bXujS3Ch1I+3bxXf8xw/QJJNb4Qwiqc1beh\nLoYIEne/b6nGF0IIIYRbkuzDlBXbLctYOXaQ+K0evxC+kGQvhBBCRDhpsxdCCJOrbZv9mTNniImJ\nCUjHr+osWLCAuXPnsnTp0qCeN9xJm70QQgiPzZv3Dq07dCSubl0SU1K55/4HOHr0qF/P0bp1a1at\nWuVy27Bhw4Ka6KdOncptt90WtPOFE0n2YcrK7ZZWjh0kfqvHXyY/P59xj/2F1KbNqRMfz0WXXcH7\n779fvv25WS8wasrT7Hrgr+gvz3B6wVbm/HqGy3plc+rUqfL99u/fz+zZs3n55Zf59ttv/VY+h8Ph\nt2OJ2pNkL4QQYaawsJAr+lzF37/ezbG/f0LJyjy+uXkCtz7wCK/8YzYnTpxg2tNPc+rlpXBJb1AK\nmrSk+LHZ7E9pzrx57wAwcco02rT/A2OXrmPC59u5pPdVXHfLUIqLi70u05w5c+jevTuPPPIIDRs2\nZOrUqcyZM4cePXqU7zNmzBjS09OpX78+nTp1Ytu2bS6PdeTIEQYPHkxqaioNGjSgV69e5dsOHDjA\nTTfdROPGjWnTpg1///vfAVi2bBnPPPMM7733HklJSXTp0qV8/+uuu44GDRpw3nnn8frrr5cfa8OG\nDVxyySXUr1+fpk2bMnbs2PJtN998M02bNiU1NZWsrCy3ZQ0boV6yrxZL/WkhhLCCyn/v3nzzTZ3Y\nra9mQ6lmo/798f9+0AmpaXrhwoU6+fK+Z28rezy3UPccOEi/9957OrHNhZplB3/ftrZQ1+05SD80\ndrzHZcvIyNArV67Ub7/9to6JidGvvvqqdjgcurCwUL/99tu6R48eWmutly1bprt27apPnDihtdb6\n+++/1wcPHnR5zMcff1z/+c9/1g6HQ5eUlOjPPvtMa611aWmpvvjii/VTTz2lS0pK9M6dO3WbNm30\n8uXLtdZaT5kyRY8YMeKsY/Xo0UOPHj1aFxcX6y1btuhGjRrp1atXa621vvzyy/U777yjtda6oKBA\nr1u3rvx9b731li4oKNDFxcV6zJgxunPnzh5/JrXlLr85X/cpZ8qdvRBChJm5/9/7FFx7l3HHXtE5\n7Yg5rxNbt26t8RhPvfQ3Cv78DDRI//3F2DhOj/07/37jdQoLC70uV/PmzbnvvvuIiooiLi7urG11\n6tQhPz+fbdu2obXm/PPPJz093eVx6tSpw4EDB9i5cyfR0dFceeWVgHEnfvjwYZ544gmio6PJyMjg\nrrvu4t1333V5nL179/LFF18wY8YM6tSpQ6dOnbjrrruYO3du+Xl++uknjhw5QkJCApdeemn5e0eO\nHElCQgJ16tRh0qRJfPXVV+Tn53v9mZiFJPswZeV2SyvHDhK/1eMHKHE4IKaO6411Ymnfvj1nvtsE\nB3afvU1rEj96m2HXDebnH76HjldUfX+zDIhP4NChQ16Xq2XLlm63ZWdnM3r0aO6//37S09MZNWoU\nJ0+eZM+ePSQlJZGUlERycjIA48aNo02bNvTr14+2bdsyY8YMwJg2eN++faSlpZGWlkZqairPPvus\n27Lu37+ftLQ0EhISyl9r1aoV+/btA+DNN99k+/btXHDBBVx22WV89NFHAJSWlvLYY4/Rtm1bUlJS\naN26NUopDh8+7PVnYhaS7IUQIszcfPUAEpa9U3XDoX0Ub13HgAEDmPTEEyQ8NADWrwSt4eAeYp8b\nRbNj+xgxYjjpzVrAzy465B39DcfJE6SlpXldrpqG9o0ePZqNGzeybds2tm/fzsyZM2nZsiX5+fnk\n5+dz4sQJAOrVq8esWbPYsWMHixYt4sUXX2T16tW0bNmSc889l7y8PPLy8jh69CjHjx9n8eLFLs/f\nrFkz8vLyKCgoKH/tl19+oXnz5gC0adOGBQsW8NtvvzF+/HhuuukmTp8+zfz581m8eDGrVq3i2LFj\n7Nq1q2ITcliSZB+msrKyQl2EkLFy7CDxWz1+gJEjb6fhnu+o8/I4OHbYSObfrCfhkWsYO+YRUlJS\neGzso/xzyhNkvDIG1a0OdYddxMgmcazLXU1CQgJj7r2bhH9PhsLTvx9Ya2L/NYkb/3gT9erV82uZ\nN27cyPr16ykpKaFu3brEx8cTFeU6BX300Ufs2LEDgKSkJGJiYoiKiuLSSy8lKSmJ559/nsLCQhwO\nB99++y0bN24EID09vTwxA7Ro0YIrrriCxx9/nKKiIr7++mveeOMNRowYAcD8+fPL79br16+PUoqo\nqChOnjxJXFwcqampFBQU8Pjjjwd9jgK/87WxP9QPpIOeEMIiXP29O3jwoL75tpE6NrGejk2urxu3\naq3//uo/dGlpaZV9i4uLq7xeUlKibxx6q05s1VarPz+leeQlXa9TN93+4kv00aNHPS5b69atyzvo\nlXXGK1PxtZUrV+qOHTvqpKQk3ahRIz18+HBdUFDg8pgvvfSSzsjI0PXq1dMtW7bUTz/9dPm2AwcO\n6KFDh+omTZrotLQ0ffnll+uVK1dqrbU+cuSI7t69u05NTdUXX3yx1lrrvXv36kGDBum0tDTdtm1b\n/dprr5Ufa/jw4bpx48Y6KSlJ/+EPf9CLFi3SWmt98uRJfd111+mkpCSdkZGh582bp6OiovSOHTs8\n/lxqw11+oxYd9GQGvTBls9kse4dj5dhB4rdi/NXNoFdUVERBQQEpKSlu75Td0Vqzdu1aFiz8L4XF\nxQzufxWDBg0iJkZWPw+lQMygJ79RIYQIY3FxcVV6vntKKcWVV15Z3ttdRC65sxdCCJOT9eytRebG\nF0IIIYTXJNmHKSuPNbZy7CDxWz1+IXwhyV4IIYSIcNJmL4QQJidt9tYibfZCCCGE8Jok+zBl5XZL\nK8cOEr/V4xfCF5LshRAiwp05cyYkzQALFixgwIABQT+vp/xRvtzc3GoXADILabMXQgiT87XN/p15\n83h+2lS2/byTenXjGTF8ONOefY7U1FS/la1169a88cYb9O7d22/HDCe5ubmMGDGCX375xW/HlDZ7\nIYQQAOTn5/OX8eNo3rAB8bF1uKJLJ95///3y7S/MmMHTYx7kr42iOHN1Z7Ze3oYznywm+4rLOXXq\nVPl++/fvZ/bs2bz88st8+62LVfB85HA4/Has2jDzTWFpaWnQziXJPkxZud3SyrGDxG/1+AEKCwu5\nqkd3dv/nHT65qAl5fS9iQp0CHrlzJLNfeYUTJ07w9FPTWdq5Bb0bJqOUomXdWGZf2JTmp47xzrx5\nAEyb9CR/OK8d6154mu2vzuSqK7ox9MYbKC4u9rpMc+bMoXv37jzyyCM0bNiQqVOnMmfOHHr06FG+\nz5gxY0hPT6d+/fp06tSJbdu2uTxW+/bt+fjjj8ufOxwOGjduzJYtWwD48ssvufLKK0lNTaVLly7k\n5uaW75udnc3EiRPp3r07iYmJ7Ny5k7fffps2bdqQnJxMmzZtyMnJKS9zxfJ9++239OvXjwYNGtC0\naVOee+45AIqLi3n44Ydp3rw5LVq0YMyYMZw5c8Zl2b///nuys7NJTU3loosuKl9+F+COO+7gvvvu\n45prriEpKSmo32WZG18IIcJMTk4OSUcO8k7nFuVLr17XJIX29eLp9pfHaZiezsVpybRKOHvOfKUU\nIxslMnfhe6SkpvLu7Ff57sq2pMfVAeCldo25aeNannxsAjNefMnrcq1bt45hw4Zx6NAhzpw5w7vv\nvltevuXLl/PZZ5/x008/kZSUxPbt20lJSXF5nGHDhrFgwQKuvvpqAJYuXUqjRo3o3Lkz+/btY9Cg\nQcyfP5/+/fuzcuVK/vjHP7J9+3YaNGgAwDvvvMPSpUs577zzOHnyJA899BCbNm2ibdu2/Prrr+Tl\n5Z31mQCcPHmSq666ivHjx/Phhx9y5syZ8ouRp556ivXr1/P1118DcO211/LUU08xderUs8pdUlLC\n4MGDueuuu1ixYgVr1qzhuuuuY9OmTbRr1678d7dkyRK6devm00WVr+TOPkxZbdWviqwcO0j8Vo8f\n4P0F87mrcWKVNdbb1YunU2o9tm7dWuMx/vbcszyTkVae6AHioqP4+3mNef311yksLPS6XM2bN+e+\n++4jKiqqyuI8derUIT8/n23btqG15vzzzyc9Pd3lcYYOHcqiRYvKy5CTk8PQoUMBYw36a665hv79\n+wPQp08funbtelZNwMiRI7nggguIiooiJiaG6Ohotm7dSmFhIenp6Vx44YVVzvnhhx/StGlTHn74\nYWJjY0lMTOSSSy4BjI58kydPpkGDBjRo0IDJkyczz1k7UtEXX3xBQUEBEyZMICYmhuzsbAYNGlRe\nkwBw3XXX0a1bNwBiY2M9/mxrS5K9EEKEGYejhDpRrvtpxaoo2rdvz6a8E+w+VXTWNq01b/9WwOCb\nb+H7HTu4Ii2xyvszEuJIiI7i0KFDXperul7p2dnZjB49mvvvv5/09HRGjRrFyZMn2bNnD0lJSSQl\nJZGcnAxAmzZtaN++PYsXL+b06dMsWrSIW2+9FYDdu3ezcOFC0tLSSEtLIzU1lc8//5yDBw+6LEdC\nQgLvvfces2fPpmnTpgwePJjt27dXKd+ePXto06aNy7Lv37+fc845p/x5q1at2L9/f5X9Dhw4UOUz\naNWqFfv27fPoMwokSfZhysrtllaOHSR+q8cPMOCGP/LOb6eqvL7vdDHrDh9jwIABPDHxSQZs2cvK\n306gtWbP6WJGfXeAfQkpDB8xghZN0vk2v+rd+29FZzhRdIa0tDSvy1W5pqGy0aNHs3HjRrZt28b2\n7duZOXMmLVu2JD8/n/z8fE6cOFG+75AhQ1iwYAEffPABHTp0oHXr1oCRLG+77Tby8vLIy8vj6NGj\n5OfnM27cOLfluOqqq1i+fDkHDx7k/PPP55577qlStpYtW7Jjxw6X5W7evDm7d+8uf757926aNWtW\nZb9mzZqxZ8+es1775ZdfaN68uduyBYskeyGECDO3jxzJd1HxjNt+kMNFJWitWX+0gGu+2scjjz5K\nSkoKj06YwBMv/Y0xh0up8/EWLvpiB3FXXcvqtV+QkJDA3Q88xOTdRznt+L1HuNaaST8f5qYbb6Be\nvXp+LfPGjRtZv349JSUl1K1bl/j4eKKi3KegIUOGsHz5cmbPns2wYcPKXx8+fDiLFy9m+fLllJaW\nUlhYSG5urss7bYBDhw6xaNEiTp06RZ06dahXr57L8w4aNIiDBw/yt7/9jeLiYk6ePMn69evLy/LU\nU09x+PBhDh8+zPTp0xkxYkSVY1x22WUkJCTw/PPPU1JSgs1m48MPPyxvgggprXVYPoyiCyFE5HP1\n9+7gwYN65NAhul58vK5fN163btpE/+OVv+vS0tIq+xYXF1d5vaSkRN960026bYMU/dQFLfRLHVro\nbk0b6Usu+oM+evSox2Vr3bq1XrlypX777bd1jx49ztpW8bWVK1fqjh076qSkJN2oUSM9fPhwXVBQ\nUO2x+/Tpo2NjY/Wvv/561uvr16/XvXr10mlpabpx48Z60KBBes+ePVprrbOzs/Ubb7xRvu+BAwd0\nr169dEpKik5NTdXZ2dn6u+++q1I+rbX+9ttvdZ8+fXRqaqpu2rSpnjFjhtZa68LCQv3QQw/ppk2b\n6mbNmumHH35YFxUVaa21ttlsumXLluXH2LZtm+7Vq5euX7++7tChg/7ggw/Kt91xxx36ySefrPEz\ndZffnK/7lDNlUh0hhDC56ibVKSoqoqCggJSUlGrvlF3RWrN27Vr+m5NDcVEhV10ziEGDBhETIwO1\nQikQk+pIsg9TNpvNsr2SrRw7SPxWjF9WvbMWmUFPCCGEEF6TO3shhDA5ubO3FrmzF0IIIYTXTJvs\nlVJRSim7UmpRqMtiRlYea2zl2EHit3r8QvjCtMkeeAhwvUqCEEIIITxmyjZ7pVQL4C3gaeARrfW1\nLvaRNnshhCVIm721WKnN/iVgHCDfbiGEEKKWTDdzglLqGuBXrfUWpVQW4PYqZuTIkWRkZACQkpJC\n586dy8fflrXrRerzv/71r5aKt+Lzim22ZiiPxC/xB/p5bZ05c4aYmJigz8u+YMEC5s6dy9KlS/12\nzD179tChQweOHz/uUzzevL+256qNsu+AzWZj165dtT6e6arxlVLPAMOBEqAukAT8P631bZX2s3Q1\nvs2CE4uUsXLsIPFbMX5fq/HnzZvHpGensPuHXcQnxDF8+HBmPD2D1NRUv5WtdevWvPHGG/Tu3dtv\nx7Q6y82gp5TqBTwqbfZCCCtz9cc/Pz+faU9P4/W33+Dk0Xwu6HQh0/8yjeuvvx6A52Y+x/RXn+HU\n3bHQMRYOlxL73yJa/dKILeu2kJCQABjLt37wwQcUFxfTt29fOnTo4FXZ3CV7h8NBdHR0LaL2jdY6\nZCvL+YuV2uyFEEK4UVhYyBW9r+Tv6//NsckxlMxvyDf99nLr/SN45dVXOHHiBNOens6pSXWhUxwo\nBY2iKR5Vl/2Jh5k3bx4AE6c8SZsL2zL2fxOZ8MlULsm6jOtuvp7i4mKvyzRnzhy6d+/OI488QsOG\nDZk6dSpz5syhR48e5fuMGTOG9PR06tevT6dOndi2zfWAq/bt2/Pxxx+XP3c4HDRu3JgtW7awe/du\noqKiKC01VuvLzs5m4sSJdO/encTERHbu3MmuXbvo1asX9evXp1+/fowePbp8lTpX7580aRLdu3cn\nOTmZAQMGkJeX53Lfo0ePcuedd9K8eXMaNGjAjTfeCMCxY8cYPHgwjRs3pkGDBgwePNjtKnyhYupk\nr7XOdXVXL6w91tjKsYPEb/X4AXJyctip91I0JgFa1oH4KOhWl1OTE5gw8TGWLVtGnXYJ0LhStyyl\nKMiCBe/nsHDhQv46528UvpLMqQfjKbq3LqdfS2bFnlzG/2W8T+Vat24dbdu25dChQzzxxBPOUxo3\nosuXL+ezzz7jp59+4vjx4yxcuJAGDRq4PM6wYcNYsGBB+fOlS5fSqFEjOnfufNYxy7zzzju8/vrr\n5Ofnc8455zBs2DC6devGkSNHmDx5MvPmzTvrPZXfn5OTw5w5c/jtt98oKipi1qxZLvcdPnw4p0+f\n5rvvvuPQoUOMGTMGgNLSUu6880727NnDL7/8QkJCAqNHj/b68wskUyd7IYS5OBwOHA5HqItheXP/\nO4+CPhh37BU1iyHm3Hi2bt1a4zGeevEZCoZFQ2qFqvY6itN3xfLvN16nsLDQ63I1b96c++67j6io\nKOLi4s7ky3muAAAgAElEQVTaVqdOHfLz89m2bRtaa84//3zS09NdHmfo0KEsWrSovAw5OTnVrgk/\ncuRILrjgAqKiojhw4AAbN25k6tSpxMTEcOWVV3LttdXfM95xxx20adOGuLg4br75ZrZs2VJlnwMH\nDrBs2TL+9a9/kZycTHR0dHmtRVpaGjfccANxcXEkJiby+OOPk5ubW+05g02SfZiyWgeliqwcO4Qm\nfrvdTvbA3sTGxxIbH0v2wN5s3rw56OUA+f0DlDgc4K45PEbRvn17zvx4Cg6VnL1NaxJtMOz6ofz8\nww64MLbq+9NjIE5x6NAhr8vVsmVLt9uys7MZPXo0999/P+np6YwaNYqTJ0+yZ88ekpKSSEpKIjk5\nGYA2bdrQvn17Fi9ezOnTp1m0aBHDhg3z6Lz79+8nLS2N+Ph4j8oF0KRJk/J/JyQkcPLkySr77N27\nl7S0tPIyVnT69GnuvfdeMjIySElJoVevXhw7dsxUcyNIshdCVMtut9Ozby9s59opzWlMaU5jbOfa\n6dGnJ3a7PdTFs6Sbr7mJhDUuNhxxUPx9AQMGDGDSE0+SMO00fFUEWsNvDmJnn6ZZQUNGjBhBerN0\n+KWk6jGOO3AUlJCWluZ1uWrqGDd69Gg2btzItm3b2L59OzNnzqRly5bk5+eTn5/PiRMnyvcdMmQI\nCxYs4IMPPqBDhw6ce+65Hp23adOm5OXlnVUzsWfPHq9jqaxly5bk5eWdVcYyL7zwAj/++CMbNmzg\n2LFjfPrppwCS7EXtWbnd0sqxQ/Djf/SJsRQMjYaBiUbbcHwUDEykYGg0YyeOC2pZQH7/YFRbN/w1\nmTpvn4LjDiOZ/1BMwlOnGfvIWFJSUnhs3GP8c/qrZMyvj7rhV+o+dIKRrW9m3afrSEhIYMyoh0l4\nzwFFFRKS1sTmFHHjH2+gXr16fi3zxo0bWb9+PSUlJdStW5f4+HiiotynoCFDhrB8+XJmz55d5a6+\nuiR6zjnn0LVrV6ZMmcKZM2f44osvWLx4scfvr6xs3yZNmjBw4EDuu+8+jh07xpkzZ1izxrjiys/P\np27duiQnJ5OXl8eUKVM8Pn6wSLIXQrjlcDj49JNcyK5bdWN2XXJX2KQNPwSSkpJYv2YdNyT0J/be\nY8Temkfjl6OZ8fDTTJs0tXy/ESNGsHPbDopOF1JwvIB/vfqv8jH2f/7znxnQvg+JYwpQCwvgg5PU\ne7yQtgea84+//sPjsng6zO3EiRPcfffdpKWl0bp1axo2bMi4ce4vFps0acLll1/Ol19+yS233OL2\nnK7OP3/+fNauXUvDhg2ZNGkSQ4YMOasPQU3vd3euefPmERMTwwUXXECTJk14+eWXAXj44Yc5deoU\nDRs25IorruDqq6+u9pihYOpx9tWRcfZCBJ7D4SA2PpbSnMbGHX1FhaVEDT1EcWFxSMZTW0l1k+oU\nFRVRUFBASkpKtXfKrmitWbt2LQsW5lBYdJrBAwYzaNAgYmJMN7lqrQwZMoQLL7yQyZMnh7ooHrHc\npDrVkWQvRHBkD+yN7Vy7UY1f0ZICsndezKqPV4amYBYiC+F4Z+PGjeU1CMuWLePGG2/kiy++oFOn\nTqEumkdkUh1RzsrtllaOHYIf/4vPvEBijgOWFEBhqfFYUkBijoMXnp5V8wH8zOq/f1GzgwcPkpWV\nRVJSEg8//DD//Oc/wybRB4okeyFEtbp06cKnn+SSvfNiooYeImroIbJ3XsyalZ/SpUuXUBdPiCoG\nDRrEL7/8wsmTJ/n++++57bbban5ThJNqfCGEx8o640kbfXBJNb61BKIaP7J6YQghAkqSvBDhSarx\nw5SV2y2tHDtI/FaPXwhfyJ29EEKYXKtWrcJ+2VbhuVatWvn9mNJmL4QQQoQBGXonhBBCCLck2Yep\ncGi3DNRyqOEQeyC5it9KS8/K798W6iKElNXj95Uke+F3ZloONdLJZy2E8IS02Qu/KlsOtWBo9O+L\np6w+TWKOsaBKZmZmaAsYQeSzFsJaZG58YRoyj3rwyGcthLVIBz0LMmO7VbCWQzVj7MFks9ksvfSs\n/P5toS5CSFk9fl9JshfChKzU4U4IEXiS7MNUVlZWqItQRXR0ND379oLVp6tuXH2aXldl+WW6VTPG\n7i+edLjLysoK2mdtRpH8+/eExJ8V6iKEJWmzF361efNmevTp6bLTmKySVj1vO9zJZy2EtUibvQWZ\ntd0qGMuhmjX22nr0ibFG4h6YCPFRxmNgIgVDoxk7cVz5fmXxW3Xp2Uj9/XtK4reFughhSebGF36X\nmZnJqo9XynKoXijvcJfTuOrG7LrkvmZ0uKv8WcpnLYTwhFTjC2ECDoeD2PhYSnMaG3f0FRWWEjX0\nEMWFxZLMhbAwqcYXIsxZucOdECLwJNmHKSu3W0Vq7C8+8wKJOQ5YUgCFpcZjSQGJOQ5eeHpW+X6R\nGr+nJH5bqIsQUlaP31eS7C1IxnCbk1U73AkhAk/a7C3Ebrfz6BNjjY5gQM++vXjxmRckkZiQdLgT\nQlQmc+OLGsmiKUIIEd6kg54Fedtu5ekY7nBg9TY7id8W6iKElMRvC3URwpIkewuw8qIpQgghpBrf\nEmQMtxBChD+pxhfVkjHcQghhbZLsw5S37VaejuEOB1Zvs5P4baEuQkhJ/LZQFyEsSbK3CBnDHTwy\nj4EQwmykzd6CZAx3YMg8BkKIQJJx9kKEmMxjIIQINOmgZ0FWbrcyY+zBnMfAjPEHk8RvC3URQsrq\n8fvKdMleKRWnlFqnlNqslNqqlJoc6jIJUR2Zx0AIYXamrMZXSiVorU8ppaKBz4EHtdbrK+0j1fjC\nFGQeAyFEMERcNb7W+pTzn3FADCBZXZiWzGMghDA7UyZ7pVSUUmozcBBYobXeEOoymY2V263MGHsw\n5zHwV/zhOkTQjL//YJL4baEuQlgyZbLXWpdqrbsALYDLlFLtQ10mYV2eJMVwmsfAbreTPbA3sfGx\nxMbHkj2wN5s3bw51sYQQAWTKNvuKlFJPAgVa6xcrva5vv/12MjIyAEhJSaFz585kZWUBv1/9yXN5\n7uvzH374gZz/vcunn+SiSzUdL+7EW/96ky5dulT7fofDgc1mIzo62lTxACQnJxtDBHtq6BILHeJg\n9Wni5xbx8gsvc88995iqvPJcnlv5edm/d+3aBcCcOXMiZ5y9UqohcEZrfVwpVRdYBjyntf640n7S\nQU8ETKSOm88e2BvbuXZjiGBFSwrI3nkxqz5eGZqCCSFqFGkd9JoCq5VSW4B1wLLKiV5Yu90qGLEH\nc9y8t3yN3+FwkLvCBj3jq24MoyGCVv7ug8Rv9fh9Zbpkr7XeqrXO1Fp31lp31Fo/HeoyCWuJxHHz\ndrudPtf0RZdquPVXmHYEdpwJdbGEEEFiump8T0k1vgiUSBs3765JgnfyYWoatI2VanwhwkBtqvFj\n/F0YIcJd2bh522oXbdthOG7+rCaJMmX/np8Pl8YbQwRXhtdSx0IIz5muGl94xsrtVsGIPZjj5r3l\nTfw1NUmwqYisnzNNOUTQHSt/90Hit3r8vpJkL4QL4TRuvjaioqP45MMVERWTEKIqabMXogZlnfHC\nqeq+IhluJ0RkkPXshRBubd68mR59erqcMyDSaiqEiGSRNs5eeMDK7VZWjh28jz/SmiTk928LdRFC\nyurx+0p64wthAZmZmaz6eGXYN0kIIXwj1fhCCCFEGJBqfCGEEEK4Jck+TFm53SrSY69pSV1/xh+O\na9pH+u+/JhK/LdRFCEuS7IUwiWCuMy9r2gthLdJmL4QJBHNJXW/PJZ36hDAHabMXIswFc0ldT88l\nd/9CRA5J9mHKyu1WkRa7t0vq1iZ+T89VdvdvO9dOaU5jSnMaYzvXTo8+PbHb7T6f3x9sNltY9jXw\nl0j7/nvL6vH7SpK9EKKKYNY0eMNut/Pw+DFS2yCEl6TNXggTCOb89TWda8Xi5cTGx1Ka09hI8hUV\nlhI19BDFhcVBb8MPZr8GIcxI5sYXIswFc/76ms7VsWNHUyZ7WdBHWJ100LMgK7dbRWLs3sxfX9v4\nazpXdHQ0Pfv2gtWnq7559Wl6XZUV9ER/Vl+DrUVnb3TRryGSReL33xtWj99XMje+ECYRjPnry45d\n3bkcDgczpz9PVr9sCiiocvf/wspZfi+XECKw5M4+TGVlZYW6CCET6bFHR0dXm+h9id/dMLqK56q4\nz2WXX8YF7S+g65bzTLFS3lm1DRfFnb0xRLUNoRLp3/+aWD1+X0mbvRBhrqaaAE86tlW3z+plq8jM\nzAx5Mg1mvwYhzEja7C3Iyu1WVo4dfo/f00lvPBlGV90+EyY/FvJED7/3Nei88VxT1DaEinz/baEu\nQliSZC9EGPJ00htPJtEpLi72alKfUMrMzOSlGS9SXFhMcWExqz5eSceOHU1TPiHMSqrxhQhDng5D\nczgcNQ6jO11wmrqJdU031K4mdrudR58Ya1yoAD379uLFZ14I2F2+rBEgQk2q8YWwEG+m1/VkGF1s\nbKxfh9oFYyrbYE7nK2sEiEggyT5MWbndysqxg/fxv/jMCyTmOGBJARSWGo8lBcYwuqdnebxPTYKV\nFG02W9Cm8zXjGgHy/beFughhSZK9EGHG20lvPJmwp2PHjqxetsqjSX1cCWZS9HbhoNow6xoBQnhL\n2uyFCEO+DkOr3O7sqt175vTny2fS81Qwp7L1pB+CP/oYBOs8QnhK2uyFsBhvptetqPIkOq7uxrP6\nZfPVV195XJZg3mmD9zUbQghJ9mHLyu1WVo4dfo+/bMrbisPQvOmJHq5V1DabzS99DGpi1osK+f7b\nQl2EsCTJXogwV9P0uq748248FEnR15oNbwXjokKIYJA2eyEsyJ/t0Q6Hg82bNxsL54RgKttALxxk\nt9uZMPkxclfYAOh1VRYvPD3LUrP2CXOQNnshhFf8cTdeZeGcCy6g61bXC+cEcuy9LzUbNakYW7cr\nuqG1Zt0X63xqLhHCDCTZhykrt1tZOXbwX/y1qaJ21blvU5cf+e6b7/hy7ZflSVFr7fex94H+/fur\n42KgyPffFuoihCVJ9kJYTNlddm3avT1ZOMeME9J44p4H7w3LjotCVEfa7IWwiOrmkvem3dvT9v6+\ng64K2th7f9mwYQOXXnYpLGwiY+uF6UibvRCiWjXdZfu73TvYY+/9ZdyT48GnP6VCmJsk+zBl5XYr\nK8cOvsXvzzH1oR5/Hqjfv8NhjBygS5zb2Hr27RXyu3r5/ttCXYSwJMleiAgXiLvsmjr3hfqCoFZu\nqQfv5FeJjddPMHP686EunRA+kTZ7ISJcoOZ4t9vtjJ04zu34c1/n7w+l8jn+29WBnHzYVGRsyIih\na3oXNny+PrQFFJZWmzZ70yV7pVQLYC6QDpQC/9Za/83FfpLshfBQIBeqqa5zX00XBGZT5QLFocF2\nmsT3Sk17gSKsI9I66JUAj2itOwCXA/crpS4IcZlMx8rtVlaOHXyLP5DTvlbXua+28/e7Esg2+44d\nO549HHH4b2Tv7mqqRC/ff1uoixCWTJfstdYHtdZbnP8+CXwHNA9tqYQIb8GaS97dTHm+zt8fjB77\nFWfLi42P5dEnxvLC07P8eoEiRKiZrhq/IqVUBmAD/uBM/BW3STW+EDVwVcUeiLnkqxvDH8pjeXKu\nnn17uexX8OknuWRmZvr9nEL4KqLa7MsopephJPrpWusPXGzXt99+OxkZGQCkpKTQuXNnsrKygN+r\neuS5PLfi89dee41/vD6brZu/BuCiLh25/+77uPvuu/1+PrvdzpW9rqSwTzQMTwKAd/KJX+ng89zP\nyczMPGv/4uJicnNzqVOnTpXjJScnG8m3eyl0joOOxjC4+LlFvPzCy9xzzz1+Lf/UGdOMvgwtYoxy\nXxRn/Hz9OJ0PnMfmdXa/f17yXJ57+rzs37t27QJgzpw5kZXslVIxwIfAEq31y272sfSdvc1mK/9i\nWI2VY4ea4w/23aqnnf/mz5/P/eMf4PiBowDUb5rK7FmvMnTo0PK3dL3iEjYd2gK7SowXMuPg1mT4\nobj8WP76/QdqlEKgyfffuvFHWgc9gDeBbe4SvRDCPX9OoFMTT8fwz58/n+F/GsHx60uNqWgXNuH4\n9aUMu+NW5s+fDxhT1W6yb4L+ifBuE+NxSTxMOgKtYkw7654Q4cB0d/ZKqSuBT4GtgHY+/qK1Xlpp\nP0vf2QvhSrDvVj09X4NzGhmJ3sXdf8oHURzdm0fWgGxy22x2uQ/rC4nacsbvd9qBHJIohL9F1J29\n1vpzrXW01rqz1rqL1jqzcqIXQpiDJzPlORwOo+rezd3/sf1HOX36tDFVrZt9sBcFZKraQA5JFMJM\nTJfshWcqduCwGivHDtXHH4ppaoOSMDXlU9X68/df05DEYA3/84Z8/22hLkJYkmQvRIQJ9t1qTQkz\nNjaW+k1T3V6ApDRLpW7dutVepHTtfgldu3b1e9nB9cQ/Wuuzxt5nD+zN5s2bA3J+b5jx4kOECa11\nWD6MogshXNm0aZPOHthbR8VE6aiYKJ09sLe22+0BP29JSYkuKSmp8vqCBQs0cUpzX33Nf5oYj/vq\na+KUXrBggS4pKdEbNmzQian1quyTmFovKGUvs2nTJrfl2LRpU9DKUblMWQOyy3+fWQOyg/qZCHNw\n5j2fcqbpOuh5SjroCVGzQEyg46v58+czesIDHNtvDL1LaZbKuPvGsmLNJ+UT6HS5LBOlFfb1m4DQ\nzKVvtk57MvGPKBORk+rUxOrJ3spjTa0cO4R//MXFxQB88803bpPY6mWryMzMdHmREsj4zTj2vsrF\nx9YiY/Ifi44YCPfvf21EVG98IURki42NJTY2ttr5ACZMfswUtRGh5uk8BkLURO7shRBBZ8Y76DJm\nqsY38+ckgk/u7IUQQRPpPcLNNPY+FEMpRWSSZB+mrDzW1MqxQ+jir7wUbG2Go9UmiQU6/mAtB+yp\nKhcfm05beuIfq///95UkeyFEjcp6hNvOtVOa05jSnMbYzrXTo09P7Ha7T8c00x10Za7G3odqTfvK\nFx9q2rGQXnyI8CRt9kKIGgWqHdtutzN24jhyV9iA0Ay1CydmGkopgk+G3gkhAiYYncQkiflOPjvr\nkA56FmTldisrxw6RGX90dLTHySoS4/dGWfz+7EMRTqz++/eVJHshhFt2u53r+/ejtI6WHuEmEog+\nFCKySTW+EMIlu91O/6xeTM9Io2NSXa7a/BOn7qhXZbY76SjmOX9VuZtpLgARPFKNL4Twu8njxzE9\nI41RGY24okE91nQ9j+wPNOqWg3DLQekR7gV/VrnLrHrCF5Lsw5SV262sHDsEJ36Hw8FSWy63tWhQ\n/lpmSgKrMttxvF9noolixeLlIUn04fb793eVe7jF729Wj99XkuyFCCNmmL0uWimU8qkm0ZKqWwNg\n7MRxXh9PZtUTvpA2eyHCgN1uZ/L4cSy1GUvBDsjqxbSZgR2PPrhvH6459BOjMhqd9fo/d/3Gx+nt\nWLTik4CdO9CCNVwtUMMWN2/eTI8+PV2uGChNK5FL2uyFiGBlHeWuOfQTx/t15Hi/jlxz6Cf6ZQW2\n5/W0mbN4cnce/9z1G6dKSjlVUso/d/3Gk7vzmPr8zICdN5AiZbhadVP6dujQoXwZYSHKSLIPU1Zu\nt7Ja7BU7yiXERLH+WAGjMhoxvVUaUyaMD9h5u3TpwrLVuXyc3o76y7+m/vKv+Ti9Hcttob1z9PX3\nH4jhajU1qwSiyr0s/spT+v7p1jvJHtSHuPg44uLjSGmeRk5OjlfHDgdW+//vL5LshTAxVx3lytzW\nogFLVge253VmZiaLVnxCYVERhUVFLFrxSdhWEXvTdl5TEvemhsCTNQBq0xcjOjqad999l+F/GsHx\n60thYRNY2ITj15cy7I5bmT9/vk/HFZFF2uyFMDGHw0F8XBzH+3UkIebsa/NTJaXUX/41hUVF0iGr\nBp62nX/11Vc8+sRYY2gb0LNvL1585oWzLnDKaghctZd/+kkumZmZVc7vbg0ArbVf+mKkNE8zEr2L\ncfcpH0RxdG+eV8cT5iRz4wsRwczWUS4c52L3JNl/ufZLsvv3rjGJ12ZCm4qfXcVJi8pqbubuPcKT\nu/NYttr1RYMrxcXFxMXHGXf0LmLj5oMUFRYRGxvr0fE8FY7fg3AnHfQsyMrtVlaLvXJHuaW/Hg9J\nRzm73c7gvn2Ij4sjPi6OwX37hKRzmy+/f0/azsdPmlBjNX9tJ7SpuAZA5b4YCTFRHvXFCPX3P9Sd\nHEMdf7iSZC+EyVXuKHfNxp+D3lEuVCMC/Gnm9OdJyClx2Xb+/LQZfp2VrqY2eH/2xYiNjaV+01S3\nFzIpzVL9dldvtxsdGm0Zm2RO/jAj1fhC+FkgqzdDVXVqtqYEb9jt9vJ2eK019dNTOPHrcVC/t513\n7NjR4/Hw1VXjd/36PJrEJdXYBu/vvhg5OTkMu+NW+FPyWU0QvHGCBW/NZ+jQoZ5/YG7Y7XZ6D+7L\n8YNHQQGZcXBrMrSpI3PyB4lU4wthAsGo5vZmKVh/CfWIgNqoPNxOv5vOsetLiU+qy5drv2TVxyvp\n0qUL0dHR9OjTE1aeqnqQSkPk3PWuj59fzE/2bzyq/YiOjmZAVi/m7j1S5XRz9x5hYLZ3Q/KGDh3K\nO2/MI+WDKLj5INx8kJQPovya6Hv07Wl0AnyvCbzbBC6Jh0lH4KdimZM/DEiyD1NWbrcyY+zBrOY2\nY/zB5E387obbnRoWw4TJjwG/t0Gv+eRT+NcJePgQfFfkcogcuJ/Q5pLz/8CzbRp53Abv66RF7uK/\n9dZbObo3j6LCIooKizi6N88viR6Mz/HU0JgqnyPDkyAn3y/n8JTVv/++kmQvhB/42tkq0Pwxl76/\n70KDxZPOdBs2bPj9zv/dxkaP9v6JMDEPNeRXtyv7VZ7QZsXi5Xyxye5V7UegJi2KjY31a8/7mj5H\nNhXBylMyJ7/JSZu9ELVkxrHw/p5Lf/PmzfTL6sn0VlWHiYV6Rj13PBlu16NPT3LbbHbZ/p71cyar\nl6zy+Fy1+Q6YuZ9HTZ8jtxwkITmRz1atMeX3IJJIm70QolwgmhTMOnVudWoabtezby/WrPwUesaD\no9KNQ3ZdPv0k1+NakdrWfgSiL4a/hsjV9DmmNE2VRB8GJNmHKSu3W5kt9mBXc9cUf6CaFMwyda43\nv//qpqr90613UlpHw62/wh8PwLQjsOOMz+UK1sJBnsTv73UA3H2OCTklrFq8MqjfA7P9/w8XkuyF\n8AOzrBAXjJ7zoRgR4Ct3nelee+VfjHrwz3BHstGzvHLvch8WqTFT7Yc36wB4wt3n+NlKuaMPF9Jm\nL4Sf2O12pkwYz5LVNgAGZmcx9fmZQf1jaMb+A6FSua264vPqxsqzrIDEvPharQsfyqlkPV0HoDZt\n+CDT5IaCtNkLYQJmqOb2d5OCP3rzB/s87uY7KKuRqLF3+c4SbMtX1+p3F061H96K5NgimST7MGXl\ndiuzxx7oP4Y1xe+PJoVgzYPvy3mqi98fnROjoqJMXTVd0+/fk3UAwjlZm/3/v1lJshciwtS27ThY\nEwQF4jyedE6M9GQI1XdMrDhBkLAOabMXIoL50r7q6Tz4tW279fd8+970V9i8eTM9+vR0uZxtbdrq\nzcRutzN24jhyV9iA39cBiITYrErWsxdC+IUnCXPtF18w7fHHajVhTyA6Enp7TKskQ+lQFzmkg54F\nWbndysqxQ2jj11pzdZ/eIV3q1l383nZOrDzlbdmiOGbn7e8/0jrUWf3/v69MmeyVUm8opX5VSn0d\n6rIIYSU1JcwWaal+mbAnUBMR1dQ50VWv/0hLhkK4YspqfKVUd+AkMFdr3dHNPlKNL0QAuJsHf+Ku\nIxwvOM3x/p38UvUeqPn2Xc13MOSOO8l58w2/rRUgRChEXDW+1voz4GioyyGE2flzfHrZsdz15l+6\nygbKp78zLgVqxrnK8x1MmfE8D426N6RND4EUrLkQRHgzZbIXNbNyu5WVYwcjfn+Og3d1LKVUlQmC\nunbt6veqd18mIvL0919WPe/pWgHhkjTL4g/WXAhmY/X//76KCXUBamPkyJFkZGQAkJKSQufOncnK\nygJ+/0JE6vMtW7aYqjzyPHjPf/jhB8Y9/BB/apLMe/2MVq6J339DVvcrWL3mczIzMz0+XnJyMv2z\nenFrg7p8kJlBVsNk5u79iazuVzDzpZe55557ztp/2sxZ9MvqyQ8nC+nXKJmeDZKYu/cIE3YcYtYD\nEynjbXxr1qzx6fPo0aNHte/v0aMHS2253H9xa2yH88lqmGRsP5xPRt1Ylqw21rR/cNS9rN/yFVFK\nMSCrF9feMoR27dqZ4vft6vlrr73GuIcfYkbbxrzXryOfHsln+c/f0C+rJ8tW53LixAlTlVee+/a8\n7N+7du2itkzZZg+glGoFLJY2eyHO5s/x6T0vvZSC7d/wdb4xwcyAxvWZdn4z1h0tcHssM6wBYLfb\nmTx+XI1t8DUNx0te9hWpiQlMz6jab2DZ6lwyMzODE5CX/D1HgQgPETnOXimVgZHsL3KzXZK9sBx/\njk/fsGEDvS7vxovtW5yd6Lbv5/2ubcha91O1xwrV+O2ymfc8TdDVJcbnDhTwWNPEsEqastiRdUVc\nBz2l1AJgLXCeUuoXpdQdoS6T2VSs5rEaK8cOUOqni9ypj03gxfYtqrZln9+MZ346UOP7QzVk7f67\n7/Jq+J+74XgTdx1hX97RgC4HHAhW//5bPX5fmTLZa62Haa2baa3jtNbnaK3fCnWZhDCD6OhoLu3c\nqdad5BwOB8tyP3Wb6Jb9doIBWb1Mt0Kew+Fg/ZavvErQNY0ucGiNIwi1hP76fAI1R4GIbKZM9qJm\nZR05rMjKsQP84/U3ar2qXU20hsnPzfBo32D3Co/yYfifq17/UVFRNE+tT+qyLcR/ZGfw+p/YfPwU\n4N+k6e/PJysryy8rG4Yrq///95UkeyHCjD/Gp9d0d9iz26V07dq1xuMEa4U8T8tdU4Iua3ooK/dj\nTUwLPl4AACAASURBVOtxYkAXjg/owjWN69Pvyx/5y3f7/JY0A/X5BGqOAhG5TNtBryZW76Bns9ks\ne4Vr5djh7Phr00nOHzPYhaJX+L///W/+MvaRgJX7uQMF/G/5Cr8kzUB8PpW//1Zb6MbK//8jroOe\nEMIztekkV9u7Q4fDwVJbbtA7uLVr1y6g5d539BgdO7oc8euVYH0+Mre/8ITc2QshfLo79HQIWJlA\nJKRAlru25ZUhcsLf5M5eCFErvtwd1tR+3r3rxVzfv19AO+4Fotz+6pgnveaFmUiyD1NWHmtq5djB\nXPG76xX+xM+H2bJ1a0A67vkj/mD1Zg/Eecz0+w8Fq8fvK0n2QgifuWv379ChA8+2aVTrde+DXW5/\n92bv0qULH3+yio/S20qveRFS0mYvhPCLip3N/N1WHcge54E6duX5+/v36smUGc9XGdJotd70wnfS\nZi+ECLlA9AoPxoQ9gSp35fH1g37bwcA+2eXNGFZdolaEhiT7MGXldisrxw7mj99fHdPcTUiT1f2K\ngEzY40+Tx4+rdv5+byfbqTjVrtl//4Fm9fh9JcleCOF3/uiY5i5h/ik92RTt/u54Mr5+0rixHi3m\n4+ru/8cffwxmOCJCSJu9ECIgarPufTiPUa+57F+hNZzo36na2L766iuvlvIVkS8i17OviSR7IcKD\nmSe+CZTqpsn9KL0tS1fn1hjb9f37BX0qYmFu0kHPgqzcbmXl2CH84vf3xDcTv9/n1YQ0wVh6t7Lq\nmjGmPT+rxj4NgNumgIy6sQGbijgcBPL7H4rvSrD4Ldkrpeoopc7x1/GEENbmLmG+8esJj9r9Q9nb\nvaZx/FZeotaM7HY72QN7ExsfS2x8LNkDe0feyAitdY0P4H5gB3Aa+AoY4WKfywCHJ8fzx8MouhAi\nkm3atEkP7ttHx0RH65joaD24bx9tt9s9el/DpHp69kXn6IKBXXTBwC569kXn6IbJ9fSmTZuCUPLf\nlZSU6JKSEpdlrC62QX1669kXnaP14IvPesy+6Bw9uG+fgJfPKjZt2qQTU+tp7quv+U8T43FffZ2Y\nGvzvSk2cec+nnFljm71SagiwAMgBNgNXANcB/wOGa60LnftdBqzVWgelEU3a7IWwDm/b/UOx9K6v\n3MXmjyWIq1N50p8BWb2YNnOW5Wb2yx7YG9u5dhiYePaGJQVk77yYVR+vDE3BXAh0m/1YYJbW+lat\n9Syt9Y1AP6A7sFopVbVRSQRcuLXb+pOVYwdrxl+x3b+m+EO19K6v3PVpcNcU8MysF/2S6L0Z528m\n/vz+OxwOPv0kF7LrVt2YXZfcFeb6rtRGjAf7nI+R8MtprVcqpboBS4AvlFIDAlE4IYSwsszMTBat\n+OSsu39/JLuKcxiUKfv3lAnjTVXzIfzDk2r8/cCDWuv/utiWCnwEtAGmAy9LNb4QItTCqRo/2MJ9\nWKO/STX+7+zA9a42aK2PAn2AjcDffCmAEEL4WzB7u0fycC0rePGZF0jMccCSAigsNR5LCkjMcfDC\n07NCXTy/8STZzwXOVUqludqotT4NXAu8Dvzix7KJalix3baMlWMHid+TNvuOHTsGfAnbUAztczgc\nrFxZuztNf61dECr+/v536dKFTz/JJXvnxUQNPUTU0ENk77yYNSsjaxniGtvstdYLlVIfAgOUUq2B\ng8BKrfXBCvs4gHsCV0whhKieu97l/+vYEfB+CdnqRgCUdXCbnpHGe/2M48/da3RwC8RUthVjK9Wa\nq7OzatVzftrMWfTL6glQtaf/HOuN88/MzGTVxysje7nhmsbmAecCPwOlFR7HgH6+jvfzxwMZZy+E\ncPLnuPpNmzbpQX16l49/H9Snd5Wx/cEaA19WnkDMGeDrHAYidAjwOPv/Ap2B2zDa71sDs4FWWuvW\ngbkEqZl00BNClPFXh7yKd+zuFp/xpoObP+4UA93ZMKLvZiNMoDvoXQ5M1Fqv1VoXaq2/w6iyP0cp\n1dSXk4ras3K7rZVjB4m/cvz+HFdf0zr0nvJXe76r2GyH84GzY6tNJ0Ff1i4IJat//33lSbJvilGN\nX9EOQAFN/F4iIYQIAU8vGqKjo+nfq6fbDm6XX5zJ1X16B2fCGq25rt9VIZn/X4QXTxfCkfpyk8nK\nygp1EULGyrGDxF85/mD2Li+7Y19qy+WRbXtdDu2LUsqj2gFP7sZdxZbVMKk8tqSYaAb9tiPsZsGr\nDat//33lSZt9WYe8kkqbGrp6XWvd2J8FrKZc0mYvhAD8N498de3j7yU04Ztt28rb8zcfP8WD3+5h\ny/FTREVFMbB3NjffPpKRt93GiQGd3bbnr/3iC6Y9/pjHc9K7i23c9/t5sFVDnr6weZWyWn3ioEhV\nmzZ7T5L9ZG8OqLWe6ktBvGX1ZG+z2Sx7hWvl2EHidxe/3W5nyoTxLFltA2BgdhZTn5/p1fC06i4a\nLrzgQoad/rXKhcA/dh5iSXo7pj4/k369enKsoIATA7q4TvbLvqJ+Ql2eat3AbQdAVyrGVqo1V/fO\nYskqGyf6dzLdLHiB7vBn5e9/bZK9J+Psg5K8hRCiNlzNI++tssVnpkwYzwPLbYBx0bDkzee4vFs3\nljrH1Fc0smVDHlqeix43lqdaN+CjQzHM3XukykXB3L1HaJ6WymNNE72ek75ibGXJLj4uzuv4qlPb\nJC2r6JlbjXf2ZmX1O3shRGBVTH41D7f7Cq3hRP9ObC8opN+XPzL9/GZn373vyuNowSm/3Y37c7hh\nbZO0J0MWRe0FeuidEEJYTsUhaTV1AhyQnYVSCofWdEyuy7LL2vHxoePUX7qZ+ks3M/67fXz0yUqU\n8unvtEv+mP/fX0vd+mvIoggcSfZhyspjTa0cO0j8oYq/uuQ69I4/0Ty1PqnLthD/kZ3JP+xn6vnN\nKLwmk5c7tCQrK4tLL73UL6MGyuJ3t969Nx0S/ZGk/TnPgSes/v33lSfr2QshhOW5a89/+bE7eWjU\nvUYV9iWtACN59/vyR+4+pyH/PnSS5XONO21/z0lfm34KZUn6PRf9EG5r0YAHlv8+r4AIf9JmL4QQ\nXqqYXKtrO3/uQAH/W77irDttf4wa8Ad/rmsf6Cl9hSGgQ+/MSpK9ECLUapMwAzlEzdNj+ytJ+2ue\nA1E96aBnQVZut7Jy7CDxR0r8vs5JX1383s7J749OfuCf/gOeipTff7BJshdCCB8Fc6remvjSs96f\nSbqs/0BhURGFRUUsWvGJ3NGbiCmr8ZVSA4C/YlyMvKG1nuFiH6nGF0KEnFmqsGtbJe+q6l+WvzWX\niKrGV0pFAa8A/YEOwFCl1AWhLZUQQrgWzCpsd/wx/K1is4K/lugV5mG6ZA9cCvyotd6ttT4DvAtc\nF+IymY6V262sHDtI/GaMP5hV2IGO318T7QSKGX//4cCMyb45sKfC873O14QQwtR87XTnj/P6q+/A\n5PHjmNoqlbtbNZTZ8CKI6drslVJ/BPprre9xPh8OXKq1frDSfvr2228nIyMDgJSUFDp37ly+GlLZ\n1Z88l+fyXJ5b4fmPP/7IX8Y+wvRWaWTUjQVg1+lintydxzOzXqRdu3Y1Hi8xMZE+V3SjoKQUgKvT\n6zPt/GYcP+Og0FHK4E07KSwqYs2aNSGP1wrPy/69a9cuAObMmRM54+yVUt2AKVrrAc7njwG6cic9\n6aAnhBBnq82EPW4Xs9m+n2WXteOCevEhXTpXRFgHPWAD0FYp1UopFQsMARaFuEymU/HKz2qsHDtI\n/BK/ze222vQdcDtP/vnNmPLD/qAPJXTH6r9/X5lubnyttUMpNRpYzu9D774LcbGEECJseJuQa5on\nf/Q3v7D2VAkr5rzvryKKIDNdNb6npBpfCBFskTruvKZpf5OXbuHL9evp2rVriEooIPKq8YUQwlSs\nMO68f6+ebnvzX92ntyT6MCfJPkxZud3KyrGDxB/s+M027tyf8Ve8iFlmy+UvPxzgie/21Wqe/ECz\n+vffV5LshRDi/2/v7oPkqs47j/+eacUaJEBjgQ0ubCPHYOPNWkiDwouxpJEtRhoJxWytK+vsVmHW\naznUZrcqtjDYQSBgMGGNIWy8MZiss8WLKWeDN454kQYJ1AKMQIN6ELDBL+yuiE3A2BCJGCHh6Tn7\nR8+Memb69fbtPvfe8/1UddE93fScp09rnnvPc865NVSduJbydeczDmJWna7rPvQe/fmLv9KxQ3u9\n7ASI9qFmDwBVxHnN96SptZf+/Secoh9sfTCVcWUZNXsAQMPq7aW/dcdOD61CO5HsUyrkulXIsUvE\n38n4k3QJ2wn0f953E1KJZA8ANVxzwzd0xYuv69Z9v2xp4lqxWKx75blOSeJBDNqLmj0A1NHqNrSb\nLv2ytuZLQ+Or+5brmhu+4X3S28jIiPr7lmnw5Gnb4774OpPyEqqVmj3JHgAa1OymOlX3m3/xdQ3t\n2Kne3t62tbXR9kU9iEHnkewDlM/nJ6+QFJqQY5eIP03x15rx/sAJp2rztu1Nv2c74k/TzoBp6v+4\nMRsfABKm3oz3LTvyiarhpyHRIzrO7AGgDbK8Rh9+cGYPAAnDjHckCck+pUJeaxpy7BLxpyn+uJbt\nlUtT/O0QevxRkewBoE0WL16soR079cAJp2reg8+w3zy8oWYPAB2QphnvSKZWavaz4m4MAGAmknw2\npPWgjWH8lAq5bhVy7BLxE3/edxO88hV/oVDQupWfVPfs2eqePVvrVn5SIyMjXtoSBckeAIAaJnZC\nXPvqCzrQv1AH+hdq7asvqL9vmQqFgu/mNYSaPQAANbRjJ8Qo2C4XAIA2SNLmSGyqE6CQ63Yhxy4R\nP/HnfTfBq9Djj4pkDwBAFVnZCZFhfAAAahgZGVF/3zINnjzzUsWd3CCJYXwAANokCzshkuxTKuS6\nVcixS8RP/HnfTfDqoYce8nJp4N7eXm3etl2HDh/WocOHtXnb9tQkeolkDwBIgYlNbfr7+71uapPL\n5VJRo5+Omj0AINEmNrUZXDCzZj60Y6d6e3s9t7AzWGcPAMispGxq4xsT9AIUct0u5Ngl4if+vO8m\ndFSxWNTWHfnJM/r8r/558rkL33uctuzIe6nhpw3JHgCQWIVCQWNjY76bkXoM4wMAEmvdyk/qH4d3\naf3J72IYn+vZAwCyplgsamt+p3aedYo+9dT/kaQpE/S+9Pc/1yN/9X2fTUwNhvFTKrS6XbmQY5eI\nn/jzvpvQcYvmzdHQWafqgVcP6JgtI5q3dUT3/mK/3nZK1Vp3nzizBwAk0pF96Usz8TefeYoe+uUb\n6jv+GP3li79SbuHvpnLNuw/U7AEADZmY9d7JBJuUfemTgKV3AIC2mdi9rnv27I7vXpeFfemTgGSf\nUiHW7SaEHLtE/MSf7+jvm9i9bu2rL+hA/0Id6F+ota++oP6+ZSoUCh1pQ/m+9FuHhlK3L30SkOwB\nAFVtuvTLGlwwXxcveJfmzOrSnFldunjBuzR48nxdddmlHW1LWvelTwJq9gCAiorForpnz9aB/oWa\nM2vqueHB0THNe/AZHTp8mATcIdTsAQBAVYlK9mb2aTN7zsyKZhbGZYwiCrluGXLsEvETf75jv+vI\n0rfXZjx3x89f08CKvo6f1Yfe/1ElbZ39s5L+laRv+24IAEC65oZvqL9vmSTNXPp2+w0+m4YmJLJm\nb2Y7JG1wzlWd6knNHgA6o1Ao6KrLLtWWHXlJ0sCKPl399RuYEd9hmbuePckeAJLHx6Y6OCJVF8Ix\ns22STij/kSQn6XLn3L3NvNdFF12kBQsWSJJ6enq0aNEi9fX1STpS18nq45tvvjmoeMsfl9fsktAe\n4id+4if+djyeuL9v3z61ijP7lMrn85NfjNCEHLtE/MRP/KHGn9Vh/Eucc3tqvCboZA8ACEtm1tmb\n2QVm9jNJZ0u6z8y2+G4TAABpl6hk75z7gXPufc65o5xz73HODfhuU1KV13RCE3LsEvETf953E7wK\nPf6oEpXsAQBA/BJZs28ENXsAQEgyU7MHAADxI9mnVMh1q5Bjl4if+PO+m+BV6PFHRbIHACDjqNkD\nAJAC1OwBAEBVJPuUCrluFXLsEvETf953E7wKPf6oSPYAAGQcNXsAAFKAmj0AAKiKZJ9SIdetQo5d\nIn7iz/tuglehxx8VyR4AgIyjZg8AQApQswcAAFWR7FMq5LpVyLFLxE/8ed9N8Cr0+KMi2QMAkHHU\n7AEASAFq9gAAoCqSfUqFXLcKOXaJ+Ik/77sJXoUef1QkewAAMo6aPQAAKUDNHgAAVEWyT6mQ61Yh\nxy4RP/HnfTfBq9Djj4pkDwBAxlGzBwAgBajZAwCAqkj2KRVy3Srk2CXiJ/687yZ4FXr8UZHsAQDI\nOGr2AACkADV7AABQFck+pUKuW4Ucu0T8xJ/33QSvQo8/KpI9AAAZR80eAIAUoGYPAACqItmnVMh1\nq5Bjl4if+PO+m+BV6PFHRbIHACDjqNkDAJAC1OwBAEBVJPuUCrluFXLsEvETf953E7wKPf6oEpXs\nzezrZva8mT1tZt83s2N9twkAgLRLVM3ezFZKetg5N2Zm10tyzrmvVnktNXsAQDAyU7N3zm13zo2N\nP3xC0nt9tgcAgCxIVLKf5nOStvhuRFKFXLcKOXaJ+Ik/77sJXoUef1SzOv0LzWybpBPKfyTJSbrc\nOXfv+Gsul/Qb59zdnW4fAABZk6iavSSZ2UWS1kv6hHPucI3Xuc9+9rNasGCBJKmnp0eLFi1SX1+f\npCNHfzzmMY95zGMep/HxxP19+/ZJkm6//fbINftEJXszWy3pRknLnHOv1XktE/QAAMHIzAQ9Sd+U\ndLSkbWZWMLNv+W5QUpUf+YUm5Ngl4if+vO8meBV6/FF1vGZfi3PuVN9tAAAgaxI1jN8MhvERRbFY\nlCTlcjnPLQGA5mRpGB9oi0KhoBVr1ukd3d16R3e3VqxZp5GREd/Nil2xWJw8oAGACST7lAq5btVs\n7IVCQcvOW6X86Ws1lj+gsfwB5U9fq6Ur+1UoFNrTyDaqFH8oBzNS2N99ifhDjz8qkj0yb8PGTXrz\nC4PSpy+WuueUbp++WG9+YVCXXHmV7+a1LGsHMwDiR80emVYsFvWO7m6N5Q+Ukny5QwfV1TdPbx86\n1FQNP2l1/xVr1il/+trSwUy5e27Vimcf0MP3bfbTMACxomYPdEASh8qLxaIe2bZVOv/CmU+ef6F2\nDm2hhg+AZJ9WIdetmok9l8tp2XmrpfvumPnkfXdo+aqBhs7QkzRU3mzfJ3XSXtR2hfzdl4g/9Pij\nItkj82762jWae9sV0j23SocOlm733Kq5t12hGwevbug9klr3r3cwc+yJJ+mouXMTMxIhJXOEBMg6\navYIQqFQ0CVXXqWdQ6ULKS5fNaAbB6/W4sWL6/6/7aj7x2lkZERLV/aXDkYmhvPvu0O6eYP0mT+W\nPn/55M/m3naFHtk2pN7eXi9tnRghmd7WqO1K2vwJoJ1aqdmT7BGUKMkh6clemnkwc+yJJ2n/yn8n\n/aevTX2h50l7cU0mLBQK2rBxU2m+gqRl563WTV+7pqGDNyCtmKAXoDjrVkmt6VbTSuy5XK7ppBxX\n3T8uleLv7e3Vw/dt1tuHDumtN9/UG6+8dOSMvpzHSXtxTSa87bbbEjN/wofQa9ahxx8VyT5g1E4b\nF0fdvxOmHMwUi6Vbh3TqoPFb3/kfiZw/ASQZyT6lJq57HFWSZpc3q9XYo1i8eLEe2TakFc8+oK6+\neerqm6cVzz6gR7c/2PGh43rx7927V8eceJL0iXdKH+uWvrhO+tH4QVwbRiIaPWiMY4SkWCzq2cJT\nQS819PH9T5LQ44+Kmn2g2IgluiRPCqs2AU63XiFdsF5zN/9lrAcozU64qzaZcO5tVzTUrjTMnwDa\nhZp9gFqpW6V9IxbfNbsodf841Yq/2hJBXTyonu13xz4S0eiSxIkh/lZHSHK5nD7auyQx8yd88P39\n9y30+KMi2QMZUe8g7o1XXtLChQs79vt2Dm3R8PDwjCF+M5ucTPj2oUN6+L7NTR2A/NHnP5eK+RNA\nkpDsU6qVulXSZpc3K201u7gnrqUmfif1rRqoOi8k6gjJ+vXrW54/kbYVKOVS0/9tEnr8UZHsA5WW\n2eVpFvdqh3oJqtMHcXV373vPSTr4h9e2ZdZ8+VLDZkYHWIGCYDnnUnkrNT1cO3bsaPk99uzZ41as\nXee6Zs1yXbNmuRVr17lCodB649osjtjbbc+ePW7u/OOdvnKL02Nvlm5fucXNnX+827Nnz5TXjo6O\nutHR0Zrv1Tdw/mQ/nX7m2VX7aXh42M2p8nsLhULd39WsQqFQMc458493lsuVHj/lpt4ee9N1zZoV\nuR1R+7+ZPkmyNHz/2ynk+MfzXqScyZl9wKKeHaG+RiauNXKWWWmJ5N4PnjNjieTEe511zjl668A/\nqeeu69W1/MgQ923f/K/60uVXxn5GW23C3c6hLTKLNGm4bZJ6fQOgE1h6B6+SvIwtqkaWhz3x+ONa\nsXpN3SVrjSyRrLb8bc63Nyo/tGV8uD2+/egnYpSm9tv0nyVpeWczS/ay+J1ENrSy9M77cHzUmwIf\nxk+76UPTfQPnJ6KEEMcw9+joqOuaNavmEPbygbWl4eTpz3/lFrdi7boj75Or/T6jo6Oub+D8mu9V\n7/lmNNNv1Yb4J0oKzXyeneiT3bt3J/I7CUwQw/jhSfNa01Z372tH7HFO3Ko3cW1Z/2o9um2o7j4H\nhUJBY25s5mueyk/erbf8Lb/1gYb3VKg3AbCZfisWi1q4cGFLs+ar9UmU/q/XJ71nn6MVq9ekYkfJ\nNP/bj0Po8UdFskfs6iWNpNVOJ5PYRwc09tDrsfyhr7Xa4T/8wb/R2Fj9EtSlm66RTl3U9tn1jR7o\nRJmHsGHjJt04eHWkWfPVDix+8pOfRIqzVp+8+dahRH0ngdhFHRLwfRPD+InTyBBvI8Opcc4Wb8QZ\nH1vq9OFep9ys0m3p+U53FSINc5ertNrhu9/9bmlo+8O9NYfWJz+nv/qhU8/MoXDNPsoNDw8751xL\nw/hLPr7UdffMn/H+3T3zp8xQb3QYPK7Z7nGWHpxz7uDBg+7gwYMV++TOO+906upK1HcSqEQtDON7\nT9qRG06yT5RGlzUlLdnv3r3bafZRM5Npz/FO33kslvaU15wnk9hdhYpJvHyJ3OTndNcep6XrjhyM\nnLvGWS43+Z71auO1nv/IGb9b/UBgad+UGOrOQ1hdfx5Co59XXN+RwcFBl+s5zsm6nKzL5XqOc9dd\nd92UPlm+eq1Tm5YJAnEi2QcoaWtNmzkTa/WsLc7Ya02U07lrYv1DPyOJTU/ic46ZPFt3rsLn9ORo\n6fYHX5zxOdXbM6HS88PDw05d1ZOcunJTYq/Vb31rzm84QdebcFcv2VtXV0N9Mjg46NQ9Z+aBXPcc\nNzg4OPV3nTtQM7ZG2t0pSfu332khx0+yD1CSvvDNnom1Oks7rtjrtVu53OQf+rb+vidHnXa+0fDn\n1H3MsVU/p0YS6cTzu3btKp3xVovfutzhw4frtmfu/OPd8PBwrLPdax1YLDr7nIY+71zPcVXfI9dz\nvHNufGSnK1ezXHLnnXcmapZ+kv7t+xBy/CR7eBVl2DUJu/fVTfZdXVPOtOPQ7KhGOz+n3nM+7jTn\n6OojG3OOaarf6s0NaKae3+oB4cGDB+seyDz22GNT509MH2k5rdd9ZNHiptqdlLN/ZBPJHt5FHZr3\n/cexZoJa1hf774uaxOL+nCbnKpx8mtPRPTPPaI/ucWecu7Tq7670s1qxnbF0edPfj1YOdBpJ9kv7\nByrPn9j5htOlf9FUu5O6bwSyhWQfoKQNZcW1gUoj4oy9k+2e0OrZehzxT85VuKtQSvan9U45o9Vv\nzY40fF1tbkArE+6mH1g0Gn+9Yfx68yeeeOKJhtrd6T33k/Zvv9NCjp9kH6AkfuE7NTQfd+y+SgpR\nz9Zbjb/qRMGuXOk252j3kdObG76u9DvKJ+TFuQKj0fivvfbaqhP0rr322rrzJ3bt2lVzAmOjOxjG\nLYn/9jsp5PhbSfbsjY/YpXVv8bS2u1lV94kvFksbzazoUe/Hlqrw8c/Etq99J/fJLxQK2rBxkx7Z\ntlVjxaJ0dI/06wOSpNy8+Rq89Ev66le/WrNNSx67W8///fN68/j3SRd8Xvr9/1ix3dv+7m8b3nMf\naFUre+OT7IEA1Up0Z/zwexrZ9cNYE9jIyIiWruyveDGeRrfPbUS1iwJ13/In2n7/vTr33HMbatP7\nT36/nj8o6YVnJDdW2snwsm9KH1o0pd0LFy4k2ccklIPtVnAhnACFPJQVcuzOxRN/q8voopQfapVL\nmilp1Io/jtUOd955Z+WNlmYfNV7mqLMfAsP4TWl2cmPW4m+GuBAOgGZUuw79o9sf1JIlS2peNCbq\nnvy9vb16+L7NU/bJd87FdgGiehcFqnTRn0pt+u93f0/64k0z9snXF2+Szl6lrrffmjISUWvP/RsH\nr44USyvqXZsiSVq9KBaaEPUowfdNgZ/ZA3FpdhldXJMX457FHsdlbKNutJSEfSPSuPyv06MiaScm\n6AGIW6FQ0CVXXqWdQ1skSctXDejGwatjq6+3Y9JezUl3P/xrPf+/n6tYo39k25B6e3urT16USmfs\ny47R8JNPasmSJRV/v6+6c7W5CuWxJU29z5r5DjMxQS9A+XxefX19vpvhRcixS52Pvx0JrJU/9LXi\nrzXp7rTTPqI9y/5t3YOLmgcMj/+1hnfuiBBxfCrF38nVDnGJ+h0I+d9/K8memj2AmnK5XGrOrqrN\nRcgPbdHI7l0N1fNr1eBvu/mmjsYj1a/BNzNXIUlyuVxb5oagMs7sAXjR7rPR8hGJZs8i4yhhtDoi\nUr5fgCQtO2+1bvraNTPakObh8E4tycyKzJzZm9k1ZrbXzEbMbKuZnei7TQDao92z2MtHJJo9i6w0\nS7/RxFMoFFpeYdDMLPU0nyHXWhVCoo9Z1Jl97bhJOrrs/n+WdEuN17YyqTH1Ql5rGnLszmUrsvVa\nQwAACYFJREFU/iiz2KPGn6YVBs1e4tfHNR7i1uheC1n6/jdLWVln75z7ddnDuZLGfLUFQPu1cgbd\nrE6cRW7YuKk0JD1tjf6bXxjUJVde1dB71KvB79395IwafBbOkNM0NySNElezN7NrJV0oab+kFc65\n16q8ziWt7QDSIWkrDOJ8H7adza5Wavaz4m5MPWa2TdIJ5T+S5CRd7py71zm3UdJGM7tMpaH8q6q9\n10UXXaQFCxZIknp6erRo0aLJJRn5fF6SeMxjHvN4xuNHH3009vefcrb9VOl5LSk9r8IjcmNHBipr\nvV8ul9NHe5do719slDbcNPX99v1Iy1cN1Gx/Lpfz/vnyOJ7HE/f37dunViXuzH6Cmb1P0gPOuY9W\neT7oM/t8wGtNQ45dIv4kxx/XCoNas9T/7PrrtH79+phbnh5J7v92y9Js/FPKHl4g6XlfbQGAZsW1\nwqBWDf7UU09tW/uRXYk6szezeyR9SKWJeS9Kutg593KV1wZ9Zg8gmeLeZpgaPCawXS4AJAxJGnHL\nzDA+Glc+gSM0IccuEX9a4m/XUrK0xN8uoccfFckeQNul6RrrScbniKhI9ikV6mxUKezYpXTFH8fW\nsdOlKf64lH+O/atXx/I5plWI/R8Hkj2Atmhmf3dUx+eIOJDsUyrkulXIsUvpiT+OrWMrSUv8cZnx\nOT63O5bPMa1C6/+4kOwBxC6t11hPGj5HxIWldwBil+ZrrCcJnyPKsfQOQKKk+RrrScLniLiQ7FMq\n5LpVyLFL6Yk/rq1jp0tL/HGZ8Tk+vjWWzzGtQuv/uJDsAbRFFq6xngTTP0f747V8jmgaNXsAbcfW\nsfHgcwwbe+MDAJBxTNALUMh1q5Bjl4if+PO+m+BV6PFHRbIHACDjGMYHACAFGMYHAABVkexTKuS6\nVcixS8RP/HnfTfAq9PijItkDAJBx1OwBAEgBavYAAKAqkn1KhVy3Cjl2ifiJP++7CV6FHn9UJHsA\nADKOmj0AAClAzR4AAFRFsk+pkOtWIccuET/x5303wavQ44+KZA8AQMZRswcAIAWo2QMAgKpI9ikV\nct0q5Ngl4if+vO8meBV6/FGR7AEAyDhq9gAApAA1ewAAUBXJPqVCrluFHLtE/MSf990Er0KPPyqS\nPQAAGUfNHgCAFKBmDwAAqiLZp1TIdauQY5eIn/jzvpvgVejxR0WyBwAg46jZAwCQAtTsAQBAVST7\nlAq5bhVy7BLxE3/edxO8Cj3+qBKZ7M1sg5mNmdl8321Jqqefftp3E7wJOXaJ+Imf+NG8xCV7M3uv\npPMkvei7LUm2f/9+303wJuTYJeInfuJH8xKX7CX9maQv+24EAABZkahkb2a/J+lnzrlnfbcl6fbt\n2+e7Cd6EHLtE/MS/z3cTvAo9/qg6vvTOzLZJOqH8R5KcpI2S/kTSec65fzaz/ydpiXPutSrvw7o7\nAEBQoi69S8w6ezP7l5K2Szqo0gHAeyW9JOlM59yrPtsGAECaJSbZTzd+Zt/rnPsn320BACDNElWz\nn8apdIYPAABakNgzewAAEI8kn9lPYWZfN7PnzexpM/u+mR1b5XWrzexHZvYTM7us0+1sFzP7tJk9\nZ2ZFM+ut8bp9ZrbXzEbMbHcn29guTcSe1b5/p5k9aGY/NrMhM5tX5XWZ6vtG+tPM/tzMfjr+d2FR\np9vYLvViN7PlZrbfzArjt40+2tkuZvYdM/uFmT1T4zVZ7fuasUfue+dcKm6SVkrqGr9/vaQ/rfCa\nLkkvSDpZ0m9JelrSab7bHlP8H5Z0qqSHVZrLUO11/1fSO323t9OxZ7zv/4ukS8fvXybp+qz3fSP9\nKWlA0v3j98+S9ITvdncw9uWSNvtuaxs/g49LWiTpmSrPZ7LvG4w9Ut+n5szeObfdOTc2/vAJlWbr\nT3empJ865150zv1G0vckfapTbWwn59yPnXM/Vf15DKYUjdg0osHYM9v3KsVx+/j92yVdUOV1Wer7\nRvrzU5LukCTn3JOS5pnZCUq/Rr/LmZ3T5Jx7TFKtydlZ7ftGYpci9H1a/zB8TtKWCj8/SdLPyh7/\nfPxnIXGStpnZsJmt992YDspy37/bOfcLSXLOvSLp3VVel6W+b6Q/p7/mpQqvSaNGv8vnjA9h329m\n/6IzTUuMrPZ9o5ru+1ntblEzamy4c7lz7t7x11wu6TfOubs9NLGtGom/Aec65142s3ep9If/+fEj\nxUSLKfbUqrPZ1HTVZtWmsu8RyR5J73fOHTSzAUk/kPQhz21CZ0Tq+0Qle+fcebWeN7OLJK2R9Ikq\nL3lJ0vvLHk9szJMK9eJv8D1eHv/vL83sb1UaEkz8H/wYYs9s349P1jnBOfcLMztRUsVNptLa91U0\n0p8vSXpfndekUd3YnXO/Lru/xcy+ZWbznXOvd6iNvmW17+uK2vepGcY3s9UqXSDn95xzh6u8bFjS\nKWZ2spm9Q9JnJG3uVBs7qGK9xszmmNnR4/fnSuqX9FwnG9YB1WpVWe77zZIuGr//WUl/N/0FGez7\nRvpzs6QLJcnMzpa0f6LckXJ1Yy+vT5vZmSoto85aojdV//ee1b6fUDX2yH3ve+ZhEzMUf6rSZW8L\n47dvjf/8PZLuK3vdakk/Hn/9V3y3O8b4L1CpRvWWpJclbZkev6QPqDRzd0TSs1mJv5HYM97381Xa\nSvrHkh6U1BNC31fqT0l/KOkLZa/5byrNXN+rGqtU0narF7ukP1LpYG5E0uOSzvLd5pjjv1vSP0o6\nLOkfJP37gPq+ZuxR+55NdQAAyLjUDOMDAIBoSPYAAGQcyR4AgIwj2QMAkHEkewAAMo5kDwBAxpHs\nAUxhZpvMbKzs9pKZ3WNmvz3+/AfN7Nvjl9MdNbOHfbcZQG2J2i4XQGLsl7RKpV28flvStZK2m9nv\nSPodlTZ9eUL8DQFSgX+oACoZdc4Nj9/fbWb/oNI++wPOuf+l8e1bzexvJB3nqY0AGsQwPoBGFMb/\n+wGvrQAQCckeQCMmkvwrXlsBIBKG8QFUZGa58bsflHSLpDckPeSvRQCiItkDqOR4Sb8Zv+9UuvrW\n7zvnOLMHUohkD6CS/ZI+OX7/Fefcyz4bA6A1JHsAlYw650Z8NwJAPJigBwBAxnFmD6ApZnaUpDUq\nbbhzkqRjzOxfjz99v3PukLfGAaiIZA+gWe+W9DcqTdyb8D/H//sBlSbzAUgQc87VfxUAAEgtavYA\nAGQcyR4AgIwj2QMAkHEkewAAMo5kDwBAxpHsAQDIOJI9AAAZR7IHACDj/j/SHNNNDhMd3gAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d2fca90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "classes = iris_df['class'].unique()\n",
    "colors = [(0,.8,1),(1,.3,.2),(0,.7, .1)]\n",
    "plt.figure(figsize = (8,8))\n",
    "plt.xlabel('P1', fontsize = 15)\n",
    "plt.ylabel('P2', fontsize = 15)\n",
    "plt.title('2 component PCA', fontsize = 20)\n",
    "for cl, color in zip(classes,colors):\n",
    "    P1 = iris_df[iris_df['class'] == cl]['P1'].values\n",
    "    P2 = iris_df[iris_df['class'] == cl]['P2'].values\n",
    "    plt.scatter(P2, P1, c = color, s = 50)\n",
    "plt.legend(classes)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three classes appear to be well separated! \n",
    "\n",
    "iris-virginica and iris-versicolor could be better separated, but it is good!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (b) Train with Multivariate Gaussian Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_train = iris_df[iris_df['class'] == classes[0]][0:35]\n",
    "for c in classes[1:]:\n",
    "    df_train = pd.concat([df_train, iris_df[iris_df ['class'] == c][0:35]])\n",
    "\n",
    "df_test = iris_df[iris_df['class'] == classes[0]][35:]\n",
    "for c in classes[1:]:\n",
    "    df_test = pd.concat([df_test, iris_df[iris_df['class'] == c][35:]])  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def classify(sample_df, valid_df):\n",
    "\n",
    "    prob = []\n",
    "\n",
    "    for i in xrange(3): \n",
    "        cond = sample_df['class'] == classes[i]\n",
    "        mean = np.mean(sample_df[cond].ix[:,0:4].values, axis = 0)\n",
    "        cov = np.cov(np.transpose(sample_df[cond].ix[:,0:4].values))\n",
    "        func = multivariate_normal(mean =mean, cov=cov)\n",
    "        prob.append(func.logpdf(valid_df.ix[:,0:4]))\n",
    "    max_prob = np.argmax(prob, axis = 0)\n",
    "    tf_number_error = [classes[i] != j for i,j in zip(max_prob, valid_df['class'])]\n",
    "    error_percent = np.sum(tf_number_error)/len(valid_df)\n",
    "    return error_percent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def find_error(sample_df, valid_df):\n",
    "\n",
    "    prob, label, flower = [], [], []\n",
    "\n",
    "    for i in xrange(3): \n",
    "        cond = sample_df['class'] == classes[i]\n",
    "        mean = np.mean(sample_df[cond].ix[:,0:4].values, axis = 0)\n",
    "        cov = np.cov(np.transpose(sample_df[cond].ix[:,0:4].values))\n",
    "        func = multivariate_normal(mean, cov)\n",
    "        prob.append(func.logpdf(valid_df.ix[:,0:4]))\n",
    "\n",
    "    max_prob = np.argmax(prob, axis = 0)\n",
    "\n",
    "    prob = np.matrix(prob)\n",
    "    for i,j in zip(max_prob, valid_df['class']):\n",
    "        if  classes[i] != j:\n",
    "            flower.append(j)\n",
    "            label.append(classes[i])\n",
    "    return [flower, label]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>flower</th>\n",
       "      <th>prediction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Iris-versicolor</td>\n",
       "      <td>Iris-virginica</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Iris-virginica</td>\n",
       "      <td>Iris-versicolor</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            flower       prediction\n",
       "0  Iris-versicolor   Iris-virginica\n",
       "1   Iris-virginica  Iris-versicolor"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Create dataframe of flowers and predictions\n",
    "[flower, label] = find_error(df_train, df_train)\n",
    "error_df1 = pd.DataFrame({'prediction': label, 'flower': flower})\n",
    "error_df1.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error:  0\n"
     ]
    }
   ],
   "source": [
    "error_percent = classify(df_train, df_test)\n",
    "\n",
    "print 'Error: ', error_percent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "100% accuracy with 0 percent error on our test dataset."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}