{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expectation maximization demo\n",
    "===========\n",
    "\n",
    "This is an implementation of EM clustering on synthetic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.misc import logsumexp\n",
    "import matplotlib.pyplot as plt\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate true data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "V = 100\n",
    "K = 3\n",
    "N = 1000\n",
    "M = 50 #expected number of words per document; generate from a poisson distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# generate the word distribution P(w | k)\n",
    "true_theta = np.zeros((K,V))\n",
    "for k in range(K):\n",
    "    true_theta[k,] = np.random.dirichlet(0.1*np.ones(V))\n",
    "# generate the prior over classes\n",
    "true_pi = np.random.dirichlet(np.ones(K))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# generate the observations\n",
    "true_k = np.random.multinomial(1,true_pi,N)\n",
    "w = np.zeros((N,V))\n",
    "for i in range(N):\n",
    "    N_i = np.random.poisson(M)\n",
    "    w[i,] = np.random.multinomial(N_i,true_k[i,].dot(true_theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = {'text.latex.preamble' : [r'\\usepackage{siunitx}', r'\\usepackage{amsmath}']}\n",
    "plt.rcParams.update(params)\n",
    "matplotlib.rc('text', usetex = True)\n",
    "matplotlib.rc('font', **{'family' : \"sans-serif\"})\n",
    "plt.rcParams.update({'figure.autolayout': True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.4599180434485888,\n",
       " 72.459918043448596,\n",
       " -3.2658688695177402,\n",
       " 68.265868869517732)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAACsCAYAAACgorNbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX1wXNWV4H9XbsvISP5oyQZ7IcGSybADDkiYGm+ZXZNI\nLc/WDt5g2WaKkNKyFZlsUvEykm0gwJAs8mADpraS2sSxU4tV+djFxuMU/CNZmIFdb01mkG2Ih91M\ngu1QmyEhGGFsbwSy7LN/3Pu6X3e/1x/qD6ml86t6pfdu36/31H3euefec64RERRFUUpB1UR3QFGU\nqYsKGEVRSoYKGEVRSoYKGEVRSoYKGEVRSkZkojuQijFGp7UUpQIREZOaNik1GBHJejz++OM55SvV\nMZHtT9e2J7p9vffwz8OYlAJGUZSpQcUJmIGBAdrbO/jhD59nYGBgorujKEoGKkrADAwMcNddnQwO\nruHUqX/LXXd1TpiQueOOOyak3enc9kS3r/eePybT+GkiMMZIWJ/a2zsYHFwDdLqUPmKxFzl06EDZ\n+qcoSjrGGKRSjLyKokwNJt00dSZ6ejZy5EgnIyP2uqbmQRYv/lPq65cCcOedt/Puu+fjeVevXj1R\nXVUUhQobIoG1w+zcuRuAxYvr6Os7CHwbOAHscedW+Bw82KdCRlHKQNgQqeIEjJ/6+qUMDz+Gtcl0\nAGqfUZSJQG0wiqKUnZLaYIwxLcASICoie1xaB3AWaPTS8sE/RLrzztvp69vkPlkCbIrnq6l5kJ6e\nvsJuQFGUgijpEMkYs09ENhhjtgAvAwZYIiIHjDFdwJCIHE8pEzpE8tbBjIzsAKC6+gHGxka5fDkK\nQFXV+9x88woaGurVyKsoZSRsiFQyDcYYsw54HUBEnnZp24FDLsspoA04HlhBADt37nbCxdpZRkd3\nAV+JX1++3EdDg9pdFGWyUEobzHKg3hjT7DQYgHnAsC9Pff7VvgQsdcdvC+yioiilJKuAMcZ0GGOG\njTFvG2O+7Ev/VQ71n/GGQM72AnaYNC4WL64DBoHH3PEhVVVfB/qAPmd32Tje6hVFKTK5DJF2YO0m\nHxljhowxQyLyBtkFxQfAaXd+FrjN/Y26tPkuTxrf/OY34+d33HFH3A/ipZeOYNe5dMY/v/LKx1ix\n4kUAVq36Ojt37mbnzt1qg1GUEvLqq6/y6quvZs+YLQ4EcMh33gi87s7fzlJuCbDFnW8B1gLNQJcv\n7ZaAchJGTc0igb0C4o69UlOzSERE+vv7pabmKvf5XqmpuUr6+/tD61IUpXi4322aHMhFgzltjPke\nsENEThlj9htjDpHQRMIE12ljzFk3NIpKwtC73BjTCpwVqwnlwUWg23fdzcWLl6ivX8q5c+cZG3sK\nT7sZGbFGYdViFGXiyCpgROR+Y0wbbkgkIk8ZY05hZ4GylfXWuRwISDucb2cvXboMjAG7XMoYY2Nj\nbjXvrvCCiqJMCDlNU4vIyynXLwAvlKRHGRgbA+giYdpZCfwIq7VcDdwbz6sL7RRl4qkoX6SqqjpE\nrgCecSmbsRrNh/HrSOSHfO5zt6uRV1HKSNGcHY0x14nIr4vVsYD6MwiYKCL/noQGswT4HvBZd/1z\nenu/wSOPPFKq7imKEkDBzo5uPczbwG63JmZtUXuYUx8uA89hvabXAP/VffIV4CtEIjUsX7683N1S\nFCWEnDUYt/5ledh10TqUQYOZObOesbFn8YdksMOlE/FrDdGgKOWnGOEahrNcl4Egm/TZsvdCUZTc\nyMfZ8bQxZgDrFR0DcB7RIiI/KEXnUpkzp4rh4U2+lE3AKFaT0ZkjRZls5DNE2upOhYSbgLf09umi\ndSjDEGnOnE9x/vwG/EbeSGQvMAOAO+74LMbMAQqLyeuPOaOzUYqSnbAh0oRuhRl0kMFVoK7u2jRX\nAZjr/vYIzCnYVUBdDhQlfwhxFcjnhz+UcryeQ5kd7m+XL60DaPWnpZQJvYmmpmUCVwhc444rBG5w\nwmZtmvCJxdbm/aBisfHV09/fL7HYWonF1qpAUqYdYQImZxuMJM8gtWGDRWWjy01n3+/Ktbi6Dhtj\nGo0xzZIS0S4T7733f4FqoNelbALeDc1/5kygs3bRSY20d+RIp+5ooCgw/iESPi/rDHk6Uq63A593\n5604b+uUPBmkZDRgiDQ/cIgEDdLcvDJvSTyeIVKQ1hONNqk2o0wbKFSD8UWlg9wj0UWd53SLWENw\nESLapfWMaPQJ500dA1506Z00NJzOVDCQ1atXc/Bgn8/IOz5NZHh4AYODa1SbUaY1eU1Tk5hBOgU8\nma2AJHYSiDlBAwVEtKutvcyFC8nT1DNmfMy5cxFELlFVdYjLl78DFDZlvXr16pwEgjfbdObMe0Qi\nPc4ZE6yP1I+A1Ro2QpnW5GODyct72q2RGRaRA9jIdY0UGNHuk09GsOtevJgw57h0abarGuDnLFr0\nn7jpplvGrXnkysDAAGvW/Dmjoze4lAvYVcVnsSuNVaAoU5dcI9rlsw6mFdjvLgVYLyKvZMk/JDbU\n5nbgv2O1l+UisscNuQYlJehUpnUwxlQDNXjbw1oj7xjwXXe9GfgYkfM53VMhtLTczvHj/0iyZ/cf\nYWMF3xtP1y1slelAMbYt8cfmbQT2YXcOCETsTFGHMQZs8O83XEcKiGhXB/h9kQAeTbnuphy8887v\nXLsJmw/8FKu5dBKNPsGtt95cck1KUSYz+QiYYRH5CEBs6MysvkhueJSaNu6Idm7hcBYu5V/tOJg/\nfzbDw56zJVgN5krs7gY/4ic/UcGiKPkMkfYDf4fdKK0FuE1ENhS9QxnjwcxGZCbJQ6RREuEyu4HR\ncQ2R8nUPaGm5g+PH78Pv2V1X9xgrVtzGqlUtvPbasZzrUpRKpxhDpK3AenecLIVwycbMmXMYHb0d\neMKlzAHuJjFMuY8ZM/bmXe94Fso1NKTPsK9YcRs9PRt10Z2ieAQtjgk6yLJNSbEOMiy0M2a2wCyf\nq8BM3+I7u7hu0aLGXNcGxUlfKNcTulCut7dXotEmqau7ViKRuZK6IG+8rgZevdFok/T29uZ9D4oy\nkVDoQjvgsDFmCBuuwRNODxVV2mVB5GOglmRXgVr8htabbsp/cV0yA0Afw8PPMDiYrIFs27aNRx99\nCv8QralpJ42N18eNud4wKx9S6330UbvWR0N/KhVPkNQJOrBOiklHrmXzORi3q8BeiUTqpbl5Zd5L\n9JPdA1aEaiDRaFPaZ5HIQonF1kpvb6/EYmuluXmlVFcvSOpTbe2ijJpJUL3RaFPO/VeUiYYiaDAP\nY5f6xzHGbATulxIGAU8mzCC9C7jA2NgIx493AfnZPvzuAT/72T9xPsVGnMlpcmyskcHBNQwObsJu\nqbKG6uoHaG5+jnPnPuLkyREuXNgJqGaiTEOCpE7Qgf0Vt7rzNqzjYgs5OD3mc5BRg5klMNtpGSvc\n+RUuVEOwFpCvNtPcvFKgIdBpsre3N82hEvp92tTauA0nElkoUC/Q4dLXxm07qaTXO0ftMEpFQRE0\nmEYROewkwMvGmK0i8pBbSFcWjKlGZAZ2FwGAB7AredcQFLYhm8Nh0NR0Q8NVwArCnSY/wS7u+xhY\nRbpLgLXhjI15wcj3kLDZbGZ0dGZaPz2N5tln7exYd/dW1XKUqUGQ1Ak6sG4Cm7Eb2G8BDmFDLuQS\ntmGL73zcAaegLsAG0+TO+9M0D792kTqTExaWIVO4htraRSltzBEbJsJ/7rfhpM8oNTXdUsB7QlEm\nJ4RoMDnvKiAi64EFWFtMPXY9DL6/gbjgVF6Q8HjAKXfdnGv7lvS3v9UowFuiD91Eo0+QzeFw587d\nbq1KJ2DXrdxzz9fYuXM3jzzydWKxF4nFXkzSfD7++BJ25W6nO74N7HXtzQP2Ar8EXgptt7GxMX6+\nbds26uuXUl+/lG3btuXyABSloshniISIPJiSlMtyf79l9m7sGAJsyIc27MrgHPkEq0R5bAb+gLer\nAPRRVXWJn/zkv7jFbsuA3EM3JIZUwQ6KNTWRNANwXd0sZs26iPWi9q8wXofdeTI5vMSqVTZ2uk5N\nK9OCILWmWAfQ7P4ekoSh+BZ33gpsDygTqoZVVc0RqPUZeWvTjL6dnZ0ikj1GbupQKNuQyrY/N20Y\nVlU1V2bM8Kal/UO3qDP0rkwy8maa8tapaaVSoQhG3vEQDUgbt1V41qxZjIx84ktJH+Fdf/31QPag\nUf6p6aNH32R4OPOQCqCqahaXL38JvwG4quqHgXkjkQif+9ztDA6uIXknytOB+RVlShIkdYpx4LQX\nSdZgtpOY6l5HnjF5a2quDtAUrs06DZyN/v7+pMVx1dULkrQeTxtatOhTadPJnZ2d0tnZGZieyWCs\nU9PKVILxajDGmJMk7Ch+jeSkiNyWoWijixtTj43N2ww8j40hcxhroBgMKhgW0S4SCTLy/j/sNDWE\nTQPnxkUSXtkX46nJjpBrqKq6H2O2YIzhi1+8i71798bz/vjH1r7iTw+L76tT00olk2tEu3w0kn0k\nbCotwL4cy3UBvyJhe+mioH2R/DaY2W5quLBp4EwOipk+072QFMVCEWww88TtYSQix4wx87IVcHn3\nYFeb+a9hXAGnLmHtLt5Cu03Ar5Ny+KeBi0GQm8CZMx/oXkiKkgP5CBhjjPky1ps6VqL+ZOT3v7+A\nndbt9KV2401Tj3cngZ6ejbz22pcYHbXX1dVb6OnxjLdjJE+Nb+Ltt+dyzz1fY2Tk3nhfdPcARUkn\nHwGzGbgH2AAcI8sCu1Iwc2a6faWurpYVK+ysTmHxbxM2mLGxc6xfv5GZM2cyf/4srEvAE9joeZc5\nf94LeLUZK2ttm0ePvkl7e0dSRDuNbqdMa4LGTUEHkyDgVNhsTaEk21lSZ3dmS2YHxxXxvqS7DSSf\n57JLpKJUIhQx4NQgdi2LiMjDRZd4GXj33fNYjcHTIGIurZg8R/IwbBfW5uMflu0msWbmfdefLhIB\nwMGulXkm6VyHUcp0Ix8BM0jytLKEZSwtdwLeHnB9JBa9jZ+eno0cOdLJyAgkfJsy8a5r29vBcTew\nrOB+KMpUo2Q7O5aCZEFgjbqLF/8p9fVLAejuvm9ca0n8q3pPnZrDyZN+/6GfA3/hu+6mtnYW11//\nHCdOXGRs7HcE+RxZjaYv6TwS6aGn58d5909RKpWS7ew47g5l2LYEkmO4LF5cR1/fQfxOhr29hS1Y\nGxgY4M/+rIOxsfku5Qzwb0gobEuIxU7T07PRt3Xsb4A/8eUxRKNvMH9+HSdP/h9shAuAn9Pb+w1d\nUKdMOcK2LcnH+DoEzHXnjdhtYctq5BVJXtxWV3dt0RwGvXqDnBCrqurTlvsnG4aTY9F4edShUZku\nUAQjb947Oxab1MVt8DclqDc9Mt7NN99EQ0PyVHjy7gG6XayiBJGPgPnIGLOZxM6OZ0vTpXBskKh7\nSRh2P0+q7aO7e+s46/WCT12N3bzeUlPzIE8+mRAYAwMDtLd3cObMB1RXPxBfnOffLtbLY7eX/SoJ\nH6ef0939jcA+5LuzpKJUBEFqTdiBDZW5iwAv6GIdZBgiBQXkXrTouoI3LMt147VU7+jq6gVp26Qk\n5+lJWkMTidTnFJtG18solQYhQ6R8fvhv49tZIMcy67COjbt8aeOOydvcvCrNptHcvKrgh5PrDzyX\nXRuT8+S2y+N4d4NUlMlCmIDJJybvUuwsUswY87ox5vlM+d2sU6vY+LuNxpjmQmPyBu0HHZSWL940\ndVAc3mx47gEDAwMBnwY7SirKtCFI6gQdwHUkdhM4ik8ryaHskPu7A/i8O28lz4BTEz2USA+zme4G\nkJznhrQhnbfH0mS6L0UpFIowRLqMnarOZ4g01wmlDnddUExekYmPwZI8nd0TOKzxNrK3MXmTN14L\nG/pM9H0pSiGECZh8ZpGi2F0A1htjtgOnROTuLNrRR8DTxphDxphjLjlrTN6wiHaQPdZuqfHab2/v\nYHAw3T1gYGCAbdu+45tK91byLssYTmKi70tR8iHXiHb5rOSdg912pB270O51EflKhvwtWKl23Amk\nD7DhMwdF5LAxZh2wRESeTiknufZpIvCmk8+ceY+33volo6O2+9XVD3DjjTfzzju/YXj4CyQcHft8\n62N0+lmZmoSt5M1Hg3kFG1N3q4icziF/KzZuDNhdyf4eG6wqa0zeyUrqQj9vk3sY4623ZnL8+H0u\nZ3KcmFtvvZlDhw5MRJcVZULJx9lxeZ517wY2uMDfH4rIXwMYY5a7GaazIvJGnnVOKMkL8mB0lPgK\n39HRLpJDOnwT+B2RSA9/+7fV1NcvHbczpqJUKiXbF8nZX/YEpBcQk3dyEjT1XFf3LgsX/mdOnhzh\nwoWdgO7eqEw/sq6DcXF4McasLX13Jjc9PRupqXkQG4ahzxltN5KI29uHFydm6dJr+fDD88B38e9l\n/eyzz01Q7xWl/OSy0O5uY8wh4AduNsg7glaWTWnCFuQ1NFyFFSAvuqPTpSnK9CanWSRnR3kQuzNj\n3FIsIqeK3qFJPosURKrxt6bmQQ4e7GNoaChpg/tixKtRlMlI2CxSztPU5aISBQyEe0Nv27YtPixS\nI68yVSlYwEyWiHaKokw+wgRMzs6OWD+iJSISBW4DnipW56YaXjyYcCdIRZkeVFREu0pAt5RVlAQV\nFdGuEkhdjKd7ISnTmXziwax3+de76w2l6pSiKFODvFbyikhedhdjTJc7bRKRh1xaB1b7afSt6p0y\nBO3dFOZBrShTnXxmkbZjg0z9Osf8rdiQDqeNMfuA7wPDWMFywAmfIRE5nlKu4meRNIC3Mt0oxjT1\nOmAjNi7MLhH5QZb8XWB9j5xwOgk0kQjX0Aq0VFq4BkVR0il4mlpEXhCRdmwYhg3GmMvGmO8ZY64L\nyb/HNwRqwUbDm0dyoNrCA+oqijJpyVnAGGNajTG7sPF4T2HjuhwgsRN9WLkW4KhvKJQ1op2iKFOD\nfIy864H9qVHs3PAnE60i8rA7P4sdYgHMJyjsPplDZiqKMvGUImRmh4gc8F2v9YJIZSizUUR2u/NW\nrJF3ubPLbMHaY95IKaM2GEWpMMZt5HXTyndjA36/7PuoRexeSWHl2oB9WKESBdaJyCvO+HuKkGlq\nFTCKUnkUImDmYo2xW7H+SB5x14FiogJGUSqPQoJ+bxSRp40xZ4H7fekCPBxSRlEUJScB4wWVGsIK\nFUVRlJzIx8g7F2tLGcYuuNuf66revDqkQyRFqTiKEQ9mD3bDtR3YtSz7M2dXFGW6k4+AmScih7Gz\nP0+hC+YURclCPgLGGGOeBI4ZY5qxy/4VRVFCyUfA3I+1vzyJdRNYX5IeKYoyZdBdBRRFKZiCjbzG\nmC5jzNu+41c5ltuRct3hHCe7wsooijI1yHeIdKuILHXH9dkKGGM2Ah2+6xYAZyzG2XIURZmi5CNg\nhvJ1DXCOjv7dHzcAH7rzU1j/JkVRpij5hGuoN8YMkXB4FF8YhlyZhzUUx+vMs7yiKBVEPgLm+ynX\n47XE6voZRZkm5CxgROTl7LmyklPAKUVRpgY5CxhnkN2PtaHsA05mCzgVwPPYNTSHgSXAYFAmjWin\nKJObXCPaISI5HSSCdu/zrnMosw5rc/myL60LGzi8K6SMKEoqvb29Eo02STTaJL29vUWrt7+/X2Kx\ntRKLrZX+/v5x58kn31TE/W7Tf89BiYEZnUDxCZhDuZbN51ABMzUp5MfX29srMEdgrztmS1PTLWl1\n5SsImptXSnX1gni9NTVXpZXr7++XmpqrMubJN99UFELFEDDfB3YBh4DtnqAp9qECZuqR648vjGi0\nyZUVgX6BhrS6xicIVvjqFYG9EoutTcofi63NmifXfIU+h8lMmIDJx8h7v1s4Z7D2l4dyLatMb3bu\n3M3IyA6gE4CRkRPcc8/XuPXWm/PY+fIl4AngvKvHq4v4LpojI/cCL8bP/W14/Th69E2Xr9Pl3YON\nBgtwPXAVkNid8+jRN4E1ST05evRN2ts7WLWqhddeOwbAmTPp8xWp+Wzb/udg+7R69eqk3UAXL67j\npZeOANDdfR+PPPJI4BOpiB1Eg6TORB6oBjPlSH67B2sgmWhra0sZIjW4eqym0Ny8SpqbVybVa89v\nENgr1dXzkoZCifKp9c6Rtra2FE2jJy2PTUs+T28jKF+wxpSq2aSWDbI5TTZtiBANJpeg3ydJrHmJ\n+j46KSK3FVvgqbPj1GNgYIA1a77E6OjT2FH2V/De4tBHLPYihw4diOdNfSvX1y9lePgLwGlXZgnw\nv1w9m5gxYyZQxaVL/4rEV3UJdl7iVeBfACtTyj8P/AF4Nqkv0EM0Oo/h4auB37n0WUQiZwAYG0to\nOVaZfwO4GVjCokX7ef/9jxkbu+Ty3OBr7zQ2EOS9wDMAVFdv4cYbP8M77/yO4eFbUvp+GruvYR/R\n6BN88MHbSc+0vb2DwcE1oc/R/+zLoeWMO+i3iDS5CvYBT4rIcedTpEMkJQ8uYoXLb9I+8YYXAwMD\n3HVXpxtGwJEjnRw82Mfo6B+wP/5nXInNwEzs8CbCpUs7gRPu+tu+PH/kzs8GlG8APg7op3E/9kFf\nXZsYG/N8c/dgJ0Jtuj1fBvwFv/3tJ8B3fZ/9c+BOX19WA53U1f0lS5cu4a23LnL8eFBd/r7DxYuj\nAf3MTtjzLOtQKkitCTpImTVKvS7WgQ6RJg2FTs966dZI2yGwVmCVwEw3TGkQmClVVXMlGm2SpqZb\nUoYQKwWi7pjtO18gcLXANW4oIa5uf9kOgfkuf23KZz0CnxL4jKvXPzSpdWVWujrXuvN6dyQPcezn\n3vk1Ak3u8Pq+0LVV7z6fLbW1i1IM10F13eT+RqWpaVngM882RMrVQF0MKNTIi41o92WsL1Ks2IJO\nmVzk+vYLywf40teQeEOvAY5jVftlwCYuX76J4eEuhocf8NX877BaybeBvwTGsMMZXF3LXH2bSf86\nbiNZA/HXO0CyNrMJOykacW34NZCr3bnXj12hz8saoc/h13psv64jWbPaxIULMzLU43HWtTcW+Onq\n1as5eLDPN/wps2aSK0FSJ+jALrLbQWKaem6uZfM5UA1mUhD09otGm9K0FJuvx/e275FIZKFEIgsD\ntIameJ7kc08zmSswy73poz7tJBrwtq935f/YfV4vMEOsIbXBaRB+bcIz1KYbWhP5eiRYe+qQIAN1\nsjF2fkgfF6bUu9fVOU/gCtefFe6+b3Lnc8RvxK6paQjUEP2LDzs7O+N52traJBJZKDNmLBBjauL9\nra5eUDJDMEWYpj4LPFh8EadUCsPDCxgcXJOkzZw58x7wP/DbN8bGGkgO2ZyqNXg2kDXufCHWpPdV\nYBbQ68sXI9iv9npX/mUSNpE9WMMvJNs0NgGfYDWC0wF1LXDlvPb+GwmtxatrG/AIVvN6FLgGq138\nCGuUDeujv15Pw1jszr3+ngD+t8sH0O2r4wQjIxedQTehIQ4NDfHoo0+5Pp6gr8/Tkl4C/j6l71ux\nG4JcDOhjiQmSOhN5oBrMpKC/v18ikbm+N+wVAte6N31HfCzf3LwqQFP5lNMA5mXQGlal5F8r3rRy\ncj7vjR4+TZ3QnjLZNBa6o07Sp7P9dXnaU4dPs+lw6UH5r3XtLMzSxxW+dM++4/U3qO/XuDJz0z6L\nRpuchujZtfzaV5Btx9P2ekJtMIWuMKYINpiiYIzpwA4wG0VkT7nbV3JjaGiIsTEh+Q37hPt0E//w\nD1Ha2zv4xS/+0VdqAHiOhK3kq1il91JAC/WkazavuLY8TgC/BGYAf+zaP4fVIvz2hg9IDy10AngT\nG1DxPeAC8FnsLNafYBfZvRlQ1zXuXpNnkWCUxOzO11z6LcBcrCb1ClZb2uru9w5gtzuWAO+7NjtJ\naFEnXP+GXB4/nxDG8PBZ7BS510dPE4y5fqYSwbOD/exnc2lv70iask61o7322p9z440309BQX/jU\ndpDUKdUBtAAd7rwLaA7Ik7f0VIpP8ixH0BvWszl4s0H+mQ9/vpskfLFakGbjaT1hZT4doIHckJI/\ntaw3A5VaV2o+T+tYHNKvlZLep3kuT4Pvs46AfB0p50F5PHtOYpGgTZ8bkC/o2a0IeT7LJNMCx0IX\nQopMHg1mA9ZIDImQmcfL3AelKFyJfRs/AXwB+3b+fUC+ERIairckP4Z9i78fUvcTWCX32yQWknnl\nL2FnmF50aZ3AD119MRLuBP6yqYv7wNo57vCV8epaDdQE9KkKu/DOX+8JrMa2FKtxfIeEC0JQ399w\n7b3h0lLzPOHuo9P3F6wG5T2TGPZ5dpDMCSIRG512bGwVyc/nNFaTeoYgN4Vkcs2XG+UWMBoys0Lo\n7r6PRx/d5K6WYL/kHt4ULMB9gGds9BsqIXnB2DKsMfc97OKzTuwQ6V5/q1gV/zGCp4QbsStyUxfU\nzcFb9WqNra/ncIezUsp4q2yXYYc+qfd7OaW8N7zzT52fIJxGEgbfVCO4xwISw50fuTSvzsfc382u\nba+/Xp49jI35h3R+A3cXduVzMD09GzlypJOREYB3M9zDOAhSa0p1YL81ze68FdgekGf8er1SVPzT\noG1tbUlTov5FXpHIlVJX9ymprV0kVVV1kjAMzxaS/HEWiF3I5lfhZ0tieta/IC51SthvNO2RhBF0\nji9f2DCqNsuQxevjlWKHPPPELga8xh2zAoY1+QzvUg2+89zn/mHbfLFDoWt9z2GvBE/Re23PlkWL\nPhO4aM9bUtDb25tTaIp8QlgEwXh9kYqJMWY7MCgih40x64AlIvJ0Sh55/PHH49ca0W5yEubj4vdC\nTvjznAf+JfACdjjh9ysywP8EbscadDeTGBpsBn7qyn+JxFCrDztksD5Adoh0u6vrb1ye+0j2PXqO\nRDjoC9ghHkAtNsiil++nrt5fYLWtGVgNxzNUf4wdQgnpfkz/0eUH6yZwFdaQ/IWUvndjBw9fSunj\naaxWtdl3T6nPxPpLRSIz+OIX/zV79+4tul9SLvlSI9p961vfQgJ8kcqtwTTjItkBW4BbAvLk/6pV\nJh3hHsmZDKCzAz7rCSnTE1A+zHg7R9JdAlbmUG/quT9/usHXmJkB9xjUl5Uh6en3lOqlXUhQrFLC\nZNBgwO5o7LTjAAAG5UlEQVQQiTXwBk5T5+pN/eqrr06oZjOR7VdK22ExTlpalnDs2Om08+7u+wB4\n9tnn0j6bNesPzpnwIrW1woULswBoaprHh26nrdHR97hwoRqAmTPPc/HiHADa2pp56623XHmorR2l\nutp6RM+fDydPngWguflTNDQsCe3vxYsjLFxYE88fjQrDwybexubNmwPv99OfnsOJE/8EwLJl/4x3\n3jkHwJ133s67755Py+9P98eyGR5+n23bHgnUKErtNZ3t/x7mTV1WDSaXgxw1mMcffzxPGVtcJrL9\n6dr2RLev9x4OIRpMPjs7Koqi5IUKGEVRSkbZbTDZMMZMrg4pipITEmCDmXQCRlGUqYMOkRRFKRkq\nYPLAGLMj5brDGNPqpt6VEmKM2eI71+deIVSkgJmIL5jbE6rDd90CICKH3XVzCdvucsd2X1pZnoEx\nZp1rZ5cvrazP3xjThnN+KvNz3+H+dvnSynbvxpgW115Z23ftXjbGvO2O74237YoTMOX8gvkRkd3Y\nBYIeGwC3xCvuGV50jDGtwMtiFyU2un9ws+tTSZ+Ba7vVtdNojGmeoOfvNxTeTRmeu6PLGPMr4CRM\nyHfvIRE5AMwr87OfLyJVIrIUWA88Nd62K07AUKYfdg6UyzO8kcQ9nnLXd2N99720kjwDETksIv/B\nXUZF5Djl/YFjjGn2vtSOuZTPI79LRK4XkVfcddm+e85X73UAEXnaPfuytJ/yvJeLyOnxtl2JAmYy\nhXxIXxpdZERkjyRcKlqw4c/mYcO4eZTsGRhj5jr7x5MuqZw/cEje7C/erRK3GW/baYye/aec373l\nQL3TXCaifU+Dfb6QtitRwED5vmCZOEviyz+f5B980XEq6lH3JoMyPQMR+Uisx/v9xhgvrmNZ2g7Q\nXqCMz90J98PYH3qr161StRfAGe//7ULNlrv9mIic813n3XYlCpiy/rAz8Dx2uALW136wxO21isjD\n7rwsz8AZ+7yx9jFgXbnadjQ6w+JGrDbRTJmeuzOqez/qD1yb5bz3D0jEcjgL3Fbm9sFqzB7jarsS\nBUy5f9hAfEy83G0+h+/N0gqcFZE3MpUvsO2NTovwq63leAatJL5U87DGzrI9fxE54Iycgh2aSRmf\n+ynsnihghwOvU97v3gu+tuZh9yIpW/vGmMaUpHG1XZErebOFfJhKuCnafdjxbxRYJyKvlOMZGGPm\nYo17uHYedunT4vn7NJglIvKMSyvbvbu2hrGG1rI+ezcc3uoz8o+r7YoUMIqiVAaVOERSFKVCUAGj\nKErJUAGjKErJUAGjKErJUAGjKErJUAGjKErJUAGjKErJUAFTQRhjThaxrg5jzBbnzJi6k3o+9Qx6\nq5tLQTHvuVCMMRv9ga+U7KiAqSyKtirSLcN/GrsM/u7x1OGWk4uI/KBY/QpgUqwENcYMYvdWnxT9\nqRRUwExijDHznIZwyBizL+Wz7xtjhtzR7NLajDH7jTG7XHqXS2/01+O0lnUuQt5WoM1pNLs8r2Hf\n5157+135IZ9n8YNY/6y1vnxDxpg5ru+XjTHXufS3w+pxbX3fRU+7OeyeU+7fq8e71w5fHUsytOMP\nvTnk/m51Zb38aRqdiMSA+5kcnvwVQ2SiO6Bk5GHgeRH5gRMiXjSxjVjNYbkxZh42RsxSV2aJiKx3\nfkRHgT3YUJ9DIvKw+7FFSbyJd2CDSR0wxpzFRjA7DGx0B8aYrcDficgzrt7Tro7trr2/9vX5eRJB\nqY4B640xx4CX3Y87qB6AW0VkqbFhKtPu2Y/L49XTSiL4kVdHWH8zaR/XiUi7q/9t7C70SoGoBjO5\nacZ59PriwADcig1lsA/YTSLSGL78H/nSdgPGGHMIK0D8gYPiuNgnbe5HGRWRX7uPGoHjAfUGvc1f\nANqxP/oHsbF0Y8B+oCmkHkh454bds58lOMGTEi/G837Op79gBc/LvutTnualFIYKmMnNMVKCXTuG\ngGMiskFENmC9rTOxAasVtGO9YTdmyPsyVhjs8qWdxMUG8ewuYYVdeMVGrDZxGGvj8eL6ZqrH+/GH\n3bOfeMhG522eSlg74vqTWs54bToafcJVKQAdIk1ungT2G2PWYwP+nAQbac3ZUg5hVf+/cvmF5B+t\ndz7k6jnr0tZjtSDBajMtxpi1bqizGzucao9XIvJ0SnvrA9rw8zo2KFHSeYZ6/P0OvOeQ5xLDCqQk\nMrTzMjYy3yFXLqnvvvxbA9rMdL9KCBquQUnCvdk7/HFAJhue0VZEDjst50kRWV1AfVuwwaumbGyb\niUI1GCWOsVH7HsKGxpzMeBrZg+4605AvV/RNWwJUg1EUpWSokVdRlJKhAkZRlJKhAkZRlJKhAkZR\nlJKhAkZRlJLx/wFhLEXhloepmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f87734443d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(4,2.5));\n",
    "ax = fig.add_subplot(1,1,1);\n",
    "\n",
    "theta_comparison_0 = (true_theta[0,] / true_theta[1,] > 2.5)*1.0;\n",
    "theta_comparison_1 = (true_theta[1,] / true_theta[0,] > 2.5)*1.0;\n",
    "\n",
    "ax.scatter(np.dot(w,theta_comparison_0),\n",
    "           np.dot(w,theta_comparison_1));\n",
    "ax.set_xlabel('density of word group 1');\n",
    "ax.set_ylabel('density of word group 2');\n",
    "ax.axis('tight');\n",
    "#fig.savefig('wsi-synth-scatter.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now do EM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def initialize(K,V):\n",
    "    theta = np.random.dirichlet(np.ones(V),K)\n",
    "    pi = np.ones(K) / K\n",
    "    return theta, pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def eStep(w,theta,pi):\n",
    "    # first term: log P(w | z=k)\n",
    "    # second term: log P(z=k)\n",
    "    q = w.dot(np.log(theta.T)) + np.log(pi)\n",
    "    #q = q - q.max(axis=1).reshape(q.shape[0],-1)\n",
    "    # log-normalize\n",
    "    output = np.exp(q.T - logsumexp(q,axis=1)).T\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def mStep(w,q,alpha=.0001):\n",
    "    # unnormalized theta\n",
    "    theta = q.T.dot(w) + alpha\n",
    "    theta = (theta.T / theta.sum(axis=1)).T\n",
    "    pi = q.sum(axis=0).T\n",
    "    pi/=sum(pi)\n",
    "    return theta,pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def score(w,q,theta,pi):\n",
    "    return (q * (w.dot(np.log(theta.T)) + np.log(pi).T - np.log(q+.00000000001))).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def scoreSimple(w,q,theta,pi):\n",
    "    eq_log_p_w = (q * w.dot(np.log(theta.T))).sum(axis=1)\n",
    "    eq_log_p_z = (q * np.log(pi).T).sum(axis=1)\n",
    "    eq_log_q_z = (q * np.log(q)).sum(axis=1)\n",
    "    return (eq_log_p_w + eq_log_p_z - eq_log_q_z).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-242017.066781 -242017.066781\n"
     ]
    }
   ],
   "source": [
    "theta,pi = initialize(K,V)\n",
    "q = eStep(w,theta,pi)\n",
    "print score(w,q,theta,pi), scoreSimple(w,q,theta,pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def em(w,K,num_its=100,verbose=True,alpha=0.0001):\n",
    "    N,V = w.shape\n",
    "    theta,pi = initialize(K,V)\n",
    "    scores = []\n",
    "    for it in range(num_its):\n",
    "        q = eStep(w,theta,pi)\n",
    "        if any(np.isnan(q)):\n",
    "            raise ValueError('why now?')\n",
    "        if verbose: print it, score(w,q,theta,pi),\n",
    "        theta,pi = mStep(w,q,alpha)\n",
    "        scores.append(score(w,q,theta,pi))\n",
    "        if verbose: print scores[-1]\n",
    "    return scores,theta,pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 -232171.835312 -143644.609126\n",
      "1 -139143.810096 -137876.259714\n",
      "2 -137869.368754 -137862.528964\n",
      "3 -137862.516564 -137862.515553\n",
      "4 -137862.515553 -137862.515553\n",
      "5 -137862.515553 -137862.515553\n",
      "6 -137862.515553 -137862.515553\n",
      "7 -137862.515553 -137862.515553\n",
      "8 -137862.515553 -137862.515553\n",
      "9 -137862.515553 -137862.515553\n"
     ]
    }
   ],
   "source": [
    "em(w,5,num_its=10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAADQCAYAAACa9N1LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt8FOXV+L9PEq4CuQGCQA0BBKXcCeA9ykUqtLXcbK1U\nrSLqq69tFRHra+X1Aoi2tb79FVBRvBW51QpWuRqgFIVwURRESQBFBSQhCaAGsnt+f8xsstnsbmaT\nTXaSnO/nM5/deWaeM2dmJyfP7ZxjRARFUZSGQFysFVAURakt1OApitJgUIOnKEqDQQ2eoigNBjV4\niqI0GNTgxQhjjDfI5jHGtDLGJNn7c4LUW2yM8daQTv2NMdkOzhtnjFkUToYxpp8xZl8wuTWlf20S\nwbMKeu9O6tvnBH2GVdDX8XXrM2rwYkt/IMlvSxaRIr/j44PUGQbU1FqiHGCqg/PCXT+YDKdy6xJO\n7ynUeZE+k2g9w/r4WzhGDV5sKRCRIv8t4Hi2Maafb8cYMwxYDZhAQcaY1caYSX77s3wtRGPMrcaY\nfLvVmG2M6WyXp9v1Ztn/9dOBWX4ygtbzO74oyLEu/jICy4wxq+3PfGPMHGPMlACdZwZcI6QOdksz\nxycrxD35n+e1dU70kzHXrp8foEvQ8hD3lG6M2WaMmWKfv8/vdwt273mRPusAOeOC9A6uDCWnkusG\nfTaV3JOT5+NK1ODFlgqGK4DFwGS//XHA6yHOXQQM99sfCiwyxiQBc4ArgGQgN0DmUKCVLbtUn0rq\nGfv8VVgt02xb10oRkeH2ZwqW8fbXeSyw0IkOxph0YJ5dZwAwzBgzJvCe7PMWAZNsGeBnOOxz0+zP\nmcYaUghaXsmt9bNuS1KANVQ0+v73noqzZx20JS0iS0QkTkTigAlAjoisCyUnzHVDPptw91TF5+MK\n1ODFFt9/Vt/2WcDxxVgvtI8JIrI0hKzFWEbI9weULiLrRKQAq6v8gYgUAsexjFQpInK7iBwIKAtW\nL9HvlG0i8pzdMr0N6B/pS2/fyzBb53QgRUR2OtRhHDBXRHaKyH6s7n9ukHvynbfOljGV8s8UoIuI\n7LCvX0SZoQksD0eBiDxpf58LpETwHCr9jYJh/87zsP9pOPjNAqns2YS6p6o8H1eQEGsFGjjD8Psj\nBfL9D4pIod0t6QekAltDCRKRAmPMdmPMUKzuj39L8AG73Cff/5r+3wMJrJfjd2x/wLm5WN2lSFlj\nX2MA1h+VUx3SgW2+k+w/PJ/h9L+nFMobwv22oUBElhhjkoHFxpgUYAYwW0SW2vvlyiu5D//frrKW\nezDC/UahWAvcF/DPKtxvFkjIZ2MT9J6q+HxcgRq82JIb2LIKwmLg51j/qYMZBH9eB0ZgdUXug3Ld\njytFpMgYcyvWZElYHNQLNG4+Q9O1MtkBrMZqVQwAbolAhwIsw+47tz/QGdgeID8v4Lwku67POK4R\nkWftMbPVxphcW0aF8jCt62pRld/IGDMLOCYiz1Uip18oGYR5NpVcO+hzq6nnE020Sxtbko21BKV0\nC3LOYuBWrNbgmkrkLbHPHeDXNUwG8u0/gCT7uJPuVmX1+htjxtp6zwVWV7FbswTL4HX278460GEu\ncKuxlr+kYz2nZCq2rpbY5w21ZTxLWet3HFYrJdGvnoQprylC3WfQlqKxJq8mUbFrHkxOapjrLiX0\nswlHbT+fqKEGL7Zsw+o2lG6+2TbsF8gel8nBGpgu8j8WiD2WlYc1EO0rexasWVGs1tRUrAH+ctfx\nF+OgnmCNHU229U6j/BIaCfgM/L7EGOMJpbMT3e16U7G6dfuAVX6tHfGT4Rvfm2vr6rXrISJP2GXH\nsSdeRGSZiMwOVh6oX5j7c3LvQiXPGhgcQs5krBb/cb/x31vCPK+hIa6bG+rZhLunUM8tyPNxHSaW\n4aGMMbNEZKrf/jishzjeHggvPccYM8nvBx2L1fROj7RMcR/GWj5yS5AWnqJElZi18OzxhbF++0OB\noSKyFkg3xvS1D02yZy9z7PP6A9jnYXdpHJXVyo0pEWGMGU7w7qyiRJ2YGTwRmUf5GaK1InK7veu/\nPGGSiHQTkXX2/gSsViB2/WF2WYGDMsVF2C361wnuUaIoUSfsLK0xJoeyfrz/gHWOiGREWxl7EPRW\nrGluHyl266+/PbaSRPnp8tQIyhQXISJLsCYVFKVWCNvCE5EuItIVa5p+qL3iehgV12BFBREptI3a\nZHu6GxF51u6WptqGD6q2zklRlAaO03V4Sb6FnSKyPcTyiXIYP79OP/JDrdWxx9zEvs52LLegAr86\nPj/AAspam0l2OZWUJfuV+V+zTkylK4oSHhFx1AhyOoZnjDG3GGPSQhiyYAo8G2QLtzBxKOWNVg7W\n2Jtv7ZnP0+B1yha9pmNNv1dW1tkuC6anK7Y//OEPMdfBTXq4SRe36OEmXdyih0hkbRanBm880A1r\n7VUXojDIbA9YDzTG+FbXz8OanZ0EHBdrPdRarHVEY7FWle+UMheioVi+fo7Lqquzoih1G0ddWrEW\nv0Y1hpYEDFiL5bxcYa2cBGkVStmaurWRlimK0nBx1MIzxkwyVjws3xYY1UOpJpmZmbFWAXCPHuAe\nXdyiB7hHF7foESmOPC3slfBD7VZYvcEYI5GOASiK4i6MMUiUJy2y65uxUxSl4eF0WUqq3corjdYh\nIvfXjEqKoig1g1ODVyF7lqIoSl3D6RjeNgLC1EgNuJbVNjqGpyh1n0jG8JwuSxngJ3wY6oivKEod\nJOJoKSKyBgchwpXaJxqtVbfIiJac+iQjWnLqk4xIcdTCM+XzTjb4qCM33XUTuUW5GFPWihYR0lul\n88IzL9SqnJvuuons/Tv44st8Ej7Pp+QHKXTqkMLAzv3qnAw36eIWGW7SxS0yAuVEgtMW3n4sv9b9\nWP6sDTp+2ajho8iOz2Z95/WlW3ZcNqNHjK51Oe1SOvNxy92c6fYFPz15ijPdvuDjFrtpnxqYx9n9\nMtyki1tkuEkXt8gIlBMJMQ3xHmuqOmkhIlw44ULe7/m+FahKoPGaxrQf1760tZZghBbxHlokeEs/\nWyZ4aBHvpUWCh5YJXs6K95Bc/D3NWsEZLzSJA06Dp3E8TiJgCeDxeGjiEZqcgZRiKGgCZwHmTIT3\n1AhOAUkxluEmXdwiw026uEWGv5yrxzqPluK0SzuUsszygpVzYl2YKvUaYwzTbhzK0+9v5bxuXnrE\nx9H9onYkJrbC4zmOt+Q4IqeJj08iLiGZuPhk4hOSy3+PTyE+IYlXF65i53eLGZ/hZcF7cfRKuIFJ\nk24sdz0ROHQI9nwCe/bAJ59ATg60bwdFRfNJ/2ofd/bYRPwvwLMA7v3kR3Tr0ZS7/nu8T+FQdwLA\n008vYt8n3/GHHm8T/3NLxj2fXE3XHk25++7AxFjBiYYMN+niFhlu0sUtMgLlRITD8CvZQKL9PR3L\n8yLmYWGiEFZGqorH45GH/hwnK1YjDz/VTo4eXSGFhVvk229z5MyZQvF6vZXKWLFivaSnTxM6DpYb\nXkDoOFjS06fJ/PnrZelSkfvvF7nySpHERJG0NJHx40VmzxbJyhIpKrJkjLpyitxGE/GC5N6AeEEG\nM1hGjPi943sZPvwBGczgmMtwky5ukeEmXdwiI1CO/Xfs6G/e6cLjfLFdy0Qk104D16DZfXQLA3sI\na7MbcXHP/6NNm1ERy/jLX1aRm/s4NFoCe38JR6aQe2Yskyf/DyNGXEZGBtxzD2RkQJs2wWWMKd7C\n2Zwu7QAb4D6znT0Drnasx+UDmnD+mu34wqHGSoabdHGLDDfp4hYZweQ4xanBKzTG3AvswFqSUml2\n8vrOtk8egPgenNuyCWNGj6mSjOJi+/GfGcvOtfPhjCXnwgvjWbHCgYB9+/g0+z12dj+f546cJHX3\nUfKS2tKhYzItvnYehf/E4f2s63kBLx86HlMZbtLFLTLcpItbZATKoeBzx/UcT1oYY+7D6s7miJV3\nos5T1UmLkpIi3slKoVXaS/RKbU9y8hVVuv5FFz3I5s2P2nuCb0ztqqv+h3feeSR8ZRG4+mq48kqY\nMsUuknJLXKqCW2S4SRe3yHCTLm6RATXgaQFWtnFjTGexMrk3aD747H/ZWdiIaenXEh8XXyUZR4/C\nvn0jOPvs33PkyGP4jF2XLg9w110jKxfwj3/A55/Db35TWhStl8cNMqIlpz7JiJac+iQjUpzO0o4F\nZgG5xph04D4RWVbdixtjZolIhUjKxpgpvlakfe0CIF3sCMbVKasuHs8p8o7M4WSza6ps7IqLYcwY\nmDz5MoYMgWee+R++/z6epk093HXXSEaNuiy8gJMnLUP38svQqFGVdFCUhojTFt40sdI1AqUBQatl\n8IwxtwJjCQgdb/vqDgdm25nMEJG1xph0Y0w/7KZQVcrEznNRHb76ag57TjTiqh/eWqX6InDHHdC2\nLUyfDnFxl1Vu4AJ55BG4/HJrUxTFMY5naSvZjxgRmWcn8qlwyO/7tcBK+3suVtCCVMoykEVaVi2D\n5/F8x4GDT/DKQcPmUZdWScbTT0N2NmzaBHERezIDu3fD/Pmwa1eVrq8oDZmwBs/uEgIUGGNWYgUA\nHU4NzdLarbC1xhhfqy+R8sY1FSuFY1XLqsXXXz/HMU9r+p97CQlxjoc/S1m5EmbNgs2boUWLKigg\nAv/1X/CHP0C7dlUQoCgNm8r+artgtbi2UubrtIbyrbBokhKkrEZHNh9++OHS75mZmSGTk3i9xXzx\nxRO88nkit10UrGEanr17YeJEWLoU0tKqpiuvvQaFhXD77VUUoCh1n6ysLLKysqpUN6zBE5EnqiQV\nK9NZkOJ8CZGM29e6CyguoHxy7jz7e6RlyX5l5fA3eOE4fPhF4pp0Zf3Xu1iYlumojo/jx+EnP4HH\nH4dLq9YTtgzdlCmwbBnEV22yRFHqA4ENk+nTpzuuG3m/zCFVmBVNt2eAU4EUe+LhdWAgVl7ZdKwx\nOVOFss6UjedFjNd7hoMHZ7BLruGn3dNpFO98ZrSkBH7+cxg5Em65pfLzQ/LQQzB6NAwZUg0hitKw\ncZqXNi3aF7YnLAYaY24BK+G23foTrLE78c2q2sELCkRkZ3XKqqrrkSMv06xZV176JJtxF0TWnb33\nXmvo7amnqnp1YMcOWLgQZsyohhBFUZzmtNjnvyylvuDE08LrLWHLlh6kdprJ4Fdu5fC9h2kc39iR\n/OeftyYp3n8fkpOrqKTXCxdfDDffXM0moqLUT2rC02JtQJpGEZFpVdKujnH06EKaNDmHlV8e5sfd\nf+zY2P373zBtGmzcWA1jB/CCHQX217+uhhBFUcC5wVtN+TGwBhE1VMTD558/Rteuf+G3bz7GPRfe\n46jewYMwfrzlCNG9ezUUyMuDBx6w1rNUadGeoij+OM1atqSmFXEj33yzlISEJE436sUHRz5geJfh\nldY5edKakb3vPrjqqmoqMG0aXHst9O1bTUGKooBzX9p+WBGPjwOLsCKmVNuX1s2IeDl48FHS02ey\neO8bXN3tapomNA1bx+uFX/0KBgwo59NfNd57D1assEIcK4oSFZz2k57FWuKx33bqf6DmVHIHx479\nE2Mak5LyI5bsXsK48yufnX34YThyBP72tzBR1Z3g8VgOt7NnQ2JiNQQpiuJPJOGhCvzCudTriMci\nwsGDj3DuuQ9x7NtjZH+Vzciu4UM2LVoECxbAli3QpEk1Ffjb3yxDd9111RSkKIo/Tg3eNmPMHCDJ\nGDOTeh7xOD//X4iU0Lr1T3h+x3yu6noVzRo1K3fOW29t4C9/WUVxcQLFxSV8/PEINmy4jLPPdn6d\nCwZeAK3h8IE82n5VxNEOrTi/17l02/ABL67b4aiZGFRG/650a9ONF//8YtX1iIEMN+niFhlu0sUt\nMgLlRILTSYvJfq5iOSJyf0RXqUOICAcOPMK55z6IMXEs2b2Em/vdXO6ct97awJ13vsEXB46Rxn84\nwMWktHmDL7+Evn2dh3rKL84nvugIN38Jj56EBw99z2vJR5Eu7eGCC6onw+t8It0tMtyki1tkuEkX\nt8gIlPNkBPWcTlrsw1qDtziIv2u94vjxNXg8RbRpM5b87/LZfGgzSyaUn6SePfsfnD7wLr/lCI9y\nmAc5xWvfJDF7tvCjH11iZUjyevGWePB6vHjOlOApsfZ9n16Pl6ScPF7/HvrYv/Xsk3D9e3Bdizz+\ns2Y9IoDXi/Uh9ksh4BU8Xi8AKbl5/P27ijJ++VE2qxa9BSL41lYLVh2viJ3ryfreOjePV4PIuH5X\nNssXLCv/Mvot1Pba5QJhZSyd+3r5hxxksbevKJScibuyWfzXv5fTRfD/XkYoGb/alc3CP79cmSoA\ntMnN45UQerz6x5dCy5DKZfxqVzavPLHAkR6VyXlp1ouhKzrUZcGM+SHr+dM29xgvh5Lx2PyQ9fyL\nw8l44ZH54SuHkBOJwYskp8VQrNBQQ4FcEbk2guu4kkBPCxFh587LOOec2zj77F/y4s4XWf7pcpZO\nKB/voHvj9iw6c4w+lJSWfQw8BzyFIc7+lbwYvMQhIT7ncpJfAB38ZH8JvGrg5iYGsQPF+K8h92kr\ndgLw+cVefuWtKOPlOLihaXyF+v74ZL34nYcbg8hYEAe/alYxUEEweQu+9XBTEBkvxsHE5hX/r4Z6\n617+tiSonBfi4PqznI3AvHyqhF+HkPFLhzJeCSvDmS/1K6fOBJUxPwIZPjk3V1NOKBnPx8EvWziU\ncfIMt7hARqAcg/NE3E7zt6YBU4BVwDZgjtM8kG7eCMhLm5//rrz3XlfxeM6IiMioV0fJax++JoF0\natRDDln/lEu3QyBpjXvK6e89UnLGKx6PSGBq2pKTJXLo9UMy4+YZ0v227sIA5IEW5eVMa4Ek9mxf\n4ZqhaHxuYr2R4SZd3CLDTbq4RUagHGogL20usB2YKvW4S3vw4CP84AcPEBeXQOH3hWz8fCOvjX2t\nwnlJPfrwwq6veJCi0rK/0paW3X9IoyblV/p4vvOQ/04+ny3+jAXfLOCNjDfo1rUbM4bP4I5Nd/D6\nuYfZdLLs/EMtoGefzo51TmnVjNdbFNYLGW7SxS0y3KSLW2RUkHPQeT2nBi8FK0T6eHuWtl50af0p\nLPwP33+fy9lnXw/A8k+Xk5mWSasmrSqce0WPluTvOsVP6EURKcAZClp8yfr1fwPAW+wlf1U+R18/\nys5NO/nniH+yssdKfvyTH7MycyV921meE79v/Htoncx7ez+niSRQnCAM6PNDurXp5ljv5MbJ0DqZ\nD4+ewOuNJy7Ow/nn/6BOynCTLm6R4SZd3CIjUA4HDzmu5zRaSiIwARiBFW9uq4jcFpGGLsR/DO/D\nD39E69bXcM45kwG4ZuE1jD1/LBP7TKxQ7620IXgOtuQaDD0v6E7HHzThlRenwXY4+vpRjr15jN2X\n7mbJ4CV8EP8BkzMmc0fGHbRv2b6CrHfuvhvzzDNcJcI7zZtjXnqJq8aOrXCeoijBiSRailODlw0s\nBJaIyIHqqecefAavqGgrH388hsGD9xEX14QTxSfo+KeOHPzNQZKaJpWv9PXXHO/YiYHt+pLfaD/J\ncYmcjD9JS29Lusd1Z9ikYSxotoAzcWf4zZDfMLH3xApr+HxIbi6/O+88/ujxYM9B8LvBg/nj5s0x\nydmpKHWRmggPdT+WL+0DxhgBxovIuqoq6DYOHnyUTp2mEhdnuUj867N/cXGniysaO4CHHmJBYjK5\nt24DoMB2OvmGb9i/bj8l7UuYdeEsRnQZQZwJ47lXUsLKq69mZFwcxuMBrNmmq3btYtWyZdrKU5Qa\nwKnBmwmkiUihHYZ9EZZvbbVwmIh7lohMNcZMkhpIxH3ixE5OnNjKBRcsLC1bsmdJ8MjGH36IvPkm\n01ukA8cqHO7foT+rJq5ydvOPP05WURFNLryQzX6tORGheMUKNXiKUgM4zksrIoUAIpJrjKm2L62T\nRNx20SRjzBhgsn28ysm5JUgibqt1dy/x8Va389TpU6zKWcWcUXPKnygC997LnrH/Q9HbLweKAaBF\nY4e5Fzdvhr/+lZnbt0OHDpWfryhKVHAaLaXQGHOvMWaoMWYKUfClFZF5WMtdKhwK2J8kIt38utAT\nsMJUQVmC7Ql+OoUrq0Bh4cbSiQqAd/a9w+AOg0ltHpDG9p134PPP+cNXkxGJzH+vHEVF8Mtfwty5\nauwUpZZxZPBEZDxWi2m8vT+hJpQJkaoxxc/QQpQTcXfs+Fvi488q3Q/anS0pgXvuoejBJ3h7TSMS\nGn8XwV0FcOedMHw4XHNN1WUoilIlwnZp7TV3/hQAqcaYGVIzOS0qJOL2G48bbru3QRSTczdvfh55\nee+QmjqS7858x9ufvc3TI58uf9Jzz0G7drxw7Me0bw+fFxXRa20vdsftZlCnQaV5LtKS0sJf7O9/\nt+JHbdsWLfUVRYmAysbwsqli/opoJOK2x/ny7Dp5WGsAq5KcO2Qi7iefXEbTpl0w5j0S0hPo374/\nbc9qW3ZCURE8/DDyr7d5/kbDsWN5NJNGPN7xcW7reRv/+f1/Qj8Efw4cgLvvtrrGZ51V6emKogQn\nKyuLrKysKtUNa/CkGrksQs2KhiFYIu4cLKOLXb7a3o9aIu5Zs14p/T7xHxMZ1zGgOztjBowcyTZv\nP/LzoWnTtWRkXMLWQ1sZNGKQszsrKYGJE2HKFOjf31kdRVGCkpmZSWZmZun+9OnTHdeNWSosh4m4\n1wLD7OUlx6QGE3EXlxSz4tMV/KzHz8oKDx6EefPgsceYPx+6dYN27TYw5IIh7Gm7hyHdhji72Rkz\noHFjuMdZ1jNFUWoGx+Gh6iP+rmVvffoWszbNYsNNG8pOuO466NaN7+6fTseO1qTqqVO9+ct1M3jw\n9IP86fY/kZmWGf4imzdbExS6BEVRaoSoeVoYY3IoG8Pzn1DIEZGMKurnSirMzm7ZAuvXw7x5/OMf\n0Ls37NiRj8ezn/aH2rM3bS8D2g8IL7SoCK6/HubMUWOnKC4gbJdWRLqISFes0FBDRcQXNWV/bShX\nW5z2nObNvW8y5vwxVoEI/O538Mgj0KIFzz8P550HvXtvYsiQIezavYu0lmm0bNIyvOC77oKhQ+Fn\nPwt/nqIotYJTT4sk35iYiGw3xgRxMq27vLv/Xbqndqdjq45WwdKlcOIE3HAD+/fDhx9Cs2bQqtUG\nBva8mF2bdzE4fXB4oQsXWrllt2+v+RtQFMURTictjDHmFmNM5xDLTeo0S3b7dWeLi2HqVHjqKYiP\n58UXYcIE2LABDh/eSN/mfdnXex+DO4YxeAcPwn//N7z2mi5BURQX4dTgjQe6AnOwln2MrzGNapkS\nbwlv7H2Dsefbzvp//Sucfz4MG4bHAy+8YO326nWSPXt20fVYV/a03cOgDiGWpHg81hKUe++FAZWM\n8SmKUqs47dIOwHbeBzKw1rTVi/BQ6w+sp3NSZ85NOhfy8qwlJBusmdq1a6FNG/j4Y+jZ8z2M6Ufe\nljwOjT5Er7a9ggucMQMaNbIMnqIorsJpC28WVnioZKzFvE/UnEq1S7nu7P/+r9V/Pf98AObPh1//\nGlasgISEDVwy5BI+KPyA3u160yg+SIal99+HZ56Bl16CuJgtcVQUJQRO/yrLhYeivFN+neYfn/zD\n6s5++im8+io8/DBgNfbeeQd69YImTWDPno30T+5P7sBcBnUM0p09ccKKgvK3v+kSFEVxKTELD+UW\nzml5Dl1SusB991muX23aANZ8w9VXW73bkSOL2bp1K90Lu/Np+qflxu9KF27fdRdccQWMGROL21AU\nxQGRhIeKo4bDQ8WC428f5+6r+vDNurWWc7/N/Plw881Wd7Zbt2y6d++OZAu7mu0qNXgiwu9uuQVZ\nuNDyqPjzn2N1G4qiOMDppAUiUm/G7fw52OcAN/7VkDv5bto0bQrAjh1w/DhccAF88gmcPLmRSy66\nhAOLDnBqxCm6pnQFYOXSpbBoEasWL+aqdet0CYqiuBxHLTy7K5tvb3nGmCtrWrHa4voPoZE0Z9Cs\np0rLnn8ebrrJGsMbPhw2bdpARocMcvvlktExw+e7x8onn+SPJ0/yTlISoktQFMX1RDJL29l2Lcug\nHs3S3vg2fDN1GsaeVf3+e8tJ4sYbYflyGDXKw6ZNm+j5fU9yeucw6ByrO7ty6VJG7txpZRrLy2PV\nsmXV1kVEWLJiBTf97ncqw2W6uEWGm3Rxi4xIaPCztFvOGHYePU5JSQkAb7xhhaxr185ah9ep0wd0\n6NCBhB0J7Gm7h4wOGaWtuxHFxQCM+PZb/jJ9epV/NN+PfuHNN3PD3r3k2ro0VBlu0sUtMtyki1tk\n+OREgtMxvEJjzL3ADqA/9WiW9n6vsPvpp5mam8tTy5bx/PPWZEVWFvzwh/DRRxu55JJLKFhWwM4L\nd5JxTobVutu1C4AlTZrwZHo6O8eN47xPPiknu8Tr5ZTXyymPp2zz2z9ZUsLGNWtYuXIlX/fty5mJ\nE8EYDn32GX/84gtH+osIH777LlmrV9d5GW7SxS0y3KSLW2QEyokERwZPRMYbY+7DmqXNqU+ztACJ\nJSV0Pv98DhywJiyuucZaoTJ6NGzYsIFRg0dxuMNhWjRpQfuW7fnzihXkpqZywwUXcHTsWGTwYDCG\n3Z9+SqfNm0sNWokIzePjOSsujrPi42kRH89Z9tbcGD564gm+OH2aM7ffDn65aYu9Xg7ZrcdwiAir\nHnmEfd9/z+k6LsNNurhFhpt0cYuMoHJefdVRPdBZWgCWtWvHmDvu4Lnn4Be/sBYaL18Oy5cLf/rT\nRu7PuJ+si7MY1GEQIsI3LVqwsndvTtxzT7kfrV9yMkv69Ss1cE3j4jAmdFxCeflllr71Fk++9BK7\nevbk2wEDwBi6NmvGH7t2daS7LFhQb2S4SRe3yHCTLm6RESjnfce1nM/STjLG7PPbPovgGuHkzgq2\n7x+RxRgz1p4ljkpZIA9ffjlHLr+c9u078MILlivZxx9bxxIS9tK8eXNafNyCvel7GdRhEMYYnh03\njh6dOtF23jyaZ2db8fOAs+LiOLdpU1o3bkyz+Piwxs7Wj3GjR7N5/nwW9OjB4Jdeonl2dkTjEvVJ\nhpt0cYuM27dnAAASoElEQVQMN+niFhmBciJCRCrdsBLnJDo51+kG3ArsCyjLBz4DrrT3+wNj7e+T\ngH7VKQuig/hYvVqkb1/r+4wZInfeKTJ37lyZOHGibE7bLIOfGSzv7n9XRESmP/WUXLR8uXxfUiKL\nly+XwTfeKM1nz5bL7rpLqoPX65XFy5fLjb/9rcpwmS5ukeEmXdwiw/47dmR3HOW0MMYskhoYtzPG\nrBKREX77Y8UvjaOdF3eViKyzk/H0x85eJiJrIy0TkdkB1xff/V93HVx0kZUn++KL4aGH4OWXr+fi\nPhdzwZMXcPXdV3P4nsOsP1rAbe+/z9YLL6S97TMrIix96y3eWreOF/74x2g/JkVRwhBJTouwXVpj\nzEzb6CQZY7J9+0ESdEeLFD9/XbByzPovgUmtZllQjh+Hf/3LMnrHjsFHH8Hll1sTFn3i+3B06FHS\nktL4yhPPr/fuZcm6daXGDsqa12rsFMXdVDZpsTrgs0YRO5etMWa43SoDK79sjfLaazByJKSkWJGd\nhg6FI0cOUlxcTOq+VNb3Xk+/9gO55qOPeOyf/2TIqFE1rZKiKDVAZQavv4jMDtGiWxuuYoiJgnz/\nLmuQ833H87AiKxdQli0tyS6nCmXJfmXlePjhh5k7F4YNg6ysTFasyCxdjnLZZZdR9O8idv/3bnbF\nncvlXi+Tli2DP/0p3K0rilKDZGVlkZWVVaW6YcfwfGNqdiLscoQyXBFd3G8Mz27RZYtIoW1gF2K1\n7gaKyLN2N3d1dcokIBm3MUZ27BB++lPIzbWis7dtC3v3woMPTqJn154MeHQA4x65hXa9H2Lr2ztp\n3KgRPP54dW9dUZQoEbW8tECGMSZU/tlqGTxjzDhgoDHmFhF5zp5cGGsv5TjmM07GmIG2MSyIRlkg\n8+dbfrPx8ZZ3RffucPbZsHHjRq7vfT3HB3r45sTnbBnwIxr/4nfw739X57YVRYkhlbXwxlGWiLsc\n0WjhxRpjjKSmCtnZkJYGv/kNtG4NkyYdoUePHrx90yZmfreGnL6vsKvtNCve3fr1sVZbURQ/otbC\nE5El0VHJvYhMITHxAUSSWb7cSkm7ceNGBl90ER+tPUazyV9xxQ+GwF+ft1YlK4pSZ3HqadHP9rDY\naoyZYoypN3HM8/MfYMSI37Nly3FOn4Y+fWD9+vV82eU80j4TSpJzGNS8G2zaBOPGxVpdRVGqgdPw\nUM9iZSvbby/efaDmVKptksnOfowbbnic0aMt19gl69bRvsUPSex5FluPbGXQpgMwfrxGNFaUOk4k\nwQMK/HxD6008PItkvv76LEaPhoU5ORzZv59/JVxE8WWFFBYX0vXFN+Gll2OtpKIo1cRpC2+bMWYO\nlsfFTOpRPDyL45w6dYqOF37L5GXL6J+RQXz2GXJ65ZDRrCtxjRrD4MGxVlJRlGriNGvZZGA7sJ96\nFw/vOJ07/57Mq+/j5/s+4qKDBxmdeQWF/ylkd+vdDNp/2pqsqCTyiaIo7sfppEVfEZlnG75FxphF\nNaxXrTFy5OP07fcox249zOWJiRRu305Gxwwat2vMtmPvM2hDDkycGGs1FUWJAk67tM8ZY/raHhfb\ngNdrUKda5c03Z/N26yIatT7DzA4d+OCDDzjvxHm0uqQVWz9/j4wul1nuF4qi1HmcTloMw4qJlwMM\nEDuhT31g9r+/4XS7pXT5++ds/+lP6dOnD6ffP82J4Sdovu8M7X91R6xVVBQlSlQWHmqS7dQ/HsuV\nrAswPlwE4brGg//3X8Q1ha9F2LhxI5deeimF/y5kb9P1DPo6zgqjoihKvcBplxZgH1Z+2no1ei93\n3k7JhQMwxrBhwwYGnzcYBHbsWcygTkMgwfHKHUVRXE5lBi/JjlHXJchWP7BnX70eD1u2bKHn6Z4k\nXtSSrXkfMuiqm2KsnKIo0aSy5kuu/bm1phWJGSJgDCePHqVr166wHZon72dnWw8DBl0Ta+0URYki\nMQsP5RaaT3kcz/jhFHz5JaMvvZTC1YWc7vw8aZ3b0apJq1irpyhKFKnM4GUTIjxUfeHbbatoe+x6\nTpw+wEU330Xxgu/Z3XoNGedp605R6hsRhYcyxswUkftrVqXaxnD04EskJLSgl+mFp91nbL2wI4PO\nvTjWiimKEmUimaUFK89E1HCYiDuqybmD8xGNGp1Fk11NSMzLYsvZHgZ1GFTd21MUxWVEavCi5lJm\njLkVCMyVMckY8xnWAuegZcaY/gAistbe7+e0LLQ2G0hJ6UThO4eIT9xLTvHX9D67dxTuUlEUNxGR\nwYtmBGQRmUfZLLCPSSLSTUTWhSmbABy3v+dieYFMoCyCS7iyoJx11v/j2jFXc+qzM3z26870atuL\nxvGNq3xviqK4E6fBA7LtiMf+20pjTFqU9QlMxB2sLKqJuEeMeJBGjY5yQ5+f0EL2sW1gG+3OKko9\nxWkLLxuYLCJdgduAJcA0YF40lRGRZ+0uaKovEXewMqLo7fHMM7+iVasWtPznlySmn2JL0R4yzgm1\nEkdRlLqMU7+pdN9YmIisMcbcJyL3mzAx4qKRiNsYkx5YRtWSc4dMxH3vvfeSmJjIrFUv8qMJF7H1\nq608csUjIe9LUZTYUp1E3E4NXqEx5l5gLfZYmF9rKyi2S1ok5GK1JMHqfq7BMlj+Zavt/YG2Lun4\nJd0OU9bZLqtAUlISd1wzhp6PDib9kZ4c//tMuqV2i1B1RVFqi8zMTDIzM0v3p0+f7riu04jH44E2\nWN3YVKzoKfh9Rox/Im77GmuBYXbMvWMisiNI2U4R2WHXL02w7bQsmB4bN25k4M4CmqaWsLN4Fxkd\nMogzkU5eK4pSFwibiLu+Y4yRNm3akP3tKL7/ye28OuktSrwlPDb0sVirpiiKQyJJxO10lnaoMSbf\n3vKMMVdWT0X3cGmXLhQ1GUji6HPZ8tUWnaFVlHqM077bLKCziKQAGcATNadS7XJp0QkKz/Sg1SWt\n2PKlGjxFqc84NXj5vrDuIpJLPcpLO+jAd8QlN+OrFl/RLKEZ7Vu2j7VKiqLUEJHO0u4A+lOP8tK2\n+uEv2Nj8A5K/StbWnaLUcyKZpTXYs7L1KS/t73O3kvmTTO3OKkoDoLIkPjN9G9ZylAIsj4cZtaJd\nLXBd0XV0Gt5JDZ6iNAAafADQg6f3s373W+w8vJMB7QfEWh1FUWqQiAKA1kdasp/Xnl3KD675AYlN\nE2OtjqIoNUiDdykYRTGppz6lQ8nZsVZFUZQaJmwLzxiTQ1mXNsXvUI6I1IuQIkv6JfDvXim0zvku\n1qooilLDhG3hiUgXOyTUdmCovfB4GLC/NpSrDe7evABPv7Y8NvWZWKuiKEoN47RLm+RzxheR7Vgh\nmOoFN1x5A58c/YRebXvFWhVFUWoYpwbPGGNuMcakVZ4Qp26RcSiDFl+1YMptU2jIgRQUpSHgKFqK\nMSYJKzRUP6zu7Qyfq1ldxhgjGddmcPbFZ/PmnW8SLqCpoijuJJJoKY5cy0SkAJhaLa1cSuexnRl9\n3mg1dorSAHAaHqqfnbhnqzFmijFmTE0rVlts+WoLGR3qxYSzoiiV4HQM71mscOn7RWQ28EDNqVS7\n5H+Xz3mp58VaDUVRagHHC4/tbq2PqISHMsbMCtjvb4wZ6z8xYu8PjVZZIBnnaEh3RWkoOP1L32aM\nmQMk2YEEqh0eyhhzKzA2oPh+O0NZkt2N7g+l+S6oblkwPTRggKI0HJyGh5qMNTu7H8vLotrhoURk\nHlamMqA0qc9W+9hse93fBOC4fUou1qLnCZQZ3EjLKqAGT1EaDk4nLcaIyDwRmSwiz9bQpMVArNBT\n/YwxU+yyJMp3n1OrWVaBWY/OIvPGTG78zY3V015RFNdTmS/tWOBarFSJP/c71B9YVgP6HBORHcYY\nX2pGsAKP1hjvdXvP+lJvnOUURQlFZevw1mCFdZ8KzMQyPkJZNzMkISYK8u0xumDkUWZ2CrCSBRVQ\nFrQgyT6HKpQl+5WV513r40DBAbKyssol+FUUxX1kZWWRlZVVpbqVxcMrxMpnMRMYZndnZwJzqGTi\nQkSejVCXJcA4+3sSsAXLAA4E1gLpwGosoxtpWWe7rCJXWB9p+9PU2ClKHSAzM7Pc3+r06dMd13U6\nS7sYK/oxWK2+xY6vEAJ7kmKgMeYWABHZDxTYXdkUEVnmC1hgjBkKFIjIzuqUVVdnRVHqNk59aVeJ\nyIhQ+3UVY4zwsPX98v2Xk/ViVizVURSlCkTdl5ayNI1rgOHUozSNl++/HIC0pLTYKqIoSo3jqIUH\nYIy5D2tMLFdE7q9RrWoJY4xoSChFqdtE0sJzbPDqI2rwFKXuE4nBa/DRUhRFaTg0+GgpiqI0HGIa\nLUVRFKU2iVm0FEVRlNom0mgpuUQpWoqiKEptE9EsrTGms+0RUS/QWVpFqfvUxCztWGPMPmCuPVur\ns7SKotQ5nLqWZYvIwFD7dRVt4SlK3SfqLTwqzsrqLK2iKHUOp760+40xKynzpfXFuxMRea6mlFMU\nRYkmTg1ejr2BZfQEK2adoihKncGpwdvmy/4F1iRGmMjFiqIorsTpGN5txpi/GWNaGWMWYXdrFUVR\n6hJOFx6PB4qwPCy2iMht0bi4w0Tcs+zPGknErShKw8HpOrw5WHkhugJX2cFAq4WDRNx97bJJxpjP\nsMcQo52I2y1UNSlJtHGLHuAeXdyiB7hHF7foESmOfWlFZIKI5IrIcKKQOtFBIm5fDopJItJNRNbZ\n+9VJzu1a3PICuUUPcI8ubtED3KOLW/SIFKcGL9s/Hh5+hiqKBEvEDZBid0ujkZxbUZQGjFODN4/y\n8fCm1ZA+x/yyjY0FK92j3S1NtTOQQQ0n51YUpX4SkWuZMWaRiExwkrXMSSJufzm+lqOILLXrdsEa\nt8u3y6ZgdVG7AKtFZK1tFNOxWm/hysYBnW1j7a+j+pUpSj0g2lnLIo6HF6VE3IWU5cNNxUqmnU2U\nEnE7fUiKotQPYhYPz2Ei7rXAMLvsmCbiVhSlOjTorGVKcIwxUwK7/0rsMMbMEpGpfvtjsXpZ6VXo\nSUVTD9+wVZfaTt0aqItfedh313FOi/qGWxYlG2Mm2dvMWOrhwxgzjBh70gRbgB5DXWL6ngSuV43V\n+tIgegwF1tgGN91vQrHWdfErr/TdbZAGzy2LkmP50oTBDU1+/wXoMVswbl87135PcmOhS+B6VYKv\nQ42FHul+186192uFILqUHqqsboM0eMTopQlCzF6aYBhj+vkHiYiRDoEL0HfEUh/A5/6Y7gJdwCXr\nS+3lYr7udH/s3yxWOH13G6rB05cmOCkxvj6EXoBe69gGbr8xJh93Bb11zeoCu7e0zQWTgo7e3YZq\n8EBfmkAdYt6686PCAvRYYIxJAvYBk4BnjTGdY6WLHwWU/XEnA3kx1AVgqIjUlCOCIyJ5dxuqwdOX\npiLp9gD9rVjufLEaO8sDfJnxCoCMGOkBlqGba48njqdsnWgseZ2yoY+g60trC2PMrb4Z0RiPPzt+\ndxuqwdOXJgARWWr/YQuQSOwmL5ZQ9tv4FqDHDBEpsj/XEoME9EHWq8ZkfWmgHvaM6Ezbxz6fWnxf\ngjwTx+9ug12HZy8zyKWW1zIF6DAMWIQ1PpQCjPOLCtNgsX+bfGBgrFu+fsEyUmL1nijRo8EaPEVR\nGh4NtUurKEoDRA2eoigNBjV4iqI0GNTgKYrSYFCDpyhKg0ENnqIoDQY1eIqiNBjU4Cm1jjEmyRjj\nNcZk+20r7WPD7GP9Auoct9MM+JeNNcZMMcYkVtfn1l+GT2515CnuxGlOC0WJNrkiMjDUMWAycBuU\nBlfII8BlyJcQyhiTDlwLLKXqpPpk+CeaUuoXavAUtyFY+VMG+JVNxvKxTfI/0edTaZcPM8aMEZFl\nxpi5fvUnicgO+9zhwFDgMPCtfdwXHGAqZflTDJZb2/3GmMVY/pkpwFQ7C94wW6c8+/pz1e2sbqAG\nT4kV6caYbL/9rSJyO2Vhu9YYY4baTvsDsAzS+AAZvhbfLOzET3bEDLHTiiZhZbnrap83ACvg6gA7\n9Wd/YCZWy3AmkGyX+7q2U4D3ReRJY0wiVhQXX5SdziIy3i7fBqjBqwOowVNiRbguLVhBFSYbYwDW\nRCB3ANDZGLPI3j/ud2y1vT/cGBOY+yBYfMQuwGIAESm0dfGxJkS54mJ00kJxJXYYpAysruMcnAds\nzQa2i8gEO52oz/CJLWMaVrDV27C6yeHIwYpE7Rsn1EgbdRw1eEqsSA+Ypd1qdw+FMsOyEOgnIgcC\nyv0RrFBS/e0xPF9CpFV2lzkn4NzXgfF2C3AYVmuwr58M32yv2HEKM4wxq7AM53g/ORIgV6kDaHgo\nRVEaDNrCUxSlwaAGT1GUBoMaPEVRGgxq8BRFaTCowVMUpcGgBk9RlAaDGjxFURoM/x9NbFHPbKCq\ngwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8771a80950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(4.5,3))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "syms = ['o','v','^','<','>','1','3','8','s','p']\n",
    "for sym in syms:\n",
    "    ax.plot(em(w,3,num_its=15,verbose=False)[0],sym+'-')\n",
    "    \n",
    "ax.set_xlabel('EM iteration')\n",
    "ax.set_ylabel('expected log-likelihood lower bound')\n",
    "ax.set_title('EM variability across initializations')\n",
    "#fig.savefig('em-sensitivity.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Word Sense Induction #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "datadir = 'books'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "import os\n",
    "from collections import Counter, defaultdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from nltk.tokenize import sent_tokenize, word_tokenize\n",
    "from nltk.corpus import stopwords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getAllInstances(datadir,instance_getter):\n",
    "    instances = [] #list of Counters\n",
    "    for filename in glob(os.path.join(datadir,'*.txt')):\n",
    "        with open(filename) as fin:\n",
    "            instances.extend(instance_getter(fin))\n",
    "    return instances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getWordNeighbors(fin,word,window_size=2):\n",
    "    outputs = []\n",
    "    def getNeighborsAux(wordlist):\n",
    "        if word in wordlist:\n",
    "            idx = wordlist.index(word)\n",
    "            neighborhood = [wordlist[i] for i in \n",
    "                            range(max(0,idx-window_size),idx) \n",
    "                            + range(idx+1,min(idx+window_size+1,len(wordlist)))]\n",
    "            return [neighborhood] + getNeighborsAux(wordlist[idx+1:])\n",
    "        return []\n",
    "    \n",
    "    for line in fin:\n",
    "        for sent in sent_tokenize(line.decode('ascii','ignore')):\n",
    "            words = [word_item.lower() for word_item in word_tokenize(sent) if word_item.isalpha()]\n",
    "            outputs.extend(getNeighborsAux(words))\n",
    "    return outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getWholeDoc(fin,word):\n",
    "    output = []\n",
    "    whole_doc = []\n",
    "    for line in fin:\n",
    "        for sent in sent_tokenize(line.decode('ascii','ignore')):\n",
    "            whole_doc.extend([word_item.lower() for word_item in word_tokenize(sent) if word_item.isalpha()])\n",
    "    if word in whole_doc:\n",
    "        return [whole_doc]\n",
    "    return []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " The fictional town of Grantville, West Virginia (modeled on the real West Virginia town of Mannington) and its power plant are displaced in space-time, through a side effect of a mysterious alien civilization.\n",
      "A hemispherical section of land about three miles in radius measured from the town center is transported back in time and space from April 2000 to May 1631, from North America to central Germany.\n",
      "The town is thrust into the middle of the Thirty Years' War, in the German province of Thuringia in the Thuringer Wald, near the fictional German free city of Badenburg.\n",
      "This Assiti Shards effect occurs during a wedding reception, accounting for the presence of several people not native to the town, including a doctor and his daughter, a paramedic.\n",
      "Real Thuringian municipalities located close to Grantville are posited as Weimar, Jena, Saalfeld and the more remote Erfurt, Arnstadt, and Eisenach well to the south of Halle and Leipzig.\n",
      "Grantville, led by Mike Stearns, president of the local chapter of the United Mine Workers of America (UMWA), must cope with the town's space-time dislocation, the surrounding raging war, language barriers, and numerous social and political issues, including class conflict, witchcraft, feminism, the reformation and the counter-reformation, among many other factors.\n",
      "One complication is a compounding of the food shortage when the town is flooded by refugees from the war.\n",
      "The 1631 locals experience a culture shock when exposed to the mores of contemporary American society, including modern dress, sexual liberation, and boisterous American-style politics.\n",
      "Grantville struggles to survive while trying to maintain technology sundered from twenty-first century resources.\n",
      "Throughout 1631, Grantville manages to establish itself locally by forming the nascent New United States of Europe (NUS) with several local cities even as war rages around them.\n",
      "But once Count Tilly falls during the Battle of Breitenfeld outside of Leipzig, King Gustavus Adolphus rapidly moves the war theater to Franconia and Bavaria, just south of Grantville.\n",
      "This leads to the creation of the Confederated Principalities of Europe (CPoE) and some measure of security for Grantville's up-timer and down-timer populations.\n",
      "[[u'the', u'fictional', u'of', u'grantville'], [u'west', u'virginia', u'of', u'mannington'], [u'from', u'the', u'center', u'is'], [u'the', u'is', u'thrust'], [u'to', u'the', u'including', u'a'], [u'with', u'the', u'dislocation', u'the'], [u'when', u'the', u'is', u'flooded']]\n"
     ]
    }
   ],
   "source": [
    "with open('books/1632.txt') as fin:\n",
    "    for line in fin:\n",
    "        for sent in sent_tokenize(line):\n",
    "            print sent\n",
    "with open('books/1632.txt') as fin:\n",
    "    print getWordNeighbors(fin,'town',2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "window_size = 5\n",
    "word = 'plant'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "instances = getAllInstances(datadir,lambda fin : getWordNeighbors(fin,word,window_size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 494,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "instances = getAllInstances(datadir,lambda fin : getWholeDoc(fin,word))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def getDataFromInstances(instances,min_count=1):\n",
    "    stoplist = stopwords.words('english')\n",
    "    counts = sum([Counter(set(inst)) for inst in instances])\n",
    "    vocab = [word for word,count in counts.iteritems() if count > min_count and word not in stoplist]\n",
    "    V = len(vocab)\n",
    "    x = np.zeros((len(instances),V))\n",
    "    for i,inst in enumerate(instances):\n",
    "        for word in inst:\n",
    "            if word in vocab:\n",
    "                x[i,vocab.index(word)] += 1\n",
    "    print x.shape\n",
    "    return x,vocab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(54, 18)\n"
     ]
    }
   ],
   "source": [
    "x,vocab=getDataFromInstances(instances,min_count=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "54 (54, 18)\n"
     ]
    }
   ],
   "source": [
    "print len(instances), x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-136.358310911\n",
      "CLUSTER 0 0.312398884841 --------------------\n",
      "['nuclear', 'new', 'power', 'assembly', 'xymos', 'house', 'idea', 'network', 'tells', 'job']\n",
      "and the soviet union and plant nuclear weapons and new terrorist\n",
      "build a small house and plant a garden if he cleared\n",
      "of mannington and its power plant are displaced in through a\n",
      "\n",
      "CLUSTER 1 0.687601115159 --------------------\n",
      "['suitable', 'feeds', 'species', 'around', 'network', 'job', 'tells', 'idea', 'life', 'house']\n",
      "this the rain feeds the plant the plant feeds the animal\n",
      "the greater availability of suitable plant and animal species for domestication\n",
      "area provided it with more plant and animal species suitable for\n",
      "\n"
     ]
    }
   ],
   "source": [
    "best_score = -sys.float_info.max\n",
    "best_theta = None\n",
    "best_pi = None\n",
    "for _ in range(30):\n",
    "    scores,theta,pi = em(x,2,alpha=1.,verbose=False)\n",
    "    if scores[-1] > best_score:\n",
    "        best_score = scores[-1]\n",
    "        best_theta = theta\n",
    "        best_pi = pi\n",
    "print best_score\n",
    "base_rate = x.sum(axis=0)\n",
    "final_q = eStep(x,best_theta,best_pi)\n",
    "for k,cluster in enumerate(best_theta):\n",
    "    print 'CLUSTER',k,best_pi[k],'--------------------'\n",
    "    print [str(vocab[i]) for i in (-np.log(cluster) + np.log(base_rate)).argsort()[:10]]\n",
    "    for i in (-final_q[:,k]).argsort()[:5]:\n",
    "        if len(instances[i]) == window_size * 2:\n",
    "            print ' '.join(instances[i][:window_size]+[word]+instances[i][window_size:])\n",
    "    print ''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 589,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "thresh = 1.5\n",
    "base_theta = x.sum(axis=0)/x.sum()\n",
    "th_0 = (best_theta[0,:] / base_theta > thresh)*1.\n",
    "th_1 = (best_theta[1,:] / base_theta > thresh)*1.\n",
    "th_2 = (best_theta[2,:] / base_theta > thresh)*1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 593,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f414b5ac1d0>"
      ]
     },
     "execution_count": 593,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD/CAYAAAD/qh1PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8FHX6wPHPJtkku5tGJ/TeQUWkl0gTpYiCCNgPEe8U\ne0HPgh4/8fTAjqJ4imfDUw8UpCnEoxwgIMSASDFKkw6BlE3ZfX5/zBBCgJTNZifZPO/Xixc7Mzvz\nPN+ZzbOz32mglFJKKaWUUkoppZRSSimllFJKKaWUUkoppfyoC7D8POOHAuuA1cDtAc1IKaWU3zwC\nJGEU8/zswA4g1ny9DqgZ2NSUUkqdFlKKeXcC1wK2AuNbm9NSgRxgJdC7FHGUUkqVQmkK/ZdA7nnG\nx2AU+dNOYezdK6WUskBpCv2FpALR+YajgeNlEEcppVQxhJXBMrcBzYEqQDpGt82LBd/UtGlT2bVr\nVxmEV0qpoLYLaFaSGfyxRy/m/2OA8Rj98g8AizEO1L4L/FFwpl27diEiQfvv6aeftjwHbZ+2TdsX\nfP+ApiUt0qXdo/8N6G6+/iTf+PnmP6WUUhYriz56pZRS5YgW+jKSkJBgdQplKpjbF8xtA21fZVTw\nHPhAErO/SSmlVDHZbDYoYe3WPXqllApyWuiVUirIaaFXSqkgp4VeKaWCnBZ6pZQKclrolVIqyGmh\nV0qpIKeFXimlgpwWeqWUCnJa6JVSKshpoVdKqSCnhV4ppYKcFnqllApyWuiVUirIaaFXSqkgp4Ve\nKaWCnBZ6pZQKclrolVIqyGmhV0qpIKeFXimlgpwWeqWUCnJa6P3M6/Wyf/9+0tLSrE5FKaUALfR+\ntWfPHlq27EjTphdTtWotnnlmqtUpKaWUFnp/GjnyVlJSRuB2HyQnZycvvvguS5cutTotpVQlp4Xe\nj5KS1uPx3AXYgHiysq5hw4YNVqdVrqWmprJkyRJWrlyJx+OxOh2lglKY1QkEk/j4hqSkLANGAtlE\nRKykYcN7rE6r3Nq+fTs9evQnO7sJXu9hLrqoHsuWfU14eLjVqSkVVHzdow8B3gJWA8uBpgWmXwP8\nAKwD7vQ5uwrmk0/eITr6LmJjr8Tl6kCfPvUYNWqU1WmVW7feOpGjRx/k5MlE0tKS2LgxjBkz3rQ6\nLaWCjq979MOBcKA70AWYZo47bTpwCZAObAU+AVJ9T7Ni6NKlCzt2JLFu3TqqVq1K9+7dsdlsVqdV\nbqWkpCAy0BwKJTOzL7/8kmJpTkoFI1/36HsAi8zXa4FOBabnAHGAA6PDWnyMU+HUqlWLoUOH0qNH\nDy3yRejUqSN2+9sYH49UnM5P6dato9VpKRV0fC30McDJfMOeAsuaBmwAkoGvC7xXKQD++c/XaNVq\nDQ5HbcLD63PDDd246aabrE5LqaDja9fNSSA633AI4DVfNwDuBhoCGcCHGEcnP/cxlgpSNWrUYNOm\nVfzxxx84HA6qVq1qdUpKBSVfC/0qYCjwb6ArkJRvWiTGHn4WRvE/hNGNc47JkyfnvU5ISCAhIcHH\ndFRFFRISQt26da1OQ6lyKzExkcTExFItw9dOZBswA+hgDt8GXApEAe8A9wNjATewExgP5BZYhohU\nmq57pZTyC/PYX4lqt5VHC7XQK6VUCflS6PXKWKWUCnJa6JVSKshpoVdKqSCnhV4ppYKcFnqllApy\nWuiVUirIaaFXSqkgp4VeKaWCnBZ6pZQKclrolVIqyGmh96MTJ04waNAIwsNdVK1aj08/nWN1Skop\npfe68afBg0fx7bfRZGdPB7bjcAwlMfErOnfubHVqSqkgofe6sdjy5UvJzv47EAtcRk7OjSxfvtzq\ntJRSlZwWej+Kja2G8YhcACE8fCvVqlWzMiWllNKuG3+aN28eY8aMx+O5Hrt9O40bp7Ju3XIcDofV\nqSmlgoTej74c2Lx5M8uWLaNq1apcf/31REZGWp2SUiqIaKFXSqkgpwdjlVJKnUMLfQVz9OhRjhw5\nct5pXq+X/fv3k5aWFuCslFLlmRb6CiI7O5urrx5DnTqNqVu3KVdeOQK32503fc+ePbRs2ZGmTS+m\natVaPPPMVAuzVUqVJ1roK4gpU/7O0qUnyc4+RHb2Ib7/3suTT/4tb/rIkbeSkjICt/sgOTk7efHF\nd1m6dKmFGSulygst9BXEihXrycwcB0QCEWRm3s7KlRvypiclrcfjuQvjGE08WVnXsGHDhgssTSlV\nmWihryCaN2+I3b4MMM5UCgtbTvPmDfOmx8c3BJaZQ9lERKykYcOG5yxHKVX56OmVFcTRo0fp3Ply\nDh+OAkKoUuUo69YlUqtWLQDWrl3LgAHDCAnpSG5uCn36tOerrz4lNDTU2sSVUn6l59EHuczMTFas\nWMHJkyfZvXs3Ho+HIUOG0Lp1awAOHjzIunXrqFq1Kt27dz/9gVBKBREt9JXA4cOHueiibhw/3gmP\npyZ2+ycsWvQlvXr1sjo1pVQA+FLow8omFVVWXnnldY4c6UdOzkwAcnJ6MHHiX9m06b8WZ6aUKq/0\nYGwFc/jwcXJyWuYb05Jjx45Zlo9SqvzTQl/BDBt2BU7na8Bm4A8cjscZNmyQ1WkppcoxLfQVzODB\ng3nxxUlUqTIYp7Md113XgGnT/s/qtJRS5ZivB2NDgBlAByALuB3YlW/6ZcA0c/n7gJuB7ALL0IOx\nSilVQoG8e+VwIBzoDkzCKOp5eQBvA7cCvYDvgMY+xqlQRIRXX32D3r2HMnLkzWzfvt3qlCq0EydO\ncOed99Gz52Duv38S6enpVqekVIXk6x79NGAt8Jk5vBeoZ75uCbwBbAPaAQuAF8+zjKDbo3/88cm8\n8srXZGQ8gc22nejol0hO/oH69etbnVqFk5OTQ8eOvdixox1ZWcOIjPyEjh2PsnLlYr0+QFVqgdyj\njwFO5hv25FtWdYw9/deA/kA/4HIf41Qor7/+JhkZnwHXIPIoWVlD+fzzz61Oq0L68ccf+e23U2Rl\nvQMMw+3+F5s2bWXXrl1FzquUOpuv59GfBKLzDYcAXvP1UWAn8Is5vAjoBCwvuJDJkyfnvU5ISCAh\nIcHHdCqHTZs2MWPGu3i9wp133kqnTp2sTkkpVcYSExNJTEy0JPa1wHvm664Y3TOnhQO/Ak3N4S+A\nK8+zDAk2jz32tDidHQW+lJCQ5yUmppbs2bPHL8teu3atOJ3VBZ4TeF6czuqyYsUKvyy7PMrOzpZ2\n7bpIRMQ4gXkSGTlauncfIF6v1+rUlLIUp+9sWAK+dnbaOHPWDcBtwKVAFPAORlfN8+b7VgH3n2cZ\nZs7BQ0S4666JfPHFPOLiYpgz5yMuvvhivyx7+PAbmTevK3C3OeYdBg5czOLFJesa+vHHH1m/fj31\n69fniiuuKNf93ampqUyaNJmfftrOZZe1Z8qUJ3G5XJbl4/V6WbhwIX/88QedO3emQ4cORc+klJ/5\n0kdvJYu/F/1v6tSpAi6BUQIXicsVL6dOnfLLsgcOHCHwoYCY/z6XXr2GlGgZM2fOEqeztjidt4nL\n1VZGj75N95CLyePxyNCh10tU1EXidN4qDkdN+de/PrI6LVUJ4cMevZWsXl9+FxISJ7DYLMQegW4y\nbtw4vyz73//+XJzOhubyl4rT2aREhSYrK0siIqIEfjHzyxCXq3lQd//40+LFiyUqqr2A21x/SRIZ\nGSMej8fq1FQlgw+FXm9q5kdebyZGDxYYx6e7snfvVr8se+TIEWRmunn++cmICA8++AQ33jj2rPck\nJSWxfPlyqlSpwqhRo4iMjMyblpqaCtiBFuYYB6GhbThw4IBf8gt2xnpqD0SYY9qRk5NFZmampd1J\nSpV3Vn8x+l1cXEOBvwjkCGwTqCpvvPFGQGJ/9dVX4nDUkIiIv4jL1V86dOgmmZmZedO9Xq80aNBK\nbLZXzF8bK8TprC4pKSkBya+i27ZtmzidNQTWCXgkJGSqtGp1qdVpqUoI7bqx1tatWyUqqo5AiEC4\njB17c8Bi167dVCDR7FbwitM5SGbNmnXWe7Zv3y7Nm18iNluoxMbWkm+++SZg+QWDL7/8j0RHV5eQ\nkDBp0+Yy/ZJUliCAZ934g5lz8ElLS8PpdBISErh7xjmdcWRm7sS4Xg3Cwh5iypQaPProo+e8Nzs7\nG7vdXq7PuCmvRIScnBzCw8OtTkVVUoG8MlYVIioqKqBFHqBPn/6Ehz8GnAI2YLd/dMEL0MLDw7XI\n+8hms2mRVxWOFvog8ckns+jT5whhYTWJixvKO+9Mp0uXLlanpZQqB7TrRimlKhDtulFKKXUOLfRK\nKRXktNArpVSQ00LvZ7m5uWzcuLFcXHHq8XjYu3cvmZmZVqeilLKQFno/Wrt2LU5nTS69NIH4+AZc\nfvn57s4cGElJSdSp04wWLTpRpUot3n//A8tyUUpZS8+68aOoqLqkp98BPAUcADoxdepEJk2aFNA8\nRIQ6dZpx4MBk4CZgGw5HHzZsSKR169YBzUUp5V961o3F0tOPABMxtkE8cL0lT4Y5ceIEx44dxijy\nAK0IC+vN5s2bA56LUsp6Wuj9KCzMBSwzh7KBZTRv3jzgecTExBAWFgb8YI5JxevdQMOGDQOei1LK\nelro/eitt14EbgF6A82oUeMU06ZNC3geoaGhfPzx+zidVxETMxinswO33HIN3bp1C3guSinraR+9\nnyUnJzNnzhzq1KnDhAkTAn7Pm/x+//13Nm/eTN26dbn00kuLnkEpVe750kevhV4ppSoQPRirlFLq\nHFro/UhEmD37X1x99Q3cfvvd7N692+qUfJaZmckTTzzDkCFjeOqpv+F2u61OSSnlI+268aP/+78X\neO652WRkPEJo6HZiY2ezZct6ateubXVqJeL1eunVaxAbN8bgdg8nMvILLrvMTWLiAkuPOSiltOvG\ncs8/P42MDDfwIB7P26Sm1uSzzz4LSOzc3FyeemoK7dv3JCFhKBs3bvR5WVu3bmXz5p243Z8CN+J2\n/5sNG7awfft2/yWslAoYLfR+lJaWgfEov8+BKXg821izZk1AYt9//ySmTVtKcvIUvv9+KH36DGLX\nrl0+LSs3NxebLRwINceEYrOFk5ub67d8lVKBo103fmSzOYDtQH1zzB00b54YkD3h6OgapKVtzIsd\nHv4Xpk5txgMPPFDiZeXk5HDRRd3ZubMrOTkjCA//jObNN7Bp0yrzQiyllFW068ZyIUBavuFTOByO\ngEQODbWfFTskJN3nomy321mxYhGjRmXRvv2TXH+9hxUrFpV5kd+3bx9XXnkdDRu2Z9iwMRw8eLBM\n4ylVWegevR+1bNmB7dtPYtzU7GdgBsuXL7jgQ7r96e9/n8azz84iI+NhQkO3Exf3EVu2rKdWrVpl\nHtsf3G43LVt2ZN++6/B4rsFu/5hGjZayZcs67Ha71ekpVW7oBVPlQL9+A1m9+iccjlA++OAthgwZ\nErDYH374Mf/5zyJq1qzKX//6EPXq1QtY7NL64Ycf6Nfvdk6dOn3jNSEqqgVr1sylbdu2luamVHmi\nhV5VWD/99BPdug0jPX07YAfcOByNSU5eRZMmTaxOT6lyI5B99CHAW8BqYDnQ9ALvexuY6mMMVYm0\nbduWbt064HAMA97A6RxM//59aNy4sdWpKVXh+VrohwPhQHdgEnC+WzROANoButseJDweD8OHD6d+\n/Yb06dPHr48oDAkJYcGCfzNlykBuvvknpk4dzpdffnh670UpVQq+/hVNA9YCp68G2gvk7xDuDowD\n/gu0Ah47zzK066aCqVGjEUeOhAI3AN8QHv4bGRkHCQ0NLWpWpZSfBLLrJgY4mW/Yk29Z8Rinndxd\n0mRU+bVmzRqOHDmI8TCTZ4FVZGeHM2XKFIszU0oVxdcTo08C0fmGQwCv+XokxuWh3wC1ASfGuYbn\nPJ168uTJea8TEhICchqi8s3vv/+OsSmrmGMigHj27dtnXVJKVQKJiYmlfiSpr3vc1wJDgduArsCT\nwODzvO8WtOsmKGRmZuJ01gIeAMYDC4F72bz5f3To0MHa5JSqRALZdfMfwA2swuivvx8Yg1EBCtJq\nHgQcDgfz5n1ISMhrQHNstoeZNu1vWuSVqgD0PHqllKpA9F43FvN6vUyd+g86derHlVdeR3JycpnE\nWbVqFX37Xk3nzgOZOXMW+oWplCqM3orQjx566K/MnJlIRsZkbLZfWLWqL0lJ62jUqJHfYmzcuJGB\nA4eTkfF3oCZbtjyC2+3m3nvv9lsMpVRw0a4bP4qJqcmpU+uARgCEh9/J1KktfLpV8IXcc8+DvPZa\nHMbxb4BVNGlyN7t2/ei3GEqp8ku7biwWEhIKZOcN22zZfr+YKCwsFJstO98Y/8dQSgUXLfR+9OCD\n9+B0jgA+JiRkMg7HIkaNGuXXGOPH34bT+RbwIjAbp/NPPPbYRL/GUEoFF+268SMR4c03Z/LPf35M\n/fp1eemlqX7tnz8tKSmJ5557mbS0TG677TpGjLgWMA4GHzhwgJiYGKKiokq83EOHDrF7924uvvji\nEj9kJDc3l6SkJOrUqRPwh6GXtt0VVWVtd2XnS9eNlSTYzJs3TyDa/GeXRo1aByz27t27pVmziyQy\nsobY7U6ZPPm5Es0/aNDVAuECMRIaGifff/99sedds2aN2O1V8tqdkDCopOn7rLTtrqgqa7uVCBXs\n2iSr15ffhYRUEZgs4BXYLxAvf/rTnwISu3PnvhIa+mxebJerqSxZsqRY87766qsCNQT2CIjA8xIZ\nWbPYsV2uOgXaXUemTp3qa1NKpDTtrsgqa7uVb4Ve++j9yOtNByZi/KqKB0azYsWKgMROSlqPx3NX\nXuysrGvYsGFDseb97rvvMG5RdPoGpHfjdh8tduz09COc3e7rS31vjuIqTbsrssrabuUbLfR+FQks\nM19nA8sC9uCM+PiGZ8WOiFhJw4YNizVvixYtMJ4f4zbHLCckpPh9vmFhLgq2u3nz5sWevzRK0+6K\nrLK2W1U8Vv8C8rvHH39cwCHQU6Ce2O3VJCMjIyCx16xZI9HRNSU2dpC4XC3lqqtGSm5ubrHmzcnJ\nkfj4FgJ1BXoLOOW1114rduxZs2YJOAV6CdSXGjWaSFZWlq9NKZHStLsiq6ztVr513ehZN36WmJjI\nm2++Sf369fn73/8e0HPcDx48yLp166hatSrdu3cv0dOZvF4v77zzDnv37mXEiBFcfPHFJYqdnJzM\nnDlzqFOnDhMmTCAkJHA/FkvT7oqssra7stOHg1ssPT2dceMmsmjRImJiqjBjxt8ZMmSI1WkppYKI\nFnqLjRx5M/PnZ5OV9QKwHYdjLKtXLynx3rFSSl2I3gLBYgsXzicr61WgAdCfnJwbWLJkidVpqUpC\nRNi8eTMrVqzg1KlTVqejyhG9e6UfhYc7ychIAWoCIPILMTEtrU1KVQoej4cRI25i6dLV2O3xhIXt\nY8WKxbRu3drq1FQ5oHv0fhQZGY7xhMUngOvwetfgcrkszkpVBh9++CHffvs7GRnbSE39H8eOPcYN\nN0ywOi1VTmih96MTJ44CszF+KPXCZruZvXv3WpxV+eXxeHjiiWdp1KgDbdp0Zf78+VanVK6JCNOm\nvUKTJpfQosVlfPDBh3nTfvllB+npAzGu5QCRoezatcOiTJU6w7oTUcvIRRf1EJvtJfM2AkfE5Wol\n8+fPtzqtcuuxx54Wp7ObwA8C88XprCWrVq2yOq1y6/XX3xSns63AKoHvxOlsIPPmzRMRkTlz5ojL\ndbHACQGR0NAp0r37QIszVmUBvdeNtXbs2CF16zaXqKimEhERK/ffP0m8Xq/VaZVbdeu2FthofjGK\nwP/JPfc8aHVa5dYllyQILMy3vmbJ8OE3ioiI1+uVO+64RyIi4iQqqok0bNhafv/9d4szVmUBHwq9\nHoz1o2bNmvHrr8n8+uuvVKlShVq1almdUrnmdDqBw3nDoaGHiIqKsS6hcs7lcgKH8oZttkNERTnM\n1zZmznyFyZMncfLkSZo0aYLdbrcoU6XOsPqLsUw8/PDD0qhRe+nYsZts3rzZsjw8Ho/cdNMt0qBB\nO+nZ83LZt2+fZbmcOHFCWrVqL05nvLRo0UaOHj0qIiJffPGFOBy1BV6Q0ND7JS4uXnbv3l2muZw6\ndUqefvpZGTNmnLz55kzxeDxlGq8whw4dkocemiRjx94un346p8j3JyYmitNZXWCK2Gx/laio6pKc\nnByATFV5gt4CwVrXXHMdc+euBv4KbMVmm83PP6+nZcvAn2LZvn0XkpPTgbuBROz2pRw6tIu4uLiA\n5uHxeHA4apGT0wa4Afg3oaEbycg4QHh4ON9//z2ffvofoqOd3H33nTRo0KDMcsnOzqZTpz5s396I\nrKwEnM73GT26I++++0aZxbyQ48eP07btZRw5cgU5Oe1wOl/hiSfG8dhjDxc63/r163n//Y8JDw9j\nwoRxlny2lLX0wSMWs9miBH7K14c6UkaOHBnwPI4fPy4QJnDUzMMrcIk8+eSTAc9l9uzZAnECbjOX\nbIEa8tJLLwU8l6VLl0p0dCdzfYjACbHbnXLy5MmA5/LOO++I0zki32dlu7hc1QKeh6p40PvRW0vE\nA1TJN6Y6brf7Qm8vM0ZMGxBtjrEBVcjMzAx4LseOHQMcQLg5Jgxw+XTlptvtRkrxK9DtdmOzxXFm\nZ8iFzWYnOzu7sNnKhNvtxuvN/1mpSm5uVsDzUKqsWf3F6HctW15q3qJ4jcAHAk5ZtGiRJblUqdJA\n4DqBdQLTxWZzyrZt2wKex4kTJ8xHDN5nnkb5iECUHDp0qNjL2LZtmzRu3E5CQuwSE1NTvvnmG59y\nOX78uFSv3kBCQl4QWCcREbdK9+4DLDkz6tdffxWXq7rAuwJrxeG4Um68cXzA81AVD3p6pbXS09Pl\n0kt7SVhYDXG56smsWbMCGn/z5s3y8ssvy+zZsyUlJUWaNr1YwsJqSGxsQ8u+cEREvv/+e4mIqClQ\nRcLDa5bokXcej0fq1WshNtsMs8tlpTid1eW3337zKZedO3dKv35XS5Mml8iNN94hqampPi3HH9av\nXy9duw6Qpk07yj33PFLie/jn396ZmZlllKUqb9CDsZXX119/zfXXj8PrvY6wsO00bZrO2rXLiIyM\ntDq1Ujlw4ACNG7fH7T5zGmZMzDDef/82rrnmGgszs1awbm9VNL17ZSV2xx33k5n5b7Ky3iA9fQk7\nd8by0UcfWZ1WqcXFxSHiBk5fzp+Ox7OV2rVrW5mW5YJ1e6uy4WuhDwHeAlZjPGy0aYHpY4A1wErg\nTSrQqUAVVWrqEaCtOWQjO7stR44csTIlv4iMjOTVV1/G6eyNy3ULLlcnrr22H127drU6NUsF6/ZW\nZcPXQj8c4zSK7sAkYFq+aQ7gb0AC0BOIBSrNY5YSExO5/vrreeCBB/B4PAGL26dPf8LDHwNOARuw\n2z8iISHB73E++OADRo4cyfTp0/2+7Au5445xrFq1kFdf7cO8ea8ze/Zblf6xeb5s78zMTCZOnMjY\nsWNZt25dINJUFdw0YFS+4fy3aLQBNfINfwb0P88yrD6m4XeTJk2y7OHgx48flwEDhktYWKTExcXL\nhx9+7PcYvXpdLhAlkCBQRRo0aOX3GKp4Srq9jx49KqGhcQJNBLoJOC25lkGVHgE8GPsO8AWwyBz+\nHWgMeAu8byIwCBh8nmWYOQcPmy0WeBcYCWQDnRk4sBaLFy+2NjE/2L59Oy1bdgA2Aa2Ag0AL3n77\nH4wfP97a5FSRLrnkEjZtqgIsBUKBt7DZHsfrPWZxZqqkfDkY6+tNzU5y5mocMLqAvAWGXwCaASMu\ntJDJkyfnvU5ISCiTrobAcgN9zdfhQF9SUoLjHutr1qwBYjCKPEAtoCVr1qzRQl8B7Nt3GLgeo8gD\n9Eek8NstqPIhMTGRxMRES2JfC7xnvu4KLCgw/R3gNQr/1rH6F5DfhYRUEZhsnu+9XyBe/vSnP1md\nll8cOnTI7Jb62rxk/wcBpyxYsMDnZbrdbrnttr9IdHRNqVGjkbz77nv+SzjAdu/eLV279hens6q0\nbHmpbNiwweqUzjJo0CCBZgJHBDwCE8Vur2F1WsoHBPCCKRvG2TSrzH8tMM60GQ9cAngwzsY5/W/4\neZZh9fryu3nz5plXgUYL2KVRozZWp+RXjzzyiFnsYwQi5NprR5RqeX/+8/3icFwpsFvgB3E665fo\nYqryIjc3V5o27SChoc8KHBT4UGJja8vhw4etTi1Pbm6u1KjRUMAu4BKbLUZWrlxpdVrKB+gFU9Zz\nu918+eWXtGvXjg4dOlidTomkpKTw888/079/f8LDw8/7nrS0NBITE+nSpQs1atQ473uKq27dVuzf\n/znQzhzzAnfddYDXXw/cGT3+sHv3blq16kJm5n5O/0nFxvbn008fYtCgQdYmV8DPP//MH3/8QZ8+\nfQgNDS16hgvIzc0lKSmJOnXqVPprGgJNL5iy2J49e2jfvivjxt1Hp07deOaZqVanVGwtW15EkyYt\nGTx4DBER1ZgzZ8553xcVFcWQIUNKXeQBYmPjgJ15w3b7TqpXD+xtlP0hJiYGjyeNMw9RySY3d3fA\nbwldHK1bt6Zv376lKvJr167F6azJpZcmEB/fgMsvv9KPGapgY/UvIL/r3Lmv+fPd6KN3uZpWiK6I\ne+65R6C6wB6z//15sdniyjzut99+K05ndQkLe0AiI8dK7dpNylV3R0lMmvSUuFwtxWZ7TFyu7jJ4\n8HVB+xhJl6tOgWNRdWTq1KlWp1VpoF031nI4YnG7U4CqAISFPczf/laNSZMmWZtYEdq3b09yci9g\nhjkmHYhFJLfMY//0008sWLAAp9PJjTfeSNWqVcs8ZllZsGABGzdupHHjxowZM6ZUe83lmc0WAfzB\n6c85PMAVV2xl0aJFhcyl/CWQp1eq84iPb0hKyjJOn0cfEbGShg3vsTqtIjVt2pTk5OUYp4dGYhw/\ndwQkdvv27Wnfvn1AYpW1wYMHM3jw+S4ZCS5hYS5yc898zmEZzZv3sjgrVRjto/ejf/3rTeBmoAEQ\nT4MGWYwaNaqIuaz32WefERFxFOOyh97AKO69d1zA80hLS+PBBx+jf/9refTRJy15UIoq2ltvvQjc\ngvFZaUZOusj0AAAVXklEQVSNGqeYNm1aEXPB6tWrGT78RoYMGaN7/wGmXTd+FBdXm9TUaIxb/WwD\nXuCjj2YxduxYizMrmsfj4bHHHuP333/n9ttvZ8CAAQGP36VLX7ZsqY/bfTWRkZ/RqdMpvv/+G0JC\ndH+kvElOTmbOnDnUqVOHCRMmFLmN1qxZQ79+Q8nIeAaIxOF4gk8/fYthw4YFJuEg4kvXjRZ6P7LZ\nXBi3CGhujrmZ+Phv2b9/v4VZVQybNm2iV69RpKVtw/ihmYvT2Zgff/yOFi1aWJ2eKqXRo//EnDkX\nA6e7Mj+na9dZ/O9/umdfUnp6peUEsOcbDivVM05La+bMmfTr14+JEycG9E6avvB6vRiX55/+/IZg\ns4Wa4/1HRFixYgUffvghW7Zs8euyS8rr9bJ48WI++ugjUlJSLM2lrHk8Xs7+27D7fduWhUOHDjFn\nzhzmzZtnyfOfg4F15yeVEZvNJdBBYIHAdAGHvPDCC5bk0r//IPNOk8ME6onLFS+5ubmW5FIc2dnZ\n0qbNZRIefqfAIomIuE06duzl95xvv32iuFzNJCpqtDgcNeW992b7dfnFlZubKwMGXC1RUR0kOnqU\nOJ3VK8SpuL5avny5OBw1Bf4l8Lk4nQ3kk08+tTqtQm3ZskXi4uIlOvpqiYrqJS1bdrT00ZOnoc+M\ntRa4zFvA1jfvK9JM6tSpE/A8MjIyBMIFNpvnxacJ1JV777034LmUxLFjx+S22/4il17aV+644x6/\n/1GtXbtWXK5GAifN9fKzREREi9vt9muc4pgzZ464XF0Fss1clkqtWk0CnkcgLVmyRPr0GSrdu18p\nc+Z8ZnU6RerZ80qx2V4zt49XIiJukKefftbqtHwq9Np141c5wHxgN8aj74Zy4sSJgETOzc3lqaem\n0L59T3r2HIjRBXL6Fgwu4CK2bdtWJrHdbje9ew8kMrIOcXGNfH6kXZUqVfjnP99g/frvmDnzFWJi\nYvya5759+wgNbc+ZG6+2wmaL4Pjx4z4tL/86T0gYysaNG0uUS27uZZzpzujG0aP7AhLbKgMGDCAx\n8StWrfqGUaOuszqdIu3duw+R7uaQjaysbqSkFH8bKYPVX4x+B7ECD5lXDKYIVJcePXoEJPbddz8o\nTmdvgeUCM8xfF6+buWwUiJJZs2aVSez27bsIdDZjvyXglO+++65MYpVGSkqKOJ3VBdaae2nvSnx8\nU/F4PD4t7+x1PlOiomrIzp07izXv6tWrxemsK7BDwCuhoU9Jp04JAYmtiueWW+6UiIgbzF9dh8Xl\nukTee+99q9PSrhurzZw507y7Y6yASyIjY8o0Xk5OjjzzzHPSs+dg8+lB68wCJhIa2tvsow8XCJch\nQ4aVWR5Gm3fnxYbxMmxY6eMdP35cJky4V3r0uEruu+9RSUtLK/Uy582bJ05nFQkLc0q9ei0kOTk5\nb9r+/fvl5psnSI8eV8mTTz4r2dnZhS4rKqr6We0OD/+zTJs2rdi5zJgxU8LDXRIW5pC2bTvL3r17\niz2vL7E3bdokV189VhIShsmsWf8M2ls0+MupU6ekX79hEhYWKWFhEXLffY+Wi3WGFnprdet2uUBr\ngS8Fpgo4Zc2aNWUW76abxovT2VdgrsDD5rGB4wIikZE3y/Tp02Xr1q2SlZVVZjmInD42sTVfoR8t\n1157bamWmZ2dLe3adZGIiHEC8yQycrR07z7AL39oHo9HUlNTz1rWyZMnpW7d5hIW9rDAPHE4rpCR\nI28qdDmxsfFntTsy8mZ55ZVXSpRLbm6unDx5ssTzlDT2tm3bxOWqLvCywL/F6Wwt06Zd+P1er7fI\nL7rKIi0trcz/hkoCLfTWAqfAznwF7xa5+uqryyRWVlaWhIaG5zuwKGI8y3WChIY+KtWq1ZMDBw6U\nSeyCBg0aJtBQ4F2Bh8Rmc8lPP/1UqmWuXbtWoqLamF1PIpAjTmdd2bFjh5+yPtt//vMfiY7ul29d\npktYWGShvyKef/4f4nS2Eng3IOs8KytLxoz5k4SGhovNFiFhYU2KHfvxx58Qm+2RfO37QerWPf8z\nf+fM+UyioqqLzRYqHTp0lz179pRVk5QP0IOxlVtERDidO29l/Hg3mzb9j1q1agUk7sKF87jzzsHU\nqfMP2rRZztq1y2nXrl3RM1YAn3/+Oa1adaRly4588MEHZ0179NEHmTnzSa69NjEg6/yJJ/7G3Ln7\n8XiOILIHmw0uuuiffo2dnJzMbbfdTVraYkSy2LJlIIMHX++H7FVlZfUXo9+d3XXznIBT1q1bV2bx\n8nfdhIU9JrVqNZbjx4+XWbxAKsuum/M503XzkMBccTiukE6dupm/0v5P4HkBl7zxxhtlEr84OnTo\nJbAs3175BzJ48OhizXtu102r83bdvP322+J03pYvRq6EhIRpN045gnbdWMvj8cjAgVdJZGQ1qVq1\nrsyfP79M4+U/GHvDDbeX6GBeRXDixAm58877/HowtjB//PGH3HzzBOnZc7A8+eSzUrt2S4HX8hW9\nt6VataZlmkNhrrzyOgkJeTEvn7Cw++TPf76v2PMX52DsV199JVFRHeXM+f0bJCqqWrk4CKkM6P3o\nrbV06VKGDx+LzXYVNtsO2rZ18N//LrzgY/lU+VatWjOOHXsGuMEc8wUxMQ+RmppiST47duygc+c+\n5OT0BDKJifmZH39c5dfuIq/Xy5VXjmT16t14vR2ABbz33usV4rz3ykJvamaxevVasW/fq8BAwIvL\nNYDXX7+JW2+91eLMlC8efPBhpk//APgXxuGsm7nzzmt48803LMvp4MGDLFq0iLCwMIYMGUJsbKzf\nY3g8Hr755hsOHDhAt27dguZ4S7DQQm8xpzOOzMxdQDWg4jxhSl3YhAl/Yfbs/yACY8Zcxfvvv2t1\nSqqS00Jvsb59h7FyZWNycqYBu3A6+7F48af07NnT6tSUUkFCb1NssTlz3uXSS7cSEuIgMvIypk+f\nrEVeKWU53aMvA9nZ2djt9tPfvEop5TfadaOUUkFOu26UUkqdQwu9UkoFOS30SikV5LTQK6VUkPO1\n0IcAbwGrgeVA0wLThwLrzOm3+5ydUkqpUvO10A8HwoHuwCRgWr5pdmA6MADoA9wB1CxFjhXKzp07\nadKkCXfccUfAY7vdbj7++GOSkpIAOHbsGE899RTr16/Pe09aWhrz58/nwIED58z/ww8/8PLLL5OR\nkXHBGKmpqcyfP5/Dhw+fNd7r9bJ//37S0tJK1Ybs7Gz27t1LTk4OAFlZWWcNn4+/YoNxm97p06dz\n5MiRYr2/tLGLs85L4tSpU/zxxx+ISN66zM7Ozpv+888/s2zZMjwezznz7t69m6SkJLxeb5FxcnNz\n2bhx43k/R4U5duwYhw8fpqgz7o4ePXrBbeDP7a0KNw0YlW94b77XHYCF+YanAyPPswwrbvxWpuLj\n483H90UL2AWcAYs9a9ass2KHhkaZt9iNEQiX6Oiq8txzz5njYgXCJSGhf978Dkc1M+cYgSh5+umn\nz4nx8MMPy5lHJUbI8OHGU6R2794tzZpdJJGRNcRud8rkyc/51IaFCxdKVFQ1cThqS0xMTZk6daq4\nXFXF4agtsbG1JDEx8Zx5/BVbRKRWrXpiPHoxRsApY8eOLfT9pY1dnHVeXF6vVx555Emx250SGVld\nGjRonbfuoqKqyYIFC6RGjUZm+1xis8XIypUrRcS462rbtp3zpjkctWT79u0XjLVmzRqx26vkfdYS\nEgYVmV9OTo6MGnWLhIdHS3h4rCQkDD7v3UizsrJk2LDR5vtiZNCgayUzMzNvuj+3d0VFAG9T/A4w\nKN/w75z5ddAT+DTftGeAcedZhtXry++MAjjZfCrSfoF4CVQ7z41dW+D0E4V+Nv8owwU+Nscly+kH\nhvfu3VugusAec9rzAmc/7/bQoUNmkf867wlF4JQFCxZI5859JTT02bzYLldTWbJkSYnyP3z4sHm/\n9JXm8v9txlttDi+W6Oia5xQHf8QWEbnrrrvMdXT68Xz/EnAWeh/20sQuzjoviblz54rL1VrgkJnP\ngwI9zGX/V0JC7ALNBI4IeAQmit1eQ0REJkyYIMZzFI6Z0yZI7dotLhjL5apT4LNWR6ZOnVpoflOn\nvihOZz+BdIFsiYwcLX/+8/3nvO/JJ58Vh+MqgUwBtzgcw+Whhx7Pm+6v7V2REcAnTJ0Eogss5/Tv\nvdQC06KB4z7GqWAygYkY1zLEA6MtjH09UMWc1grohbGJxpjj2gLdmDt3Lj/++CNwHVDPnHY3kH7W\n0hcuXAhEAUPMMZ2AtnzxxRckJa3H47krL3ZW1jVs2LChRNlv376dsLAmQA9zTE2gGdDNHB4IVOG3\n3347az5/xAb4+uuvMfZRWptjbgS8LFiw4ILzlCZ2cdZ5Saxd+wPp6aOAGmY+9wK/mlN74fXaMfa3\nqmH8ud5DTk4mAP/973qMQ2lVzGkTOXTo7K65/NLTj1Dws5aYmFhofitWrCcj4zbACdhxu8ezatW5\n62rFivVkZo4DIoEIMjNvZ+XKM+/z1/aubHwt9KuAq8zXXYGkfNO2Ac0xPjXhQG/gf+dbyOTJk/P+\nFfVBqRgigGXm6+x8r62KnWsOpwIbMHYENprjjgM/0rFjR2rWrIlxTN1tTlsOOM5aeteuXTG+37eZ\nYw4Cv9C1a1fi4xueFTsiYiUNGzYsUfb169cnK2sXZ3oBBUgB9pvDO8nJOUB8fPxZ8/kjNkCbNm0w\n1lGqOWY94OXyyy+/4DyliV2cdV4STZo0wun8L2e2+XcYRRhgDzabG1gMnO6b/xabzQ5A06Z1C0xb\nhtPpumCssDAXBT9rzZs3LzS/Fi0aEh6+jNM7o6Ghy2nW7Nx11bx5Q+z2M+8LC1tO8+Zn3uev7V2R\nJCYmnlUrA8kGvIlR8FcBLTB2Fceb04dgnHWzHvjzBZZh9S8gv7PZbGZ3Q0+BegLRsnbt2oDEvuWW\nWwrEdpnDfcwughgZNuxqc3wfgWpSvXpjERFJT083+4nrCvQScEifPn3OidGr1+ViHAdIEIiTBg2M\nh0uvWbNGoqNrSmzsIHG5WspVV42U3NzcErfhhRdeEoejlsTEDBaHo6YMG3adOJ21JSZmsDidNeWt\nt945Zx5/xRYRCQmJFahmrh+nNG/estD3lyZ2cdd5cWVnZ8vllw+WqKg2EhNzhURGVpWIiCrmuqwl\nTz89RUJD4wSaCBiPSHzppZdEROT48ePictUWaC7QTWw2lyxYsOCCsYzjQU4z7/pSo0YTycrKKjS/\nEydOSMuWHSU6urNER/eWOnWanfeJaEeOHJEmTdpLdHQ3iY7uIQ0atDrroef+3N4VFfooQeu1adPm\n9IaQo0ePBjT2F198IT179pTRo0dLTk6OPPjgg1K3bl3p27dv3ns+++wzGTlypDz33NkHsbKzs2XA\ngAHSunVrmTFjxgVjzJ49W0aMGCHTpk07a/yBAwfkq6++kpUrV5bqsXNbtmyRuXPnyrZt20REJDk5\nWebOnSu//PLLBefxV2wRkZEjR0qrVq1kypQpxXp/aWIXd50Xl8fjke+//16+/vprOXTokGzbtk3m\nzp0rW7ZsERGRjIwMufvuu2XMmDHn7IBkZmbKSy+9JJMnT5adO3cWGeunn36SJ554QmbMmCEej6dY\n+WVmZsqSJUtk4cKFcurUqQu+LyMjQxYvXiyLFi067wFbf27vigh9lKBSSgU3vamZUkqpc2ihV0qp\nIKeFXimlgpwWeqWUCnJa6JVSKshpoVdKqSCnhV4ppYKcFnqllApyWuiVUirIaaFXSqkgp4VeKaWC\nnBZ6pZQKclrolVIqyGmhV0qpIKeFXimlgpwWeqWUCnJa6JVSKshpoVdKqSCnhV4ppYKcFnqllApy\nWuiVUirIaaFXSqkgp4VeKaWCnBZ6pZQKclrolVIqyGmhV0qpIKeFXimlgpwWeqWUCnK+FHoH8AXw\nX2ABUP0877kfWGP+e8rn7JRSSpWaL4X+z8BmoDfwAfBEgelNgLFAN6ArMBBoX4ocK6TExESrUyhT\nwdy+YG4baPsqI18KfQ9gkfl6EdC/wPTdwBWAmMN2INOn7CqwYP+wBXP7grltoO2rjMKKmD4OuK/A\nuIPASfP1KSC2wPRc4BhgA14ENgI7S5emUkopXxVV6N81/+X3BRBtvo4GTpxnvkjgn0Aq8JfSJKiU\nUqp0bD7M8wBGgX8GGA30Au4qsMxFwHfAC4UsZyfQ1If4SilVme0CmpVkBl8KvQOYDcQDWRgHXg9h\nnGmzEwgFPgH+l2/5j2GcgaOUUkoppZRSSimlyqVgvNAqBHgLWA0s59xjDkOBdeb02wObml8U1b4x\nGNtqJfAmvnUFWqmo9p32NjA1UEn5UVHtuwzj73EF8CkQHtDsSqeotl0D/IDx93dnYFPzqy4Y7Suo\n3NaWBzhTvK8HXi4wvQnGhjldLFZS/i+0uhbj7CIwNsjcfNPswA6M00/tGBulZkCzK73C2ufAOCYT\naQ5/jPHhq0gKa99pEzD+mJ4LVFJ+VFj7bMCPGH93AOOBloFLrdSK2nYpQBxn/x1WNI8ASRifv/xK\nXFsCea+bYLzQKn+b1gKd8k1rjVEIU4EcjC+u3gHNrvQKa58b4+pntzkcRvnfXgUV1j6A7kBnYCYV\n79cKFN6+FsBRjB2wRIyi+EsgkyulorZdDkabHBjbTqh4dmJ8oRX87JW4tpRVoR8H/FTgXyzFv9Dq\nH1SMC61iONMmAA9n1mkMxoY47XxtLu8Ka58Ah83XEwEX8G3gUvOLwtoXj/EL9G4qZpGHwttXHeOL\n7DWMna5+wOUBza50CmsbwDRgA5AMfF3gvRXFlxh1saAS15ayKvTvYnS75P+XSvEutPoIo2hUhAut\nTnKmTWCsT6/5OrXAtGjgeIDy8pfC2nd6+B8YRWJEAPPyl8LaNxKjGH4DPIpxGvHNAc2u9Apr31GM\nHalfMIrJIs7dKy7PCmtbA4wv6IZAI6AWxvYMFiWuLYHsulkFXGW+vhLjIFB+NmAesAnjxmkV4adW\n/jZ1xehPO20b0ByognGQqzfGtQUVSWHtA6NLIwLjwJebiqew9r2GUfguB57HOAbxQUCzK73C2vcr\nEMWZg5i9MPZ+K4rC2haJsYefhVH8D2F04wSLcl1bHMBnGEf4v+XMwYP7MQ7iDcfo412GcZR5OcYG\nLM9sGGebrDL/tcA4E2W8OX0IxoGS9RhfXhVNYe27BOOPaXm+f8OtSdNnRW2/026hYh6MLap9l2P0\nb68DXrIiwVIoqm33Y5zcsQJ4j6Jv91JeNeLMwdhgqi1KKaWUUkoppZRSSimllFJKKaWUUkoppZRS\nSimllFJKKeWb/wedj0ZBhg1MCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f414baacfd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter(np.dot(x,th_1)/x.sum(axis=1),\n",
    "        np.dot(x,th_2)/x.sum(axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 467,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(172, 58)"
      ]
     },
     "execution_count": 467,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 468,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 58)"
      ]
     },
     "execution_count": 468,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_theta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 514,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f414d1a2c50>"
      ]
     },
     "execution_count": 514,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD/CAYAAAD/qh1PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4U9UbB/Bv0qx7k3TRQdl7o+yyKbsCQkGmILJVlrJE\nf8hUkQ0ylI2AyJQ9BaHsLRukbJC9W7qTfH9/JJQU2tKmaQvlfJ6nj01yz7iRvrk559z3AIIgCIIg\nCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIgCIIT+QPYkcDzHwI4BGAfgC7p2iNBEATBab4GcBLWYG5P\nDeACADfb74cA+KRv1wRBEITnlKkoexFAMwCKl54vanvtKYBYAHsAVE9FO4IgCEIqpCbQrwRgSuB5\nV1iD/HNhsF7dC4IgCBkgNYE+MU8BGO0eGwE8ToN2BEEQhGRQpUGd/wIoCMADQDiswzZjXz4of/78\nvHTpUho0LwiCkKldAlAgJQWccUVP23/bAOgK67h8XwBbYJ2onQPg9suFLl26BJKZ9mfo0KEZ3gdx\nfuLcxPllvh8A+VMapFN7RX8VQGXb74vtnl9v+xEEQRAyWFqM0QuCIAhvEBHo00hAQEBGdyFNZebz\ny8znBojzexe9vAY+PdE23iQIgiAkk0KhAFIYu8UVvSAIQiYnAr0gCEImlxbr6AVBEN5pR44cwaFD\nh5AtWzY0btwYSmXGXlOLMXpBEAQnmj17Hr78chAslg+hUh1F9ep5sW7dUqcFe0fG6EWgFwRBcBKz\n2Qy93h3R0UcAFAYQA4OhHFauHI+6des6pQ0xGSsIgpCBIiIiYDabARSyPaOBQlEU9+7dy8huiUAv\nCILgLEajEQUKFINS+SOAGAC7YDZvR8WKFTO0XyLQC4IgONGWLStRsuRmKJUyPD0/xooVC5A/f4rT\n0ziVGKMXBEFIAxaLJU1W24gxekEQhDdERi+ptPfm9EQQBEFIEyLQC4IgZHIi0AuCIGRyItALgiBk\nciLQC4IgZHIi0AuCIGRyItALgiBkciLQC4IgZHIi0AuCIGRyItALgiBkciLQC4IgZHIi0AuCIGRy\nItALgiBkciLQC4IgZHIi0AuCIGRyItALgiBkco4GeiWA6QD2AdgB4OV9spoCOAzgEIDPHe6dIAiC\nkGoqB8sFAdAAqAzAH8B423PPTQBQGkA4gLMAFgN46ng3BUEQBEc5GuirANhs+/0ggHIvvR4LwB2A\nBda9DcXmsIIgCBnE0UDvCiDU7rEZ1uEci+3xeABHYb2i//OlYwVBEIR05OgYfSgA40v1PA/yuQD0\nBJAbQB4AvgCaO9iOIAiCkEqOXtHvBfAhgOUAKgI4afeaDtYr/GhYg/89WIdxXjFs2LC43wMCAhAQ\nEOBgdwRBEDKn4OBgBAcHp6oORSrK/QLgPdvjjgDKAjAAmAWgD4CPAUQBuAigKwDTS3WQFEP3giAI\nKaFQKIAUxm5HA70ziEAvCEK6uXz5MtatWwe1Wo2WLVvCy8vLaXUvW7YcP/44BWazGV991RldunRy\nWt0vE4FeEAQhAf/88w9q1AhEbGxTKJVhMBr34fjx/fDz80t13evXr0erVl8gImI6AA1kuTumTRuM\nDh3ap77jCRCBXhAEIQHVqzfE7t1BsI4iAypVP3zxhQKTJ49Ldd2NGrXBhg31YB3BBoA18Pf/BQcO\nbEl13QlxJNCLFAiCIGR6Dx48AlA07rHJVAR37jx0St2SpAXwxO6Zp9BqNU6p21lEoBcEIdNr3Lge\nJGkogNsAQiDLExAUVM8pdX/7bW/I8kgAPwEYD1nuj2HD+jqlbmcRQzeCIGR6sbGx6NGjH+bPnwOL\nRYMqVcphxYrFTpuQPX78OKZOnQ2TyYzPPmuPSpUqOaXehIgxekEQhESMH/8zhgyZioiI3lCrz8LL\nayvOnDkMDw+PjO5aiohALwiCkAhXV1+Ehe0EUAQAIEktMGlSXXTr1i1jO5ZCYjJWEAQhEbGxUQCy\nxD22WLIgMjIy4zqUjkSgFwThndC8eStIUgdY8y0uhEr1Jxo2bJjBvUofYuhGEIR3QnR0NPr3/w7r\n1/8FLy8vTJ06Ev7+/hndrRQTY/SCIAiZnBijFwQh1U6ePIkOHb5A69ad8ffff2d0dwQnEIFeEIQ4\np06dQuXKtbFgQW4sXVoWjRt/grVr12Z0t0ASEyZMRsWK9dGgQUucOHEio7v0VhFDN4IgxOncuQfm\nzcsO8n+2Z1ajTJnJOHp0e4b2a/DgEZgwYS0iIoYDuAKDYQSOHduHAgUKZGi/MoIYuhEEIVViYmJB\nGuyeMSAmJsapbcTGxsJkenl7iqT98stsRET8DqAhgJ6IimqDZcuWO7VfmZkI9IIgxOnSpS0kaSSs\nm8dtgSz3RPfuzkm3Gxsbi3btukKSDNDp9OjWrTfMZnOyyiqVStjvXaRQmODiIsJXcol3ShCEODVq\n1MCff85DhQoz8f77IzF+fD98/nlXp9Q9fPhPWLXqGszmBzCb72DRomOYMGFyssr2798LstwKwGIo\nlT9Ckv5EmzZtnNKvd4EYoxcEIV1UqFAPhw/3AfCB7ZnlqF17MbZtW/nasiQxb958LF26AV5ebhg2\nbCAKFiyYpv19UzkyRu/o5uCCIAgpkitXVvzzzxGYzdZAr1IdQa5cWZNVVqFQoFOnDujUqUMa9jDz\nElf0giCki6tXr6J8+eqIiioLwASD4Sz++WePU7bze5eIO2MFQXijPXz4EJs3b4ZCoUCDBg3g7u6e\n4HEmkwkuLi7Pg5pgRwR6QRDeag8fPkRQUDvs27cNOp0BkyaNQ9eunTO6W28UsY5eEN5RFosFU6f+\nisaN26JXr364f/9+RnfJIW3adMHBg/lgsYQjImIfvvpqKPbs2ZOmbZ49exZlywbAyys36tQJwu3b\nt9O0vYwgrugFIRPo1as/5s7djYiIHlCrj8DXdwvOnj0Co9GY0V1LEb3eExER5wF4AwBcXAbi++/d\n8e2336ZJe48fP0aBAiXx+PF3IOtBpZqFAgW24syZQ7a1+28ecUUvCO8gs9mM6dOnIiJiI4D2iI2d\njCdP8mHjxo2JlgkLC8O5c+cQFhaWfh1NBk9PXwDHbI8s0GqPw8fHJ83aO3LkCMzmgiA/B5APJtNI\nXLt2Ezdv3kyzNjOCCPSC8JazWCywfjvW2T0rJ5pmYN269ciaNQ/8/Zsga9bcWL16Tbr0Mznmzp0M\nWW4LWe4EgyEARYtGoV27dmnWnsFggMl0G0Cs7ZknMJvDodfr06zNjCCGbgQhE2jevD02bnyMyMh+\nUCoPw939Z/z77zF4e3vHO+7x48fIkaMgIiI2APAHcASyHIhr1/6Fl5dXhvT9ZSEhIdi5cyfc3d3R\npEkTaDSaNGvLYrGgbt0gHDgQgYiIWtDrl6FDh9qYOnV8mrWZWmLVjSC8o6Kjo/G//w3H1q27kSOH\nHyZPHplgZsejR4+iVq3OCA09Hvecq2s5/PXXtLdytyVniI2Nxdy5c3HhwmWUL18GLVu2fKOXdYpA\nLwhCku7fv49cuQojKmo/gMIALkKnq4ArV84ia9bk3aUqZCwxGSsIQpK8vb0xbdpESFIVuLnVgCRV\nxOTJ40SQz+QcvaJXAvgFwHsAogF0AXDJ7vXyAMbb6r8JoD2Al5Naiyt6Qcgg169fx8WLF5E/f37k\nzp07o7sjpEB6Dt00A9AIQCdYZ3S+BRBkV+c/AD4CcBlAVwC7AJx/qQ4R6AVBEFIoPYduqgDYbPv9\nIIBydq8VAvAQQF8AwQDc8WqQFwRBENKJo4HeFUCo3WOzXV1eACoDmAKgDoDaAGo62kFBEAQhdRzN\nRx8KwP7eaiUAi+33hwAu4sVV/GZYr/h3vFzJsGHD4n4PCAhAQECAg90RBEHInIKDgxEcHJyqOlIz\nRv8hgI4AKgIYDOuuvQCgAfAvgLqwTtD+CWA2gE0v1SHG6AVBEFIoPSdjFXix6gawBvyyAAwAZsE6\nVDPKdtxeAH0SqEMEekEQhBQSN0wJgvDG2bZtG8aNmwGS6Nu3K+rXr5/RXXqriUAvCEKaeJ7vfsuW\n3ciVKyuGDfsWvr6+ry23bds2NG7cDpGR1i/4kvQNVq6ch8DAwLTvdCYlAr0gvANCQ0OxYMECPH36\nFPXq1UP58uXTvM2ePfth3ry9iIjoBZXqCHx81uPs2SNwc3NLstwHH7TE5s0fwDq6CwALUbv2Kmzb\ntjLN+5xZiRQIgpDJhYaGolSpyhgwYCeGDg1FjRqNsHr16jRt02w2Y8aMabaMl21hMk1EaGghbNiw\n4bVlrRdz9jFJAXGBl/5EoBeEt8hvv/2G27eLISpqOczm0YiMXIqePb9J0zYtFgssFgvswwWpg9ls\nRmhoKC5fvoyYmJcznFj169cNkvQNgPkAFkKr7Y/+/bulaX+FVzm6jl4QhHR26dIl/PXXVkRHA0A4\nAD2AAggLe5JmbV6+fBl16jQBoAbgB6AnlEpvaLUHcetWJfj45IBK5QFJIrZuXYtSpUrFK1+1alV4\nebniv/++B6kHYMb1647t3hQWFoYxYybgwoXrCAjwR7duXd7Y7f6EFygIQvLs2bOHer0XtdoOBKoR\nKETgPHW65mzdulOSZTds2MC2bbuye/evePny5RS1my1bfgI5bD/tCLizWrW6XLduHWU5K4FLBEhg\nEbNmzUeLxRKv/IIFC6jX1yNgsR13mgaDV4rPPyoqisWKladW+zGBGZRlf3bp0jPF9WQGAFI89iU+\nDgUhA9y+fRshISGIiorC0KE/oFSpGqhbtylOnTqV4PGffdYf4eHTEB09D8BOAMWhVpdBkyZGzJ49\nOdF2Fiz4HS1afI5Fi97H9OkySpeujOvXryerj7t27cKtW/cBLLS1+QhAdtSvXwOPHz+GUlkTQD7b\n0R/jwYM7r+xBGxoaCrM5D16M0+dBVFRYisfpd+zYgRs3lIiO/h1AN0REbMZvv81BeHh4iup5V4lA\nLwjpiCS6du2FvHmLo0yZ+vDxyYexY//CiRND8PfftVGlSh1cu3btlXL3798D8L7tkQJABfTo8QWW\nLJmb5P6mQ4aMQUREZwCnYbGoEBbWGLNnz01WXzdt+gtAbwABsAb0SQBu4NChf5A/f36QBwA8Hzba\nDUmSYTQa49VRp04dKJUrAawHcA1a7ReoW/fDFO/gFB0dDYXCHS8+MAxQKlWJzg0I8YlALwjpaPHi\nxVi8eB+io68iPPwywsJCERm5BEBtkD0RE9MYa9a8ull3nTq1oNUOA/AMQAhkeSbq1n19rsCHD+8B\n+ANAMQC3YLGsRWho2GtKWXl4uCL+NhPXAGiQJYsvKleujM6dW0CSisPNrTb0+mZYsWLRKwG8cOHC\nWLduKfLnHwwPjypo1AhYsmROstq3V716dWg0Z6FUjgNwCFptZ/j7V4GHh0eK6xLSV0YPdQlCuuvf\nfyCBH2zj1SSQhcBF2+9RdHEpzfz5y7BLl568f/9+XLmwsDA2atSSKpWWkuTOiRMnv1L3+vXrWb16\nI1at2pBr1qyhxWKhSiUTuGLXXk0OGzYsWX198OABjUZfAq0I/I+AJzUaD27bto2PHz/mxo0bOXfu\nXG7YsIG3bt1y2nuUmIsXL7J27SbMn78027f/jKGhoWne5psIDozRZ6SMfr8EId3NmjWLslydQJQt\n8DamUlmAwHwCJQnUJrCKanVP5stXkhEREfHKvzzZ+dymTZsoy34E/iCwhJKUjatXr6ZKpSXwNC7Q\na7XtOH369GT39/79+6xbN5BGow9z5izCpUuX8fz588ySJSddXWvRYHiPFSrUZGRkZKreFyH5IAK9\nIDjf0aNH2aHD52zXriv37NmTqrpiY2PZoEFz6vV56epakVmy5OQPP/zIevWaUanU230AWGg0VuRf\nf/2VrHoDA1sQ+M3uyn0Ra9ZswubN21OSmhA4QmAuDQZvXr16NVXnUKVKIBWKSbZ2TNTpGnPs2HGp\nqnPPnj2cOXMmd+7cmap63gVwINCLdfSCkIQjR46gRo0PEBExAIAGK1c2w7p1i1GrVi2H6lOpVFi/\nfhmOHz+OsLAwlC5dGkajEd263UeOHAURf27xxXi3yWTC6dOnAQAlSpSAShX/T9c6Nh5r94wJSqUS\nCxfOQJ8+32LLls7w8fHGtGmbUr1H7NWrV0HWtj1yQVRUTYSEXEqyjL2IiAjs3bsXCoUCVapUwahR\nEzBu3CwAtaBQjMLnn7fGuHE/pqqPwpsjoz8YBeG1mjf/lMAkuyvlBaxevZHT27FYLKxfvyklKYjA\nOqrVXzFPnmIMDw/n06dP+f77lWkwFKTBUJClSlXh06dP45XfsWMHJcmbwHQCMylJvsn+NpBSQUFt\nqVb3IGAm8ISyXJ7z5s1LssylS5c4ePBQ9ur1FbNlK0ijsRKNRn/mzl2MWq0bgTu29/cBJcmbly5d\nSpO+ZwYQQzeC4FyNGrUhMMcu0K+mv3+9NGkrMjKSAwZ8x0qVAvnpp5/z3r17JMkvvuhju1HKTMBM\nrbYDe/To+0r54OBgNm78MT/8sA3//vvvNOkjaZ2kLV26KnW6LFSr9ezatVeicwckef78eRqNPnRx\n6UugAoEeccNTanV3qtVZ7d5f0s2tPPft25dm/X/bQQR6QXCu9evXU5azE1hDYBNlOR/nz1+Yrn2o\nXPkDAmvtguEaVqnSIMFjo6KiGBMTk+Z9slgsvHXrFh8/fvzaY9u3/4wKxfe2vtcnsN7uXFbSxcWL\nwGLbB9kqurr68smTJ2l+Dm8riDtjBeFVq1evhp9fAej1nmjSpA1CQ0NfX8imYcOGWLDgZ5QqNQHv\nvTcSU6cOQfv27dKwt68qU6Y4tNplAMwAzNDplqFs2eLxjomKikKTJm2g17tClo3o2bOfLRFZ2lAo\nFPDz84O7u/trjw0NDQeZzfaoLIAZAKIBREGS5qJ9+4+QPfsQKBQa+Pr2xZYtq1+b/lh4e2T0B6Pw\nDjh69CglyYfATgJ3qdW2Z6NGrZxW/927d1mz5oeUZU/myVMiRatGLBYLQ0NDkxz2IK1r6MuVq0G9\nPg/1+twsXz6AYWFh8Y7p1WuAbXw/gsBDynIFTp36q0Pn5GzLl6+gLOclsJfAfiqVvnRxMVCjcWOD\nBs0ZFRVFkoyOjs7gnr4dIIZuBCG+0aNHU6X6ym6o4AElyc1p9ZctW50qVR8CdwmsoV7vlazli6dP\nn2aOHIWoUkk0GLJw3bp1SR5vNpt55swZnj17lmaz+ZXXixWrbPswe36e8xgU1M7h83K2GTNmMWfO\n4syWrTCHDx/JO3fuxM1BCCkDMXQjCPF5eHhAo7mAF38bIdDrXz/ckBzh4eE4fvwgTKZxAHwANEZ0\ndCXs2rUryXJmsxl16jTGf/99A5MpHM+erUerVh2TTDamVCpRrFgxFC1aNMHUvDlz+kGpPBT3WKM5\niDx5sr1yXEbp1q0Lrl8/jZs3/8WQIdZtCL29vTO6W+8MEeiFTK1t27bIlesOZLkJVKoBkKSmmDZt\nbKrrvXDhAlauXAmz2QLghu1ZM0ymSzh//nySZe/evYunT5/Bur2eAkBFqFQVcPz4cYf7M2XKT3B3\nnwSDIQgGQ134+e3C4MEDHa5PyFzEDVNCpibLMo4e3YWFCxfi0aNHqF17LSpUqBD3ekhICC5evIjC\nhQsjf/78yapz9eo1+PjjLnBxqQVroC4F6yYgrlAoYl47kejh4QGLJRLAeQCFAYTBZDqDbNkcvwIv\nWLAg/v33GP7++2+oVCoEBgbCYDA4XJ8gOEtGD3UJ77gJEyZTknzo5laPkuTNmTPnvLaMxWKh0ehF\n4CCBKbRuAHLUNtGYg2q1F3fs2PHaembPnkdJ8qXB0Jp6fQF269Y72f02mUy8e/cuY2Njk10mMc+e\nPeP27du5d+9ep9QnpD2IyVhBSJ7r169TkrIQuGabvAyhTuceL2OkvUWL/mCJElVYrFglAi62Nd81\nCWyONwFapEiFZPfh1KlTXLhwIXft2pXsMvv376eHRzZqtZ7U6z25cePGZJd92Y0bN5g9e0G6ulak\nwVCc5crVYHh4uMP1CekDItALQvLs2bOHbm7+8e7IdHUtzuPHj79y7KpVqyjLOQlsIvAXFYpsBFoT\n+DDeXbMKxff85JNuadbniIgIurv72W7eIgHr9oK3b992qL6GDVvSxWWorS4zdbrmHDr0+0SPf/To\nEc+ePftKRk0hfUGsuhGE5ClUqBBiYy8COGh7JhgWy13kzZv3lWNnzPgDERE/AggEUBfkFOh0uwBs\ngnUHpkEA+oL8CfXrV49X9smTJ1i7di02b96MqKioVPX52rVrMJn0ABrbnqkClaoozp4961B9ISGX\nYTYH2h4pERVVH2fPJpycbMaM2fDzywN//ybw88uH/fv3O9SmkDFEoBfeSd7e3liy5Dfo9Q0gyzlh\nMLTA6tVL4Orq+sqxOp0GL7bMA4AnqFKlAt57rzKAUbBeYBkA9MOGDTvijrpy5QoKFXof7dpNRcuW\nI1C6dNUU3ZX7sqxZsyI29h5e7Pp0DzEx55EzZ06H6itfvhQ0mnkALAAiIMt/oGLFUq8cd+7cOfTp\nMwjR0UcRFhaCp09nomHD5jCbzY6eivAOyehvQILAiIgIXr58OcmNMw4dOkRZ9iLwE4GxlCQv7ty5\nk2XK1CSw2m74ZxzbtesaVy4wsDmVyp/iEnhptZ9y4MDvUtXfX3+dSVn2patrEGU5OwcPTnyo5XUe\nPXrEMmWqUZKyUqv14EcffZLghOyKFSvo6tok3jCXTpeFd+7cSc2pcNOmTQwMbMFGjVqnax766Oho\n9uzZn9mzF2GRIhXSLMtnWoEYoxeEtPHPP/+wU6fu7NjxCx48eJAk+eeff9oSns0lMIWy7MVDhw7F\nlSlUqDyBfXYBcjY/+qg9zWYzR44cy3Ll6vCDD1rw1KlTKerL2bNnuXz5ch47dizV52U2m3n16tUk\nx/mPHTtGWc5ml0p4L/V6z1St0lm3bh0lyY/APAIzKEneKZqUTo1u3XpTkuoROElgDWXZ2ynvZXpB\nOgZ6JYDpAPYB2AEgsQXIMwH8lMhrGf1+CUKqrV+/ng0atGLTpu144MCBeK916dKLOl0bAjEEnlKW\nK3PKlGns0+cbynJFAhupUEyi0eiT6l2fEhMZGemUHDJDhvxgW4panXq9V6pW+5Bk1aoNCCyx+xD8\nJd1SNri7ZyNwOa5tpXIghw8fkS5tOwPSMdA3AzDX9rs/gNUJHPMZrB8EIxOpI6PfL+Edcv78eZYt\nW4Ourlnp7187XTa2CAsLY82aDalWG6hSSezY8QuazWYaDF4ErsYFGo2mG8ePHx9XLiIiggMGDGL1\n6h+ye/c+DqXsjYqKYrNmbenioqGLi4adO/egyWRK1fmEhITw77//dniVjz1r6uXldoF+Bps0aZvq\nepMja9YCBPbHta3Vtue4canbCjE9IR0D/XgALe0e//fS65UBzAHwKcQVvZDBwsPDmSVLTgITCNwg\nMJJZsuSKy5qY1h48eBAv26Srqy+BC3aBpiMnTZpE0npDVo0aDajTfURgJTWaTixevEKKc8z36/c/\nSlJDAuG07gJVhePGTXLqeSXXkSNH2KpVRwYFtePmzZtJkitXrqQs57DloZ9PWfbl9u3b06U/v/22\nwDbkNooq1Wf09c2T6P0TbyKkY6CfBetas+eu4cUKHj8AmwFIADpABHohg+3bt48KRcF4k4lATm7Z\nsiVN242IiGCTJm3o4qKhRiNzyJDvabFYOHToD5TlkgT+oFI5lO7ufrx16xZJ8sqVK5QkX9twj3US\n12AokeIdl0qXDiCw1e58FzEwsEVanGaS/vnnH9tE9njbWLwfV69eTdJ6f0KNGh+yVq2gdJ8QXbx4\nMf39q7B69YC3bkNypOPm4KEAjHaPlbCu0QKA5gC8AGwEkBWADOAcgAUvVzJs2LC43wMCAhAQEOBg\ndwThVbGxsYiKikJsbCzIBwAiYP3nGAbgKa5evZpoWbPZjB9+GI1Vq7bAxycLxo0bhvfeey9F7ffs\n2R9btkTDbH4Es/kxxo2rj6JFC2Do0P8he/as+PPPlfDx8cTw4Xvh5+cHALD+HStgvzG4QqG0PZ98\nuXJlw8mTB2A21wEAqFT7MiSb5cSJ0xER8TWAvgCAyEgfjBjxM5o0aYKgoCAEBQWle5+uXLmCL77o\ni2fPmsJiMeCDDz7C1q1rULly5XTvS3IEBwcjODg4Q9puBmCe7feKADYkcpwYuhEyxIgRP1Gl0lGl\nklimTHUCegLlCQwlUJYqlTf3798fd/zt27c5Y8YMzpo1i/fv32ePHn0py9UIbKNCMZVGozevXLmS\nrLYtFgt79OhLwI3AEbur6ils3/6z15atXLmubRJ3AzWaL1i4cJkUT6heuXKFXl45KUnVqFB4ElAy\nS5Yc3L17d4rqSa02bToTmGz3Hmzi++9XT9c+vOyzz3pTqRxk16f5rFw5MEP7lBJIx6EbBYBfAey1\n/RQC0AZA15eO+xRiMlZIZ9Z9XgsQ+I+AmWp1bxYtWp4ajRtVqvep1Wbnxx93jtvZKSQkhO7ufpTl\ntpTlVrYA6W4r/3wcvWvcOPrrTJ8+k7JcntZcODPtJl07ctCgIa8t/+zZM/bq1Z/+/vXYqVMPPnz4\n0KH34ebNmzQYfAgsImAhsJFGow8fPHjgUH2O2Llzp22Hrz8IrKUs5+ecOfNSVefrduR6nRYtOhCY\nYRfog1m8eOVU1ZmeINbRCwL57beDCAyz+0O+Rnf3bDxy5AhnzpzJLVu2xAsWjRu3oVI5Ku54F5dB\nVKmM8SZMdbp2nDp1arLatwaSWQSOEfAh0IoKRTXmzVs8WZtpP7dv3z6OHj2av/32m0NLJM+cOUOj\nsVC8uQk3t8oMDg5OcV2psWXLFlap8gHLlavNefPmxz0fGxvLK1euJHtV0dOnT1mvXlO6uGhoMGTh\n9OmzHOrPn3+utG1teJDAWcqyP7//fpRDdWUEiEAvCOSvv/5KSQqkNcMkCSxmkSLlEz2+fPk6BDba\nBcSlzJfvPcpyCQLz6eLyLbNkycG7d+8mq/2+fQdSre5ou4q+QYWiJUuUKPvKPq9JmT17LmU5G1Wq\nvtTra9Lu0iLHAAAgAElEQVTfv1aKV97cu3ePWq0rgZu283pESfLlv//+G++4CxcucMaMGVy8eHGS\ndwg7U0hICLNnL0hZzk6NxsCRI8e+tkxQUFtqtZ8SeEbgDGU5p8MrdaZPn0k/v0L08srDgQMHJ7g9\n45sKItALgnUNeYUKNWkwlKPR2JQGg3fc3awJGTHiJ8pydQL3CNyiLJfnpElTOGfOPDZq1IZduvTk\n9evXX9vus2fPGBj4EV1cNARUVKvz02BoRE/P7Dx//nyy+2+xWCjL7gTO8HlmSYOhMlesWJHsOp4b\nOXIsZTknZbkD9fr8/PLLgfFeDw4Opix7UZI60GCoxeLFK6RLquKiRctToXg+dv8fZTn3a++MdXPz\nI3CdwCMCXxF4j1Wr1nrnNhWHCPSCYBUbG8uNGzdyyZIl/O+//5I81mQysVu33tRoZGo0en755QCH\nrvCCgloTaEIgisAdqtWF2bNnzyTXaFssllcCVWxsLJVKld0SS1KWO3DmzJkp7hNJHjhwgLNmzUpw\nyKZAgdJ8ka/HQkkK4uTJkx1qJ7ksFguVShfb+/R8DqQ7f/755yTL5c37HoGVBEoR6EZgBdXqQDZo\n0DxN+/umgQj0guC46dNnUpLcqFSqWKlS3UQDtMlk4tGjR3no0KG4IG0Nzp607jb1YpVNixbtE21v\nxozZ1OlcqVSqWKFCrXhDQxUr1qFa3dt29bqNsuz1ypCLM3h4ZKf9XbrAUH777SCnt/MyP78CBNba\n2gynXl+S69atS7LMtm3bbENRJW3DYiQQRa3WPdUJ1t4mEIFeEByze/du292S5wjEUK3uzVq1GjMm\nJobbtm1j3brNWLNmEy5btozlytWgXl+QRmNxFi5chg8ePODVq1dtgX563NUx0Jpt2iScv2XPnj22\n9v4lEEuVqg9r1GgY9/q9e/cYENCIWq2Bfn4FuGnTpjQ576ZN21Kj6UggksB5ynJubt26NU3asrdn\nzx4aDN50c6tNvT5PvFVQSZk3bx51utJ2gT6GWq2nU9IyvC0gAr0gOOann36ii0t/uyvbBwS0BJRU\nKFxpzbK4iCqVB9XqFgRMBCxUq3vyk0+6MSwsjCqVRMCbQAsCNalQuCY67jx69GiqVH3t2ntErdaY\noj4fPHiQs2bNSlXqgCdPnrBOnSZ0cVFTktw4ZcovDteVUnfv3uXmzZt59OjRZC+ZjIiIYJ48xalW\n9yGwkTpdc9aq1SjVSy7fJhCBXhAcM3fuXOr1te1W6mwhkJfAx3ZX6SRQjdY14c8fb2POnEWYO3dJ\nurnlpErlSo2mIrVaX3bp0jPR9ubPn0+9PsCuvb+YLVvBZPd3zJiJlOXstknWguzevW+qzt9kMr01\nwfLu3bv85JNurFChLr/6auA7t7UhHAj0itcfkmZsfRaEjBcTE4Nq1QJx9mw0YmPzIjp6LYA/AcwH\nUBXA57YjW0GpDIPFsgaAC1SqeiDPwWxeCkCGTtcRQUEl8cUX3VC9evVEWrOmZ6hRowFOnQoDWQgW\ny0YMGNADZrMF+fLlRfv27aFSJZyh5MmTJ/D1zYWYmDMAcgJ4ClkujoMHN6NEiRJOfFeEN5FCoQBS\nGLvFVoJCopYvX4HcuUvA2zsPPv/8K8TExDit7tjYWHz//SgEBrZEnz4D8fTpU6fV7QiNRoM9e7Zg\n3rw+6NXLD5LkBaAagG4AhgCYDWABJGk3ihV7CFnOC70+H/T6CzCbhwHIAuu+q6Nw7tyNJIM8AKjV\nauzatQkLF36DyZMD0Lp1C4wf/yd+/NEFvXotQMOGLWCxWEASv/wyA5Urf4AGDVri2LFjePDgAdTq\nLLAGeQBwQ2xsDixfvjzR9s6cOYMSJSrBYPBChQq1cOXKldS/aYKQDBn8BUhIyq5duyhJWQnsIBBC\nSaqX6uEBe02atKEs1yfwBzWaTixSpGyy0gZHRUXx5MmTvHbtmtP6QlqHAypWrEOlUkW93pPlylWn\nwVCSstyJGo0HixevxDp1mnHDhg20WCwMCQnhuXPn2K5dVwIlCOQkUJhAPlaoUDtFbT969IgajcE2\nL7CaQFUqFN788su+HDlyjC3T5WoCU6jXe/HkyZP09s5lmzew2IaZ3ClJeTl27MRX6n/69Ck9PbNT\noZhJ4A6VyjHMkaPQO7f+PLOAGKMXnKV//4EERtiNRZ+jr29+p9R97949ajRuBCLiVqgYjWVfO6l4\n+fJl5shRiEZjYep0XuzY8QunjStXrRpIlaqPbe36CUqSLydOnMgZM2bw6NGjiZbr3/9rAgEEom1B\n90tWq5ayBFnXr1+3fahuJeBnW3a4g1ptERqNvrSmUni+G9IAfvfdEJ48eZIeHjkIuBDITiCYwDka\njd6v1L9r1y66uvrb/b8kDYb8PHfuXNwx586dY+XK9Zk9e1G2aPFpilI1COkLDgR6MXQjJMjd3RVq\n9XW7Z67BYHB1St2rV6+GyWTBiyzZCigUWpjN5iTLtWnTFbdudUZY2L+IirqMZcsOYunSpU7p04ED\nwTCZRgBQA3gPZnMLkES3bt1QpkyZRMvduPEAwMcANLAOm36MW7fuJavN//77D/36fYPBg0fC09MV\nwP9sPx8CCEB09FRERloA2L8vZiiVSpQsWRJDh/aHRtMRwAQA7QFUQHh4NEJDQ+O14+7uDpPpFoBI\n2zNPEBv7EG5ubgCAhw8fokqVOti/vyFu3lyKNWu0CAz8KMWpkQUhIRn9wSgk4d69e/T1zUuNpiMV\niu8oST6vvaElOebNm09Jyk2gDIHmBLZRpfqauXIV4bNnz5Is6+qaldYdopx/c4+XVy4Cu/g85YBe\nX52LFi16bbnRo8dSkj6Iu6JXq/vxvfcqs2TJaixVqhpXrVpFkrxz5w6/+moAW7bsyN9/X8SbN2/S\n0zM7XVz6EphEnS4bXV2zM34ythXMkaM4ZbkwgSVUKMbQYPDmhQsXSJIXL16kTudGwIvWrfEeUqls\ny0aNWsXro8ViYYsW7anXV6BCMYh6fUn26NEv7vXVq1fTaKxn166JGo1bklkuT5w4wTlz5vCvv/56\na1brZBYQQzeCM92/f5+jR4/moEGDk8wVkxKFCpUn8Detian6ECjMAgVKJeuGlzJlatjlR4mgXl+R\nCxcudEq/1q5dS1n2pix3psFQmZUq1WFMTAzPnDnDmjU/ZKFC5dmjR79Xkn5FRkYyT54SBLJQochB\ng8GXanV+AhtsY+h6Dh48hFmz5qVK1ZvADMpyMdauHUi1+jO74LqXPj55qdd7UaEYRmACZdmX69ev\n57x581mnTjM2b96ep0+fjtd+9+7dCXxpV889SpL7K+dnNpu5aNEiDh06jCtXrowXnNetW0eFoojd\nUs9HVCq1iSZhs27F50u9vj31+mJs1aqjCPbpCCLQC2+6AgXKEthpF5h+YPfuXyWr7Pnz5+ntnZuu\nrmUpy9nZokV7p2Ud/OWXGfTxyU+jMRs/+qg1o6OjeevWLbq5ZbV9uOylTteYTZvG38B6+PCRlCR/\nAmtsgd2d9htPA8OoUGgoSR/ZPXeJarWBCsW3ds/9Sy+vPDx79iw/+6w3P/3082RtcTdjxgzKckO+\nuFN0F31986Xo3JctW0al0pdAMwITCZShQiExNjb2lWNjYmKo1RoInLX7wC381m3H9zaDCPSCMz14\n8IDLli3j6tWrU3RTyuzZc1m2bC1Wrhz4yr6s06ZNpywXIrCKwCzKsleSk50ve/bsGffv388zZ868\n9iryyZMnPHz4MG/evJnkccuXr6As5yNwgMApynI5jho1nr/99hv1+lZxwxnAzwQkjh8/Pm7FSpEi\n/gR22wXsnLaJ0eeP/0elUk+ttoPdcw+oVsuUZW8CSwkcoCxXZZ8+3yT7fXguPDycRYuWoyx/QLX6\nS8qyT9xwUXItWrSIen0TAqMJ9CAwh0qlJsGUxQ8fPqRG4xpvYtfV9SMuWbIkxX0XHAMR6AVnCQkJ\noadndhqNjWg0VmfBgu8na4OImTNnU5YL2oYuFlOSfOKlAbBYLJw79zdWqhTIOnWacu/evWnS/x07\ndtBg8KaraynqdB4cOzbx3aGaNv2EwGy74LWN779fnX/88QcNhg9sz9UkkItAQ+p0tVizZiOazWaW\nL1+b1q3y1tK6OqY2rWkQ5hEYRcCTWq0rDQZvWnc12ktJqsf27T/j9u3bWbp0DebPX4b/+98wmkwm\nh841IiKCs2fP5tixY/nPP/+kuPzt27dt31ymEzhGrbYt69RpkuCxFouFuXMXo0Lxs+1bxAHKshcv\nXrzoUN+FlIMI9IKz1KkTRKVyHJ8vf9RoOnLgwO9eW65EiSoE/rILmhPZtm0XHj9+nCdOnHA4mKVE\nbGwsXV19aF2uSALXKctZefLkyQSP79SpO4GvCZynNXXuXFar1oChoaHMlasIlcpyBHIQGE/gEwIl\nqNXmYvPmbVirVh1a96MNIOBDjcaTnTp1oVLpQZXKl5Lkxo0bN/LYsWOsVq0BCxeuwL59v33j1rCf\nPHmSlSrVY65cJdi+/WdJbpISEhLCAgXep1KpptHoxTVr1qRjTwWIQC84S8GC5V4aa57Jli07vrZc\nqVI1bOPVz8v9SE/PXDQYCtBgKMDSpasyNDQ0Tft+584d6nRZXhpeCOLy5csTPP7bbwcT0NG6Hj0L\ndTp37tu3j6R1QlqplAhcifvQA+rTut79cwIygcN8sYNTLh48eJCPHz/mqVOnGBoayocPH3Lv3r28\nfPmyQ+djsVg4duxEenvnoadnTg4aNOyN2BEpMjJSTMJmAIhALzhL1669qFI1o3XZ4CPKcnlOn/76\njS9WrlxJWc5mG6YYS5XKlWr1h7Su6DBTq23Pnj37JVnHxo0b2b37Vxw+/Hs+evQoxX1/cUW/zRaA\nb1CW/RK8oj98+LCtv9dsx/5OH588ca+bTCbbjlHP7D44WhPIQmsed/eXPlCacunSpXHlg4ODbel4\nK1Cn8+KgQSNSfD7z5y+kLBchcJzAOcpyuQTvgBXeDRCBXnCWjh0/p1KZjdZUvWoajX7J3vN006ZN\nbNr0E378cWe+/34VAt8TaEfgUwJjWKVKg0TLTpnyC2U5D4Ex1Gg6MmfOwsnePNre9u3bbWP0panT\neXLMmIQDozVr5Sd2wdpCpVIdb/L5ww9bUattQeAkgQW2q/jNtg+v3HyRzfIsZdkn7o5Ti8VCDw8/\nWlMUkMBdynIuHjhwIEXn0qBBKwIL7fq4ieXKpSzNgpB5QAR6wRmioqLo4qIl8MR2JRtJozGAI0eO\nZK1aTejvX4+zZs1J1tf2Dz9sRsDTdoU/mYArAwLqJnqLvaurL4HTcUFNkppyxowZDp3H48ePefDg\nwSS3Ety5cyf1+vwEHtva3EoPj2zxzi08PJwdOnzB7NmLsHTp6gwICKRO9xGBUwRGEpCp02WlVmvk\nb78tiCsXGhpqy1Fvn3qgDRcsWJBQVxLVrl1XKpX26Sh+Yb16zVL+hgiZAkSgF5whPDzcFuh7Eggk\n0IeSVJEajQeBOQTWUK8vkqxNKvz96xFYYhekfqZKlYPu7n48cuTIK8drtUYC9+OO12q7c+LEtB2m\n6NVrACXJj25uNWgweL825054eDg7duzObNmK8P33q3LXrl28fv36K5tqWywWennlpHUpKWndBDt7\nguedlAsXLtDV1Zcq1Rd0cfmKer2XQ6trhMwBItALzmAymWy347clsI5AW9suS/ZXlbuYP3+Z19Zl\nDfSr7MrNsA3jLGaePCVeOb51647U6ZrYruqXUKVyo7d3XubJ8x4XL067tdpnzpzhtm3bGBISwgYN\nWlCWPejnVzDVaR8OHDhAd3c/Go3FqNW6c9So8Q7Vc+3aNf7000/8/vsfGBISkqo+CWkjJiaGAwZ8\nx0KFyrNKlUAePnw4TdqBCPSCM5w4ccI2nPH8lngzVarstlUmL9aaFyxY7rV1LV68hLKc2xbsFxPI\napskDaeLi+aV4yMjI9mlSy9my1aY3t55qdVWtI2Nb6ckZeO2bdtSdC7R0dE8fvw4Q0JCkjXUVLdu\nEDWazgTuEdhOWfZOdFlmcj179ownTpx4pzawfhd16dKTslybwF4Cc6jXp839BRCBXnCGkydPUq/P\nR+vdoNa7QiUpNyXJncA4AvMpy7k5b978ZNW3ePESFiniT6XSJ27iUqGYzsKFy8Y77tq1a1yxYgX3\n7NlDi8XCvHnfJ3DI7sNlHLt2TXx7vpfdunWLefOWoMFQmLKcjY0atUzwtn57KpWWQGi8oaOff/45\n2W0K7y7r38fNuH87Gs3nnDBhgtPbgUhTLDhDsWLFUKRIDmi1HQCshlb7KYoXz4M9e/5Gy5an0aDB\nRvTu/TEuX76EOXPmwGQyJVlf69atcO7cAQwY8Dm02u4wGgvD23sUVq1aGHfMli1bULRoWXTqtAD1\n63dE69YdodfrAdyKO8bF5Rbc3AzJPo9OnXrjxo3GePbsHCIiLmH79geYPn1GkmX0eg8AF2yPCJXq\nAtzd3ZPdpvDuUqu1AF6kiFYqn0Kj0WRch94QTv+kE5wnLCyMvXsPYLVqjdi794B4Syu/+24E9fpC\nBAZTr6/JgICGyb7j9fbt2zxz5swru0l5ej7fPMOaKMtgKMkRI0ZQln0IDKNK1ZMeHtl4/fr1ZJ9D\njhzFCJyw+0YwmR07fpFkmQULfqcs+9HFZQBluQFLlqyYYM4XQXjZ6NHjbXmcplOl+oo+Prl5//59\np7cDB67oxebgQoqEh4fD3d0bJtNlAFkBmGAwlMK6dVMREBDgUJ0mkwkajRZkDAAXAIAsd8X48WVR\npkwZLF++Cnq9hC5dOiFHjhzJrrd+/Y/w998lYDYPBxANSWqEUaMao3fvXkmW279/P4KDg+Hl5YV2\n7dpBkqRktbd27Vp89lk/hIY+RK1a9fD77zPiNvcQ3g1LlizF6tV/IWtWTwwc2Bd+fn5Ob8ORzcEd\nDfRKAL8AeA9ANIAuAC7Zvd4GwJcATABOAeiOVz+FRKB/S+zevRtHjx5F7ty5UalSJeTJUwzR0Q/x\n/J+Pq2sDzJ/fDaGhz7B8+Ub4+Hhg8OAByJMnT7LbKFKkHEJC2oH8CsAlyHJ17Ny5Bnfu3MGqVZvg\n7e2BPn16wdfXN9l13rx5ExUr1satW7GwWJ5CklywdetKVKlSJWVvQDKcOHEClSrVRWTkcgBFodEM\nRO3az7BxY+IbdguCIxwJ9I5qBmCu7Xd/AKvtXpMAXASgsz3+A9a90V7m9K80gvONHTuRspyLWm1P\nGgxl2KxZOxYtWo4q1TcE/iMwkzqdkbVr16VOl4/AfCqV39HDI1uS6YHNZjP/+OMPDh8+nCtXrmRI\nSAhz5y5KrdaTGo2ev/wyg7/+OtN2l+xEqlQ96OubJ0Vfhc1mM/Pnf49KZTfbSohlNBp9krXJSUqN\nHz+eGk0vu2Gix9RqDYkef/nyZe7evTvJXZwyA4vFwl9/ncmaNZuwefP28fapFRyDdFx1Mx5AS7vH\n/9n9rgDgbfd4GYA6CdSR0e+XkITQ0FD++uuvthunTtqCVyT1+oJcu3Yta9RoSI3GQIXCQKWyMYGu\ntOZ/2W1brdI53oqD1atXs2nTT1ivXiNOmDCBDRu2sG1t9z/q9SXYs2c/WiwW3rt3L2783ssrN4Gj\nditg2nHSpMTTDb/s1q1b1Om8+GJTDtLV9YM0ybZoTaUQaNfWQXp65kjw2OHDf6JO50U3t4rU6724\ndetWp/fnTTFixE+U5ZIEllOhGE2j0YdXrlzJ6G691ZCOgX4WgEC7x9eQ8AqeXgA2JFJHRr9fQiKO\nHz9OtdqT1myOZWnN1HjcFigbcM2aNZw2bTpVqiwEOtldxS4hUJUAqVb34JgxY0iSc+f+RknKRaAE\ngZIE/G0fChF8nvVRq3XnrVu34vXDaPThi2RjpErVhz/99FOyzyMsLIxqtUzglq2OaOr1ReLlx3eW\niIgIFitWnrLc0DaRmzXBG7z++ecfynJ2ArdtfdpBo9E7XdI3Z4QsWXIROBP3/1Ct7s5Ro0ZldLfe\nanAg0KtSWsAmFIDR7rESgOWlx2MAFADwUWKVDBs2LO73gIAAhyfzBOcKDAxCbGwhADsBaADMA9AN\nwBBYLIdRvvxs9O5dHyZTbQDF7EoWgXU55K/QaJaiWbMDAIARIyYgMrI5gLOwfu7vBdAH1lE+APCA\nWp0FT58+jTd51bp1K/z+e2dERo4CcBEazUJ8+GFwss/DYDBg0KDvMHZsNURFNYNOtxfVqhVLkzF6\nSZJw5MhOLFy4EA8fPkStWqvh7+//ynEXLlyAi4s/rBPZABCAmJhYPH78GF5eXk7vV0azxiWF3WOl\n7TkhuYKDgxEcHJwhbTeD9a8fACri1av2WQCmIOkJg4z+YBQSoVSqCXxnd6V+k4BEP78CcXuD5spV\nwnbzVF7b0M5dKpV1mCVLXtauHcRjx47F1Zc9e1ECHQkMt9X3lMDzVMZ3qFSOZY4chV7ZjCMmJoZf\nfTWQefK8x1Klqjt8Jb5lyxaOHDmSixYtYkxMDHft2sU1a9ZkyJ2qJ0+epCxnpTXFMQlsoIeHX7rk\nl7979y7nzJnDOXPmxJvr2LVrFxs0aMW6dT/i2rVrEywbGhrKS5cupXjDlKFDf6Asv09gFRWK8TQa\nvR3Oyy9YIR2HbhQAfoX10mwvgEKwrrTpCqA0ADOAHXY/QQnUkdHvl5CIHDkKEChIa3IxC4FBLFq0\nfLxjJk6cTFkuSuAz2zCMlg0btnhlfTxJfv/9KGo0+Qjkt31omOni0oYGQzbq9VlYvnxNXrp0Kc3P\ny2QysV69IBoMRenqGkij0YeHDh1yejtnzpzhwoULuWvXrgTTLkycOIVarRuNxmJ0dfXlnj17nN6H\nl126dImentkpy60oyy3o5ZWT165d4549e2x7184gsICynJ0rV66MV/bXX2dSqzVSr89FL6+c8T7E\nX8disXDKlF9YtWpDNmnyMU+fPu3sU3vnQKRAEJzh9OnTtklMLQF3urpmf+VGJYvFwpkzZ7Ny5Q8Y\nGNg8yYD5fIckH5/8VCjUVKlklitXI01uJknK/PnzqddXJRBju5pezAIFSju1Det8hA8NhtbU6wuw\na9feCR537949njhxItk5/i0WS6p2c2revD2Vyu/jvqW5uAxm27Zd2KJFBwJT7L69rWCFCnXjyr34\nBnLR9vof9PXNK3aWykAQgV5wlvDwcO7YsYN79uxx6kRhdHQ0nz596rT6UmLEiBFUKr+1C2p3Kcue\nTqs/KirKlmb5nK3+UOr1eXjw4EGH64yJiWHHjl9QrZao1Rr5zTdDHAqylSoF0pqJ9Pm5/8mAgMZs\n3vxTAlPtnl/J8uXrxJX7/fffaTS2snudVKv1Dm0GIzgHRK4bwVlkWUZAQACqVKkCFxcXp9Wr0Wjg\n6urqtPpSoly5cpCkFQDuACBcXKaiZMnSaNOmM7JmLYhSparh8OHDDtf/6NEjKBQ6WCelAcAIF5cS\nuHXrVlLFkjRkyA9YujQEsbE3EB19EGPHzkb27MURFNQW//333+srsGnYsCZkeQyAhwDuQ5bHoWHD\nmujduzMkaTisU25LIMu9MWDAZ3Hl8uXLB7P5IIAntmf2QqeTMuz/ofD2yegPRuEdNHToD1SrZep0\nWVioUGlWrVqPWu0ntqvwBTQYvFOUT8ee2Wxm1qz5CMyNW0svy16pWjdesmRVAtttcyX1CHxMYBdd\nXAYzW7YCyR76MZlM7NatN9VqiWq1xB49+sZNAG/fvp116zZjQEDjV8bnSfLLLwdSlrPRza029Xov\nbty40eHzEVIPIteNILxeeHg4wsLC4OHhAVk2wmIJA6AFAOj1H2PatPr49NNPHar7zJkzqF+/Ke7e\nvQmNRovFi39D48aNcejQIRw+fBg5c+ZEo0aNoFQm78t0vXrNsG1bAMiWAIoDuIvnq6JdXathxYoh\nqFu3brL79/xvznYbfbKdOnUKt27dQsmSJZEtW7YUlRWcy5EUCI6uoxeEt5Zer4der4fFYoGLiwoW\ny30AOQAQCsVtW3pkxxQvXhw3bpxHWFgYDAYDlEolfv11Jvr3Hw6L5UOoVHNQp85SrFz5e7KC7aRJ\nP6BSpVqIidmDqKgoADGw/tkSZARUqoT/hA8cOIDFi1dAr5fw+eddkStXLgApD/DPlSxZEiVLlnSo\nrPBuy+AvQMKbKDIyktOmTeO33w7ihg0bnFLngwcP2LFjd1aqFMg+fb5hRERE3Gvffz/Kllp2DLXa\nlixSpCwjIyP5119/MTCwBT/4oGWKd7WyFx0dTY1GT+CCbTgninp9Ee7YsSPZdfz333+cPn06y5at\nRkmqSWA+tdpPWLx4hQTXtW/YsIGS5EPgB7q49KGbW1ZevXrV4XMQ3iwQq26Et1lMTAzLlatBSWpA\nYBhluQB//HFMquqMjIxkgQLvU63uQWAddbqPWLNmw3grV1auXMkePfpw1KjRDAsL45YtWyhJvrax\n9jmUZV+H89E8ePCAGo0r4+fbacqlS5emuK7Y2FiOGTOejRt/zK+/HsTQ0NAEjytRojKBNXHtKZX9\n2a/fQIf6L7x5IAK98DZbt24dDQZ/vtir9gZVKl2C2/+ZTCY+ffr0tUsNg4ODaTSWtQu0MdTpvHjj\nxo1Ey9Su3ZTAfLslhXNYv35zh87JYrEwX76SVCpH29bvb0/1BO3r5M1bisABu/6PY7duvdKsPSF9\nQSyvFN5WJHH69GmYzTKAKNuz1rw3Dx48wA8/jESXLj3xxx9/YNGixTAYPJElix/y538PFy9eTLRe\n65i0fRom699JUmPV1r8l+z8NF4fzsygUCmzbtgbFiq2CUikhS5ZPsWrVohTl6k+p9u2bQ5a/BHAM\nwFbI8ni0bt00zdoThKRk8Ofi2+3Jkyd8+PBhRnfDKWJjY9mwYQtKUk4CBWjNmrmbanUPVqgQwKJF\ny1GrbUNgEnW64lSp3Gz5dSxUKH5mvnwlE607KiqKhQqVpkbzOYFVlKSmrFOncZLfBDZs2EBZ9iPw\nO63PhUMAABDlSURBVIEFlKSs3Lx5c6rPMz3y2TxvZ9Cg4cyRoxgLFizL5ctXpEu7QvqAGLrJ/GJj\nY9miRXuq1XpqNEbWrdsk3uRiRrh79y4DAhpRlj2ZN2/JFCcfmz59OmW5BoEo21DDSKpUnqxWrT7b\ntm1Lrbai3dDLXQJquzQGFqpUOj579izR+h89esRu3XqzatWG/Prr75K1B+z69etZs2YT1qoVxE2b\nNqXofAQhLUGso8/8Ro4cgx9/3IqIiDUAVNDp2qJr17yYPHlMhvWpTJnqOHWqAkymrwHsh17fBWfP\nHo1b0veye/fuYenSpYiNjUVQUBAmTZqGKVN8AXxtOyIE7u61ERMThZiYUjCZVAA2wXp35ipYd648\nDaAogGOQ/9/evYc3VacJHP+2adIkvVmBUihQgUpF8FLxAiiIVCsqXUFQRwadWVFnxDousCoqVqA7\nFJBZLiqjqzgPjuK1lZnR4aJiFwQK6gC6KLQiLBUUGFtKS1pp03f+OKGEEpKSpmkT38/z5HnS/s75\nnd/bpG9Ofpdz7MOpri73e+pgaykvL2fmzNns3r2PwYMzKCv7nq+//pYrrriI3NypWK1W35Uo1UQw\nbyUYCG39wRiSsrLGCrzuNtD2gWRkDGuz9lRVVUlUlNVtAFUkNnasLFu2zOP23333nXTs2F2s1jvF\nYrlfYmI6Sm5urtjtgwSqBRrEZHpKrNZkgffEuEFHssBsgR6u1aHZAjESEzNC7PZOrdI1UVdXJzk5\nUyQuLkkSE1Nk3rzm39lKRKS6ulp69uwnFstvBJZIRERHiYy8U+BtsdlGSWZmtl4YTPkFHYwNf336\npGKxfMzx1zoqqoi0tNQ2a4/VanWdYZS5fuMEdnPWWWd53H7GjHx+/DGN2toEjh0bwNGjeRQVfc7N\nN/fFau1JbGw63bu/Q3Q0GDc1SQY+AOZj3Hp4FfBXIiJm0Lt3JVu2fMLYsae9t43fZszI5+WXP6Wq\najMVFSvJzV3Mm2++1ez9V69ezaFDSRw79jxwHiLJNDQsBcZSU/Mm69cXU1ZW5qsapQJCE32ImTnz\nCXr02ERc3JXExQ0jKelNFiyY1arH3Lx5M0OH3sQFF1zFzJn5OJ3OxrKoqCjy8/Ox24cRGfkYdnsW\nF13UweOyfKfTSWHh+4iYgd7AEmAlP/5YwbJlS9ix41OKi9+lpGQLI0bcQHT0ExjdNU4iIxuAyxvr\nEhmAiIk+ffq0SswFBX/H4cgDUoH+OBwPU1Cwotn719fXExFxfIWtcPIi9EgiIvyfyaPUmdJLIISY\nxMREvvhiI2vXrsXpdDJkyBDi4uJ87+inHTt2MHz4TRw9OhfoybffPk5VVTVPP/37xm2mTHmIiy/u\nz8aNG0lJGc/48eM9Ls3/7LPPqK6OwrghWRRGX3sy1133AACpqSe+mbz00iJqau5j5cpu2GxxjBx5\nA+++uwiH43rAhs02mxEjhrVa3B06JAKlwFAATKZSOnXy/C3Fk8zMTCyWKURGzqahYQAREfuBiYjc\njNX6ZzIyLjjtGIZS4aStu7pUM+Tl/ZeYTJPdxgRKJDGxm191rVmzRuLjB7vV1SBRUZ2ktLS0cZvy\n8nIZPXq8dO6cJpdeeo1s27ZNRIyFR5MnTxWz2SZRUdEybtyEM76tXVMlJSVy/vmXi8lklm7d0mXD\nhg2NZZs3b5aYmI5iNudIdPSvpEOHbl4XWXmya9cuufHG26R//ytlwoQHZPz4e+Syy66VnJwpXmcJ\nKeUNOr3y5+X119+QtLRLpFu38yU3N69V5mnPmTNHzOb73JLzP6RTp55+1XXkyBFJSjpHIiPnCnwp\nZvMk6dfv8pPaPWjQtWKx3Oe6bPCLkpCQfNK9XZ1OZ0BuhFJXVycpKedKRMRCAYfAuxIXlyQHDx5s\n3KakpETmzJkj8+fPb5P7yyrlCZrofz5Wr14tNltXgQ8FPhe7/TLJy5sd8OPs27dPEhO7uu7M9KLY\n7WmycOGzfte3a9cuufrqmyQl5TzJzv7FSYn18OHDYjbHCNQ3frDExWXL22+/HYhQTmlHTEx3tw8w\nkYSE4bJq1aqAH0upQMKPRK999CHqtdcKqal5BMgEwOFYwNKlv2PatEcDepyuXbuyZcsGZs36A//8\n5zpuv30Wt912q9/19erVi6Ki9zyWRUdHI+IEyoFOQAMiP7TossGnk5iYSF1dBfA9xqUWHNTV7aJj\nx44BP5ZSbU1n3bSC+vp68vJmc801N3P33RP54YcfAn6MuDg7kZHu9X5PbGzgEyIYg6QvvLCIgoKl\njUm+traWCRNySErqRZ8+A/jggw887ltTU8OECTl069aXjIyhFBcXn/Y4VquVyZP/k5iY4cDTWK1j\nSEuzkJmZGfCYEhMTeeKJJ7DbBxMd/SAxMYMYNepaMjIyAn4spdqaroxtBePGTeAvf/l/HI6JREVt\noHPnv/L1158HdHbM7t27ycgYTHX1HTidZ2OzLWL58lfJysoK2DG8GT/+XgoKDlJbOxcoxW6/m40b\nP+TCCy88abuxY+/i/ferqa2dAXxJTMxDbNtWTO/evT3WKyK89dZbrFu3iZ49uzFx4v3YbLZWi6Oo\nqIitW7fSq1cvsrOz293qWqWa0pWx7YDD4RCTKVqgyq2f+VopKCiQ+vp6mT7999K370AZODBL1q9f\n36Jj7dmzRx57bJpMmvSwFBcXe922trZWtm7dKqWlpY0rMvfu3Sv33JMj2dnj5E9/WnpGKzVjYzsK\n7GuMMSpqksyeffIYQUNDg5jNNoHDjdvZbPfIc889d+bBKqVERPvo25lTL3M7dWouixcX4XDMBfZw\n3XU3s3lzEf369QNg586dlJSUkJ6e3qyFQKmpqcyaledzu7179zJkyPVUVERQX1/ByJFZLFw4m4sv\nHkRl5V04nZfy0UdzKSvbz5NPTvVYx8GDB5kz57/Zt+8Q2dmZWK0xVFfvA4z7h5rN3xEbe/JZekRE\nBBaLjbq6A0CC8VeJPIDdbvfZZqVUeGjrD8ZWM3bsnWKzjRB4T6KiHpfk5F5y+PBh6dChh8BOtzv/\nPCzTp88QEZEFC54Vmy1J4uNvEJstSRYtWhyw9lx99U1iMuW5jntU7PYrZdy4cWK13uU262SnxMd3\n9rh/RUWFdOnSW8zmHIHnxW4/X0aNuk3s9hSBPLFY7pLu3dOloqLilH3nzVsgdntvgT+IxfJrSU3t\ne9o7IymlfEPP6NuH1157iRkz8lmz5hnOOSeFefPWkpCQgNlsAY40bmcyHSE6+mz279/P1KlPUlu7\nhZqaVGAPjzxyCbfeOprk5OQWt+err77C6Zzv+smOw/FvlJWtoKEh3W0rG05nvcf9CwsLqay8kLq6\nZwBwOLJYtSqD994rZMWK1Zx9dl9++9sFHq9vM2XKQ6SlncOKFWvo0qUXDz00v1VX8iql2pe2/mAM\nusWLXxC7vafA82IyTZXExK6yb98+2bRpk8THX3LSnO74+Ivl008/Dchxhwy5QUymfFfdNWK3D5H8\n/HyJje0k8KzAh2K3XyU5OVNO0+7FYrP92q195WI22/Tqi0q1AfR69O3f8uXLeeONv9GhQzyPPjqJ\nHj16UF5eTo8e6Rw9WggMAdYSGzuGsrLS014F8kzs2bOHq67K4sgRK05nOVlZV/POO6+wfft2pkx5\nikOHyhk9+nqmTXsUk8l0yv579+6lf//LqKrKAy7AZsvjlltSePXVF1vcNqXUmfFn1o0m+nZi9erV\njBnzSxoaLERGHqOwcJnHK0D6q6amhu3btxMbG0t6evoZTyPctm0bDz74OAcOHOLGG4czZ85MLBZL\nwNqnlGqeYCb6SGAxcCHwE8ZlCHe5lWcDTwL1wMvASx7q0ETfxE8//cSBAwfo3Lkz0cYF2ZVS6iTB\nTPS3ACOBu4ErgMeAUa4yM/AVcCngANa7tj3YpA5N9EopdYb8SfT+XgLhSmCl6/kmjKR+XF/gG6AS\nqAM+4fhFvZVSSgWdv4k+Hvd5gsb94yLdyirdyqo4vlpGKaVU0Pmb6I8A7pOhI4EG1/PKJmVxQIWf\nx1FKKdVC/i6YWo8x4Po2MBD4wq1sB3AukAgcxei2edpTJdOnT298PmzYMIYNG+Znc5RSKjwVFRVR\nVFTUojr8HYyN4MSsG4B/BwYAscCLGIOvuRhn+kuAP3qoQwdjlVLqDOk8eqWUCnPBnHWjlFIqRGii\nV0qpMKeJXimlwpwmeqWUCnOa6JVSKsxpoldKqTCniV4ppcKcJnqllApzmuiVUirMaaJXSqkwp4le\nKaXCnCZ6pZQKc5rolVIqzGmiV0qpMKeJXimlwpwmeqWUCnOa6JVSKsxpoldKqTCniV4ppcKcJnql\nlApzmuiVUirMaaJXSqkwp4leKaXCnCZ6pZQKc5rolVIqzGmiV0qpMKeJXimlwpwmeqWUCnOa6JVS\nKsz5k+htQAGwFngf6Ohhm0lAseuR63frlFJKtZg/if5+YBswFHgFmNakvBcwDhgEDASygAta0MaQ\nVFRU1NZNaFXhHF84xwYa38+RP4n+SmCl6/lK4Nom5XuB6wFx/WwGavxqXQgL9zdbOMcXzrGBxvdz\nFOWjfALwH01+dwA44npeBSQ0Ka8HyoEI4GngH8A3LWumUkopf/lK9EtcD3cFQJzreRxw2MN+VuBl\noBKY2JIGKqWUapkIP/aZjJHgZwC/AIYADzSpcyXwETDXSz3fAL39OL5SSv2c7QLSzmQHfxK9DVgK\ndAF+whh4PYgx0+YbwAS8Dmx0q/8xjBk4SimllFJKKaWUUqpdas5CK4BOQAlgCVK7WiISeB7YAHzM\nqWMO2cBmV/k9wW1awPiKEcAOrAfSg9iuQPAV2x0YXY6fAH/Ev67OtuQrvjEY789NwO+C27QWa877\nEuB/gPxgNSqAfMU3Cfg/V9nHQJ+gts6LyZxYJXs7sMDDNtcDWzBm8oRCor8FY3YRwBXAcrcyM1CK\nMf3UjPEPlRTU1gWGtxgBLgU+A/bTjt5szeQtNhvGmJPV9fMyjA/uUOItPhPGCVUcRlLZAZwd1Na1\njK/3JcBvMBLlrGA1KoB8xfdnIKO5lQXzWje+FloBOIFMoCJYjWoh95g2YSS94/piJIpKoA7jrHBo\nUFsXGN5iBOMDeRSwM5iNChBvsdVirO6udf0cRegt/PMWnxM4D2MtTCeMxH8sqK1rGV/vy8HA5cAL\nhN43MfAd3wDgcWAdMNVXZa2V6CcAXzZ5JOB9oRXAhxiLrUJFPCdiAuOfJ9KtrNKt7HQxt3feYgTj\njOm7oLYocLzFJsAh1/MHgRiM92co8fXaNWCcOW7B+PrvCF7TWsxbbF0weg9yCM0kD75fu9cxvrEM\nB64CbvJWWWsl+iUY17dxf1Tie6FVqDnCiZjA+Hs2uJ5XNimLI3S+qbjzFmOo8xVbJDAP41vmmCC2\nK1Ca89oVAilANHBXkNoVCN5iG4sxBvh34FGMKeChFBv4fu0WYpwU12GMeXrtxglm18164EbX8xsw\nBmVDnXtMA4Ev3Mp2AOcCiRjdG0Mx1haEGm8xhjpfsb2AkQBHc6ILJ5R4iy8e+F+M96YARzHOGkOF\nt9iewejquAaYjTG+8kpQW9dy3uJLwOglicH4xjIcY5ysXbABb2H0KX3IiYHJSZw6yPUtoTEYG4Ex\nG2O969EHY6bGva7ykRiDsJ9hXPUzFPmK8bh2NfLfTN5iy8BIfB+7PUa1TTP95uu1uxdjVtE64FlC\nq5ujue/LXxGag7G+4rsDI7esA55qiwYqpZRSSimllFJKKaWUUkoppZRSSimllFJKKaWUUkoppVSb\n+BcY2oKX/xwnNAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f414c8a7710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter(np.dot(x,best_theta[0,:].T),\n",
    "        np.dot(x,best_theta[1,:].T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "scatter(best_theta[:,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesk Algorithm #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from nltk.corpus import wordnet as wn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stoplist = stopwords.words('english')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def0 = set([word for word in wn.synsets('plant')[0].definition().split() if word not in stoplist])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def1 = set([word for word in wn.synsets('plant')[1].definition().split() if word not in stoplist])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lesk_counts = [[len(set(instance).intersection(defset)) for instance in instances] for defset in [def0,def1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f876586cf10>"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEWtJREFUeJzt3d9rXGd+x/HPN2hFHUo6lhxc2oLjSXu7rIQKvljIUMkx\nvcjFRl7vZViI/A+43sYXJYI6rIzrvQtsrIWmd0lZd2+ttbRMYO+qrM1CIbDeSQNlKSaSp2nBG0r5\n9kKPpPFoRvPjOTPzHc37BYPOmec853zn6Oh8fM48MzZ3FwAA0bww6gIAAGiFgAIAhERAAQBCIqAA\nACEVElBmduuYtpX0WCtiWwCAyZAdUGZ2VdJym7ZFSZvuvi6pnOYBAOgoO6Dc/a6kWpvmsqSlNF1L\n8wAAdDQ1yJWnK6d985I+GuT2AAAnx1AGSZjZvKRP3f3RMLYHABh/wxrFt+juN4a0LQDACTCQW3xm\nVnL3epq+6u630/Siu281Lct3LQHACebu1k+/IkbxXZa0YGZvNzy9mdqWJK2Z2WMz25XUMozcPfzj\n3XffHXkNJ6FG6qTO6A/qLPaRI/sKyt1/KumnTc8tpJ+bkmZytwEAmDx8kwQAICQCqkuVSmXUJXQ0\nDjVK1Fk06iwWdcZhufcIswsw81HXAAAYDDOTj2qQBAAAg0BAAQBCIqAAACERUACAkAgoAEBIBBQA\nICQCCgAQEgEFAAiJgAIAhERAAQBCIqAAACERUACAkAgoAEBIBBQAICQCCgAQEgEFAAiJgAIAhERA\nAQBCIqAAACERUACAkAgoAEBIBBQAICQCCgAQEgEFAAipkIAys1vHtC2b2aKZrRSxLQDAZMgOKDO7\nKmm5Tdu8JLn7Vpqfy90eAGAyZAeUu9+VVGvTfEXS0zRdk7SUuz0AwGQY9HtQJUm7DfOzA94eMBE2\nNjb0+uvLev31ZW1sbJzYbWKyTQ1hGzaEbQATY2NjQ9/5zlt69mzvrd9f/vIt/exn/6RLly6dqG0C\ngw6ouqSZNH1a0s6AtweceHfu3E1B8ZYk6dmzvecGGRaj2CYwkIAys5K71yV9LGlB0pak85IetFp+\ndXX1YLpSqahSqQyiLADAgFWrVVWr1ULWZe6etwKzy5LuSvqBu/8kPbft7gtpekV7AyTK7r7eor/n\n1gBMkubbbadO/e3Qb/ENY5s4GcxM7t7XWz3ZAZWLgAJ6t7GxoTt37kqSrl27OpSgGMU2Mf4IKABA\nSDkBxVcdAQBCIqAAACERUACAkAgoAEBIBBQAICQCCgAQEgEFAAiJgAIAhERAAQBCIqAAACERUACA\nkAgoAEBIBBQAICQCCgAQEgEFAAiJgAIAhERAAQBCIqAAACERUACAkAgoAEBIBBQAICQCCgAQEgEF\nAAiJgAIAhERAAQBCIqAAACFN5a7AzJYl1SWV3X2913YAAFrJuoIys3lJcvetND/X1D4nqZbaa83t\nAAC0k3uL74qkp2m6JmmpxTK30s+yuz/M3B4AYELkBlRJ0m7D/GxjYwqkz81st2k5AACOVcQgCWvb\nYFaS9FjSiqR1MztfwPYAABMgd5BEXdJMmj4taaepfUXSB+7+lZnVJV2WdLt5JaurqwfTlUpFlUol\nsywAwChUq1VVq9VC1mXu3n/nvUEPC+6+bmbXJT1w90dmVnL3upldd/fbDcuvNI/kMzPPqQEAEJeZ\nyd3b3mk7tm9uOJjZivYGSBwMIzezbXdfSNPXU/tMm2HoBBQAnFAjDahcBBQAnFw5AcU3SQAAQiKg\nAAAhEVAAgJAIKABASAQUACAkAgoAEBIBBQAIiYACAIREQAEAQiKgAAAhEVAAgJAIKABASAQUACAk\nAgoAEBIBBQAIiYACAIREQAEAQiKgAAAhEVAAgJAIKABASAQUACAkAgoAEBIBBQAIiYACAIREQAEA\nQiKgAAAhEVAAgJCmcldgZsuS6pLK7r7eon1e0nlJM63aAQBoJesKKoWP3H0rzc+1WOwdd78nqdSm\nHQCAI3Jv8V2R9DRN1yQtNTaa2WVJ/ypJ7n7b3R9mbg8AMCFyA6okabdhfrapfUHSrJnNmdn1zG0B\nACZIEYMkrEP7l/tXTun9KgAAOsodJFGXNJOmT0vaaWrfkfR5w7J/Kele80pWV1cPpiuViiqVSmZZ\nAIBRqFarqlarhazL3L3/znuDHhbcfT3dwnvg7o/MrOTudTM7L+myu99O7b91939pWofn1AAAiMvM\n5O6d7rS1lHWLr+HW3aKkurs/Sk2bqf1zSfV0a2+mOZwAAGgn6wqqkAK4ggKAE2tkV1AAAAwKAQUA\nCImAAgCEREABAEIioAAAIRFQAICQCCgAQEgEFAAgJAIKABASAQUACImAAgCEREABAEIioAAAIRFQ\nAICQCCgAQEgEFAAgJAIKABASAQUACImAAgCEREABAEIioAAAIRFQAICQCCgAQEgEFAAgJAIKABBS\ndkCZ2bKZLZrZSoflruduCwAwObICyszmJcndt9L8XJvlliRdzNkWAGCy5F5BXZH0NE3XJC21Wc4z\ntwMAmDC5AVWStNswP9u8gJnN7V9hAQDQrSIGSViH9pkCtgEAmDC5AVXXYQCdlrTT2MjVEwCgX1OZ\n/T+WtCBpS9J5SQ8kycxK7l6XVDazsvZu/c2kwHrYvJLV1dWD6UqlokqlklkWAGAUqtWqqtVqIesy\n97zxC2l4eU1S2d3X03Pb7r7QtMwPJH3X3R819ffcGgAAMZmZ3L3TW0Gt+446HAgoADi5cgKKb5IA\nAIREQAEAQiKgAAAhEVAAgJAIKABASAQUACAkAgoAEBIBBQAIiYACAIREQAEAQiKgAAAhEVAAgJAI\nKABASAQUACAkAgoAEBIBBQAIiYACAIREQAEAQiKgAAAhEVAAgJAIKABASAQUACAkAgoAEBIBBQAI\niYACAIREQAEAQiKgAAAhEVAAgJCmcldgZsuS6pLK7r7eon0lTb7q7u/kbg8AMBmyrqDMbF6S3H0r\nzc81tS9K2kzBVU7zAAB0lHuL74qkp2m6Jmmpqb3c8FwtzQMA0FHuLb6SpN2G+dnGxqZbfvOSPsrc\nHgBgQhQxSMI6LrB3K/BTd39UwPYAABMg9wqqLmkmTZ+WtNNmuUV3v9FuJaurqwfTlUpFlUolsywA\nwChUq1VVq9VC1mXu3n/nvUERC+6+bmbXJT1w90dmVnL3elrmqrvfTdOL+wMqGtbhOTUAAOIyM7l7\nxzttrWTd4nP3h6mARUn1hlt4m+n5JUlrZvbYzHYlkUQAgK5kXUEVUgBXUABwYo3sCgoAgEEhoAAA\nIRFQAICQCCgAQEgEFAAgJAIKABASAQUACImAAgCEREABAEIioAAAIRFQAICQCCgAQEgEFAAgJAIK\nABASAQUACImAAgCEREABAEIioAAAIRFQAICQCCgAQEgEFAAgJAIKABASAQUACImAAgCEREABAEIi\noAAAIU3lrsDMliXVJZXdfb3XdgAAWsm6gjKzeUly9600P9dLOwAA7eReQV2R9PM0XZO0JOlhD+1j\nYWNjQ3fu3JUkXbt2VZcuXRpxRXHk7JvmvpLarqvddo5bx2uvzevevQf64ov/0Llzf6zl5b/WJ5/8\nSrXab/TkSV3f+Ma03njj2/r1r/9dn332b/r6669l9gd65ZU/1fvvr7XdfvN6f/jDv2v7ut977z39\n6Ef/KEl6441v63e/+299+eWOvvpqV0+e/Jfc/1dnz57WSy+9rDNnZvXaa/P65JNftX2dtdpnevLk\nf1K/syqXy8/1Oa7/jRt/ry+++E+dO/dn+uY3X9G9ez/X73//fzp37k/0/vtr2t7e1q1bH+jZs691\n7txZff/73zuyrub90NjevO/btXX6XRepn+Mz6t971LoGyt37fkj6saRvpelFSWu9tKfnPbL79+/7\nqVNnXfrQpQ/91Kmzfv/+/VGXFULOvmnuOz1d8unpl1uuq912Oq1Desmla2l6xqUX0/yZpmWW08/D\n56emZtts/1padr//GZ+eLrV83Tdv3mxY9mi/w1pebKjzsOajr3O5of1Mm/W27j89XWrxug/7vfDC\nH7ao78Xn1nXz5s22+2F6+uW0jU5tx/+ui9TP8Rn17z1qXd1I5/j+Mqbfjn4YQHN+fEC1bfcxCKiL\nF99MB4Wnx4d+8eKboy4rhJx9c7TvhbbraredbtYhvdkwfSHNNy/zag/bb9X/QsvXPTPzaod+bzbV\ndbTm51/nq039jltvc/9W++ZCh333fF2dX8+FLtqO/10XqZ/jM+rfe9S6upETULm3+OqSZtL0aUk7\nPbZLklZXVw+mK5WKKpVKZlkAgFGoVquqVqvFrKzfZNsLRs1JWknT13V4O690XHvTOgYZ3tnG+dJ6\n0LjFxy0+bvENR9S6uqGMKyjb698/M1vR3gCIg2HkZrbt7gvt2pv6e24NgzaRb052iUESDJJot+8Z\nJFGsqHV1YmZyd+ur76jDYRwCCgDQn5yA4pskAAAhEVAAgJAIKABASAQUACAkAgoAEBIBBQAIiYAC\nAIREQAEAQiKgAAAhEVAAgJAIKABASAQUACAkAgoAEBIBBQAIiYACAIREQAEAQiKgAAAhEVAAgJAI\nKABASAQUACAkAgoAEBIBBQAIiYACAIREQAEAQiKgAAAhEVAAgJAIKABASAQUACCk7IAys2UzWzSz\nlTbtK+mxlrutUapWq6MuoaNxqFGizqJRZ7GoM46sgDKzeUly9600P9fUvihp093XJZXT/Fgah4Nh\nHGqUqLNo1Fks6owj9wrqiqSnabomaampvdzwXC3NAwDQ0VRm/5Kk3Yb52cbGdOW0b17SR5nbAwBM\nCHP3/jub/VjSB+7+MN2+u+ju77RYbl7Sd939Rou2/gsAAITn7tZPv45XUG0GP+y6+z1JdUkz6bnT\nknbarGaxVThJ/RcOADjZOr4H5e7rLR73UvPHOnxf6bykB5JkZqX9/mZ21d1vp+mxHSQB4GToNPI4\nLXOr1z5F67POW+nn0OocpKxBEu7+UDoInrq7P0pNm+n5JUlrZvbYzHYlcTuvYONyEI/LSSGiLj7K\ncaQ94gm1TZ1DPRY7jTxOz12VtNxLnwh1Jitm9htJvx10jQ119PxRo26Pz+zPQaUrqq3GARHuvpB+\nbrr7jLv/efr5iz5fTLgTas5OL7DGsTiIx+WkkLYT6iTbxUc5jrRHPKEe0z7sE2qnkcdy97upres+\nA9BPnZK04u5/4e6/GHB9kvr7qNH+Mt0cn0P9JomMz02FOqHm7vQCjcVBrDE5KQQ9yXbaD63ar2jv\n/eF2fQahnzql4R+Lx448LrBPrn63OZPOR9cHUFMr/XzU6Hvq8vgc9lcd9fpizqfpaCfUrJ1eoHE5\niMflpBDxJNtpP7Rqj3hCbdc+7GNRkvoZmDWKwVw9b3P/jpakWRvOe/4dP2rUcHdtXtJ26rPTrk+j\nYQdUPy9GCnZCzd3pBRuHg1gaj5NC1JNsp/0QZSRsz3WO4FjsduRxbp9cPW8zveWwf5t8R8P7YoSO\nx1+6+/Dp/riFbvpI+R/U7UcvL+aRdPiBXzO7aGaL+7dYBqynnW5mXfXpuYhihvk3r2+///5BnL0/\ni66zzz5F6OskKw3s+Oy0HxrbG/+RFO2EeqTOQR2LHXwsaSFt57mRx+5e76VPwDprOvxH/ayGU2c/\nHzXq+m+78IAq6ER18GKCn1D72um9aPo2jmZhDuIB1DmQk0Lm730UJ9lO+66xvZzarVWfAeu1zk3t\n7cuhnlDTPyYXrPXI44VU82VJC2b2trv/5Jg+0ercSoN4JOnLYdSpLv627ehHjbr/23b3oT0kzWnv\nfr0kXZf0rTRdaljmasP0Ynr8UZpf2+8TsM6WfYZQ60ra/krDc9sN05e1d1vq7YbnltPjb4b4u++n\nziN9Rvl7b9E+N4zjs4t916p9qPsuo86hH4s8hvd71977tLuSHqeff9WuT6tH1lcd9SP9i7MmqeyH\nt0a23X3B9j439c/phcxIuuxpaHrqft7d/yFwnUf6YLwc93s/pn3oxycwCYYeUAAAdIP/URcAEBIB\nBQAIiYACAIREQAEAQiKgAAAhEVAAgJAIKABASP8PWXNPnFGwlMUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f87658937d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter(np.array(lesk_counts[0])+0.05*np.random.randn(len(instances)),\n",
    "        np.array(lesk_counts[1])+0.0*np.random.randn(len(instances)))"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}