{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "\n",
      "np.random.seed(42)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Problem 1"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Generate points on a circle."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "angles = np.random.random(10000) * 2.0 * np.pi\n",
      "radius = 1.0\n",
      "x = np.cos(angles) * radius\n",
      "y = np.sin(angles) * radius\n",
      "plt.scatter(x, y, s=1)\n",
      "plt.gca().set_aspect('equal')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAH0CAYAAABYRNWcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4lNXB/vHvMzOZmeyEfUeUCuKCouICAhYBtaKtrVVR\nXPurWq0U37dWqq1SFd9qtbVVwfat4G4X3Fe0NQXFBRdcUREMe0CWkHWWzDy/P86EwDsIJJnkzHJ/\nrisXyZCJNyPM3HPOec5xXNd1EREREdmBx3YAERERST8qCCIiIpJEBUFERESSqCCIiIhIEhUEERER\nSaKCICIiIklUEERERCRJSgvCggULOPXUU+nbty8ej4f7779/t99fUVGBx+NJ+pg/f34qY4mIiEgL\n+VL5w+rq6jjkkEM4//zzOe+883AcZ6/u99JLLzFs2LDtX5eVlaUyloiIiLRQSgvCSSedxEknnQTA\nBRdcsNf369y5M927d09lFBEREWmDtFiDcPrpp9OjRw9GjRrFvHnzbMcRERHJeVYLQnFxMbfffjv/\n+Mc/eOGFFxg3bhxnnnkmDz/8sM1YIiIiOS+lUwwt1aVLF6ZNm7b96+HDh7N582ZuvfVWzjnnHIvJ\nREREcpvVgrArRx55JPfdd98uf2/QoEEsX768gxOJiIhknv32248vv/yy1fdPizUIO1qyZAm9e/fe\n5e8tX74c13X10YKP66+/3nqGTPrQ46XHTI9Zen7oMWv5R1vfUKf8Msdly5YBEI/HWblyJUuWLKFL\nly7069eP6dOns3jxYl555RUA7r//fvx+P4ceeigej4dnnnmGe+65h1tvvTWVsURERKSFUloQFi9e\nzLe//W0AHMfh+uuv5/rrr+eCCy7gvvvuo7KykhUrVmz/fsdxuOmmm1i5ciVer5fBgwczZ84cJk+e\nnMpYIiIi0kKO67qu7RB7y3EcMihuWigvL2fs2LG2Y2QMPV4tp8es5fSYtZwes5Zr62umCoKIiEgW\nautrZtotUhQRERH7VBBEREQkiQqCiIiIJFFBEBERkSQqCCIiIpJEBUFERESSqCCIiIhIEhUEERER\nSaKCICIiIklUEERERCSJCoKIiIgkUUEQERGRJCoIIiIikkQFQURERJKoIIiIiEgSFQQRERFJooIg\nIiIiSVQQREREJIkKgoiIiCRRQRAREZEkKggiIiKSRAVBREREkqggiIiISBIVBBEREUmigiAiIiJJ\nVBBEREQkiQqCiIiIJFFBEBERkSQqCCIiIpJEBUFERESSqCCIiIhIEhUEERERSaKCICIiIklUEERE\nRCSJCoKIiIgkUUEQERGRJCoIIiIikkQFQURERJKoIIiIiEgSFQQRERFJooIgIiIiSVQQREREJIkK\ngoiIiCRRQRAREZEkKggiIiKSRAVBREREkqggiIiISBIVBBEREUmigiAiIiJJVBBEREQkiQqCiIiI\nJFFBEBERkSQqCCIiIpJEBUFERESSqCCIiIhIEhUEERERSaKCICIiIklUEERERCSJCoKIiIgkUUEQ\nERGRJCoIIiIikkQFQURERJKoIIiIiEgSFQQRERFJooIgIiIiSVJaEBYsWMCpp55K37598Xg83H//\n/Xu8z0cffcSYMWMoKCigb9++3HjjjamMJCIiIq2Q0oJQV1fHIYccwp133kl+fj6O4+z2+6urqxk/\nfjy9evXinXfe4c477+S2227jjjvuSGUsERERaSHHdV23PX5wcXExd999N+edd943fs+sWbOYPn06\nGzZsIBAIAHDzzTcza9Ys1qxZkxzWcWinuCIiIlmlra+ZVtcgvPHGGxx33HHbywHAhAkTWLduHStX\nrrSYTEREJLf5bP7HKysr6d+//0639ejRY/vvDRgwwEYsEdnBpk2b+Pe//02nTp3417/+xapVq3j+\n+ecpKiqiW7duVFVV4TgOXq+XWCzGyJEjmTBhAj169GDEiBGUlZXZ/iOISCtYLQh7WqMgIqkXjUaZ\nOXMmf/zjHykpKaGiomKH3y0CXCAOOInPQ0Aw8XV94vsKqa5ex7p1NUANUAhEgQgVFet5+OF5gBeI\nJO4bBfISX/uBusR/wygrK+OKK67g2muv3WlEUUTssVoQevbsSWVl5U63bdiwYfvv7coNN9yw/fOx\nY8cyduzY9oonktE+/PBDxo8fz8aNGxO3FGL+yccxL9RhtmypBfIxL/4OUIuZeYwDBZiC4KN5NrIw\n8XkjUIp54S9JfJ8DFCe+r+l78nZIVIMpIDWJn91UNvLYunUrN954IzfeeEvivk0lwqWkpITbb7+d\niy66CI9HV2aLfJPy8nLKy8tT9vOsLlKcPXs2v/jFL9i4ceP2dw0zZ85k1qxZrF69OjmsFimKJIlG\no0yaNImXXnopcYsfCAANmBf3EOaFOY55US4GYk33Tnxu3s37fD5KSkq46KKLcByH733vexx55JH4\nfK1/L7F8+XKeeOIJXnnlFT7++GPWrl1Lfn4+DQ0NmFGGPJoLSv0O92wqIzWYUYhGTLGoBmD06NH8\n/e9/3z4tKSI7a/NrpptCtbW17vvvv+++//77bkFBgfub3/zGff/9991Vq1a5ruu611xzjTtu3Ljt\n379t2za3Z8+e7llnneV+/PHH7rx589ySkhL3jjvu2OXPT3FckYz0+OOPu5i37C7ku1Cww+e44Ev8\n6k3c5ncB1+fzuT169HBnzpxp+4/wjT7//HP30EMP3eHP53WhOPFnDOz2zwe4s2fPtv1HEEkbbX3N\nTOkr7quvvuo6juM6juN6PJ7tn1944YWu67ruBRdc4A4cOHCn+3z00Ufu6NGj3WAw6Pbu3dv9zW9+\n881hVRAkB/3nP/9xy8rKXHBcCLrgSbwg5iVeNJ3EC2ieC7her9c999xz3draWtvRU+ovf/mL6/M1\nlYNAoiA0FYWmkuAkCoUpE2eddZbt2CLWtPU1s92mGNqDphgkF9TX13PsscfywQcfYKYKIpipgugO\n31UENODzOfzlL3/hggsusJDUvk8++YSjjz6a2to6zGPkYqYimgQw0xQxfL44ixcv5tBDD7URVaTD\ntfU1UwVBJA0sXLiQ0aNHY9YPBDGloGntgEvTXHx+fj6vvvoqRx11lL2waaympobevXtTW1uLWcPQ\npC7xawBwyctzmT59OjNmzOj4kCIdRAVBJEPNnz+fiRMnYl78m6bR8zAjBQ7go6ysiCeffDJRHqSl\nvv76a4466ii++uorTPEKJX6nIPGrB6jlBz/4Af/4xz+sZBRpLxm9k6JIrvnggw9wHAfHCTBx4omJ\nWwswxaAYiDJ06FDi8RiuG2HLli0qB23QrVs3VqxYgeu6uG4DkydPTlwq6WKKWC0Q4J//fAHH6YTj\nODowTiRBIwgi7ayxsZGDDjqIzz//AlMEYokPL2bePMw111zDLbfcYjNmznn44Yc599xzMVMRYcx+\nEDWJX70UFcFbb73F0KFDbcYUaTWNIIikqVmzZuE4Dnl5xXz++eeYd61NOwnCmDGjcN0QruuqHFhw\nzjnnJEYWannwwTk0b+Bk9ouorY1z4IGH4jg+zjnnHLthRSzQCIJICkWjUXr06MHWrY2Y1fQRzMI4\nB6jjtNNO48knn7SaUXZv1qxZ/OQnP6F5h8iaxO+Y3R/ffvttjjzySGv5RPaWFimKpIGXXnqJE088\nkeZtiqF5x8J6amtrKSws/Mb7S/qJxWIMGDCAtWvXYkpeGDMdEQca+K//+i9+97vfWc0osjuaYhCx\naOHChTiOlxNPPI3mF49iwOGUU8bgunW4rqtykIG8Xi9r1qzBdV0eeuivTbcmPjzcfvs9OI6XgQMH\nWkwp0n40giDSCocccggfffQlze8qmy5TrGP16tX07dvXaj5pH42NjRQXFxMK+TCjQw00nRlRVuZj\ny5YtdgOK7EAjCCIdaPjw4ThOPh99tAKzxiAfiLHvvt1pbNyG67oqB1nM5/PR0NCA69bwi19cSfMh\nUzVs3boVx/HTvXt3vZGRrKARBJG98K1vfYsvv6zAlIKmDXcCjBgxjLfeestqNrHr7rvv5oorfoop\ni/WYvx95BAIRQqHQ7u8s0o40giDSji699FIcx8OXX67DlAMztPyjH/0I1w2pHAiXX345rhvnww/f\nTNxiFqqGw2EcJ8Dhhx9uM55Iq2kEQWQXHnnkEc4553xMKfDTtLHRpEkTefrpp+2Gk7Q2b948fvCD\nH2Iui2w+E2LIkH4sXbrUZjTJMRpBEEmhZcuW4ThezjnnUszcchCAgw8eiuuGVQ5kj77//e/jujHm\nzPkTZufMOBDms88+w3HyuPbaay0nFNk7GkEQwWxw5Pf7Me/2Iph9+vMYMKA3FRUVVrNJZjv55JN5\n4YVXMNNTDTTtqbBq1Sr69etnN5xkNW2UJNJGY8aMYcGCtzCXLPqAfLp2DfD1119bTibZpE+fPqxb\ntw6zmLEBCOLxRIjFYpaTSbbSFINIK7366qs4jj9RDgCKgEZWr/5U5UBSbu3atcTjcRwnhFnXEkp8\nHeD888+3HU8kiUYQJCeVlJRQUxPHrDOoBYJccsn5zJ4923IyyQXLli1j//0HY6a0zN8/CFFVVUVp\naandcJI1NIIg0gKTJ0/GcbzU1NRgrlk3m9+4boPKgXSYb33rW7hunGOPPQSzJsFMM3Tq1ImJEyda\nzSbSRCMIkjP8fj/RaB6mGBQBtVRUVDBgwADLySTXOY6DGU2ow5zlUcPHH3/MgQceaDeYZDSNIIjs\nwVlnnYXj7FgO/Awbth+u66ocSFpwXZfy8ucwaxPM8dIHHXQkgwcPtppLcptGECSrBYNBwuEYTRse\neb0xampqyM/Ptx1NZJfMaEIe5lJbs0lXKFRHIBCwG0wyjkYQRHbhiy++wHH8hMNhzJNskBNOGE1j\nY6PKgaQ113V5+ul5ia8iQIBgsJhf//rXNmNJDtIIgmSd/v37s3r1NqCapn3xa2trKSwstJxMpGVK\nS0upro4kvnLIz4f6+nqrmSRzaARBJCEej+P1+li9ejXmnVcx48ePw3VdlQPJSNu2beOaa36G2a45\nRkNDA47j48MPP7QdTXKACoJkhXnz5uH15hOPxzC7IcIDD9zN/Pnz7QYTaaNbbrmFhoZtmHU0hUCM\nYcNGcemll1pOJtlOUwyS8fbff3+WLVuF2Sq5BMepTRQFkexywAEH8Nlna2g6erxr16B2/ZRvpCkG\nyWl+v59lyypoOkdhxIghKgeStZYuXcq9996O2Viplk2banEcH5FIZE93FWkxFQTJSOFwGMdxEnsb\nRAEvv//9bbz11lt7uqtIRvvxj39MQ0PTAtwQECQQKOXxxx+3nEyyjaYYJONUVFQwcOCBmHMUIpjF\nW3UEg0HLyUQ6lsfjSTwnFgBxJk0az9NPP207lqQJTTFITrn66qsZOHAIphw0EAh4cN2YyoHkpHg8\nntiO2Rw89swz/2LMmDG2Y0mWUEGQjHHEEUdw221/xBxuE+EHPzidUChkO5aIVR9//DE33XQdZvfF\nEAsWLMDj0VO7tJ2mGCQj9OrVi8rKjZhdEQP85jf/xa9+9SvbsUTSxrZt2+jUqSvmMl9TnPV8mdva\n+pqpgiBpz+xNb05fhADvvruI4cOHW04lkn4ikQiBQEHiqxjgob6+VtuL5ygVBMlqphyAeVfUSHV1\nNcXFxTYjiaQ98++mgKb9EpYufZchQ4ZYTiUdTYsUJSs1NjbiOH7Ai3mii+G6rsqByF5wXZeiIg+m\nINRzwAFH8Pbbb9uOJRlGBUHSTigUIi8vH/PXMwbU47pxy6lEMktNTQ1+vw+zPXOIo44ayXPPPWc7\nlmQQFQRJK1u3biU/vwjzzidKIBDQtJJIK4XDYQYN6oUZiWvklFPO4corr7QdSzKECoKkjbq6Ojp3\n7ooZNSjC63V0GaNIGy1btoyxY48FioEQf/rTHGbNmmU7lmQALVKUtGAu0eqBeacTJxiEhoYG27FE\nssY+++zDypWbABfw8MorTzJu3DjbsaQdaZGiZLxIJEKnTp0wBy6F6NatWOVAJMUqKir4/vdPxBSE\nGCeccCLPPvus7ViSxlQQxKpoNEog0LRNcgF5eV42btxoNZNItvrnP//JNddMxWzN3MikSd/jhRde\nsB1L0pSmGMQqxwlgnqzyCQQiWnMg0gHOPfdcHn74USAIePngg9c45JBDbMeSFNNGSZKxHMdH0/7x\nPp+PaDRqO5JIzhg8eDBffLECc7YJbNmymrKyMruhJKVUECQjmZ3eHMx8qPaMF7Fh4MCBVFSsxmxG\nFiISqSMvL892LEkRFQTJOIFAgEjEjzlbwYvrNtqOJJKzSkpKqKmpwyxJi+O6MduRJEV0FYNklKFD\nhxKJ5GHKgUflQMSy6upqzDqgPCCO43gtJ5J0oYIgHWbKlCksXbocqAP8hMO6lFEkHZh3meHEVwV4\nvSoJooIgHeS5557joYf+DkSAACtXLsPv99uOJSIJzVMLtcTjcPnll1vNI/apIEi7i0QinHLKmYAf\nCPDgg/9L//79bccSkf+jtrYWc7hTIffcM4c77rjDdiSxSIsUpd05jgdzxYKfUaOOYOHChbYjicg3\nWLFiBfvtdxDQABRSVbWW0tJS27GkFbRIUdKauZzRBQo47rgjVQ5E0ty+++7LtGmXYg53qqdTp862\nI4klGkGQdtOjRw82bqwGokBM/+9EMsixxx7LG28sxkwN1uvfbwbSCIKkpbvvvpuNGzdjyoGPxkZd\nziiSSRYtWgQ0AvWAlx49elhOJB1NBUFSLhQKccUVUwGzKvq9997QZVMiGci8+8wDPGzcWMu9995r\nO5J0IE0xSMqZMxaCQIwzzpjE3//+d9uRRKSVFi9ezIgRYzCLFvMJh6t0iXKG0FbLklbKysqoqgpj\nnkx0xoJINpg0aRLPPvsfzCZncf27zhBagyBp46GHHqKqKoQpBz49iYhkiWeeeQbHqQXyAdhvv/3s\nBpIOoREESYlIJEIgEMRc0ujnq68+Z5999rGcSkRSyXHyMXuaRFi4sJxRo0bZjiS7oREESQuBQICm\n/Q4uvPAclQORLLRw4cuYg51iHHfcWMtppL1pBEHa7IwzzuCf/3wac710rf4fiWSxAw44gM8+WwXU\n061bNzZu3Gg7knwDjSCIVXV1dfzzn89iDmGKEo1GbUcSkXa0dOlSzDojL19/Xc11111nO5K0ExUE\naZOioiLM1IKPGTN+ic/nsx1JRNpZbW0N4APC3HzzTNtxpJ2oIEirHXHEEZiT38L4fPDrX//adiQR\n6QCFhYVMnDgWMGuPevbsaTmRtAetQZBWqa+vp7CwELPLGrhuxG4gEelwPp+fWCwKFDN79m1ccskl\ntiPJDrQGQaww5QAgn1//errVLCJiRyhUDxQAIS699ArbcSTF2qUg3HPPPQwcOJD8/HyOOOIIXnvt\ntW/83oqKCjweT9LH/Pnz2yOapMCpp54KeDHDi9XMmDHDciIRscHn83HhhWdiRhIbdeZKlkn5FMPf\n/vY3pkyZwqxZsxg1ahR33303c+bM4dNPP6Vfv35J319RUcG+++7LSy+9xLBhw7bfXlZWRl5e3s5h\nNcWQFhzHj3lC8FFfX0l+fr7tSCJikXlO8AAun332IYMHD7YdSUjDKYY77riDCy+8kIsvvpjBgwfz\nxz/+kV69ejFr1qzd3q9z58507959+8f/LQeSHsyRr3HA5fjjD1c5EBEWLvx34jMvQ4YcaDWLpE5K\nC0IkEuG9995jwoQJO90+YcKExNni3+z000+nR48ejBo1innz5qUylqTIli1b2LixATO1EOLf//73\nnu4iIjlg1KhRBIMOZn8EP8cdd5ztSJICKS0ImzZtIhaLJd5lNuvevTuVlZW7vE9xcTG33347//jH\nP3jhhRcYN24cZ555Jg8//HAqo0kKdOnSBTN6AF99tcJuGBFJKw0NDZipxwZee+0D23EkBazvatOl\nSxemTZu2/evhw4ezefNmbr31Vs455xyLyWRHV1xxBWZhoovHE9JZCyKSZOLEb/PSS68Dcfr378+q\nVatsR5I2SGlB6Nq1K16vlw0bNux0+4YNG+jVq9de/5wjjzyS++67b5e/d8MNN2z/fOzYsYwdO7Y1\nUaWF7r77fzHvDuqJRmO244hIGnrxxRdxHHOF0+rVm6ivr6egoMB2rJxRXl5OeXl5yn5eyq9iOPro\noxk2bBj33nvv9tv2339/zjjjDG6++ea9+hnTpk3jmWee4csvv9w5rK5isOK0007j6afnA2EKCvKp\nq6uzHUlE0tTbb7/NUUeNBqJAXM/ZFrX1NTPlUwxXXXUVU6ZMYcSIERx77LHMnj2byspKLr30UgCm\nT5/O4sWLeeWVVwC4//778fv9HHrooXg8Hp555hnuuecebr311lRHk1Z6+umXE5/5VA5EZLdGjBgB\nhBNfBXjzzTc5+uijbUaSVkp5QfjhD3/I5s2buemmm1i/fj0HH3wwzz///PY9ECorK1mxonmBm+M4\n3HTTTaxcuRKv18vgwYOZM2cOkydPTnU0aQXz/60BKGLOnN1fqioiAmbBeteuvYAwxxxzjEYRMpTO\nYpDdcpxSoBpAj72I7DWPx4PrFgB1PPLII5x99tm2I+WctNsoSbJHaWkp5rLGEh588EHbcUQkg8Ri\nMcBMSU6ePMVuGGkVFQTZpXg8TnU1gAtUc+6551pOJCKZxHEcSkpKgCLA4a677rIdSVpIUwyyS/n5\n+YRCXiDC/PnPMX78eNuRRCQDmWlKM5qg5++OpSkGSbnGxkZCoaazMKIqByLSat27BzEvNXlMn66j\n4TOJRhAkSUlJCTU1cSDK3Ll/5vzzz7cdSUQymOMEMZc+Orhu3HacnNHW10wVBEniOAHMtsphXFe7\nJopI2/h8PmIxc0T8a6+9yMiRI21HygmaYpCUOuiggzA7oDXypz/daTuOiGQBc5CTC8QZNUonPWYK\njSDIThynKPGZFhSJSOo4jpP4zM9nn33I4MGDrebJBRpBkJT5/ve/jxk9qGPChAm244hIFtm2bRvg\nByIMGTLEdhzZCxpBkO2aG74P141azSIi2cc8xwQwhzhFbMfJehpBkJRYuHAh5h+uh4MOUrsXkdR7\n4YUXMAug4wwdOtR2HNkDjSAIsOPoQRGuW2M1i4hkL3PJYyPmkkeNVLYnjSBIm8XjcSAfcPD5Qrbj\niEgWu/jiczGjCDF+/vOf244ju6ERBKFTp06JBUQlfPTR64lLHUVE2ofZa8WsQdBzevvRRknSZo5T\ngLl6IaZdzkSk3XXu3JmtW7cCJUSjm/H5fLYjZSVNMUibPPbYY5gNTDyMG/dt23FEJAcsX74cyAPq\nGDZsmO048g00gpDjmhcnBnBdrT8QkY6hy6rbn0YQpI1KEr+GraYQkdxyzDHHYBYrOtxxxx2248gu\naAQhh51//vk88MDfAS/z5z+hY51FpEOZxYpRwNVzezvQIkVpteYhvjztaiYiHc48B3kBP65bbztO\n1tEUg7RKNBqlaXivT5/utuOISA6aO3cuZpF0lBtuuMFuGEmiEYQcNWjQIJYvXwV4qa7eSHFxse1I\nIpKDmkcytSdCqmmKQVrFcXxADNA/ShGxxxSEQszVDFW242QVTTFIizU0NGAOZgoSDAZtxxGRHDZv\n3jzMVVRRvvOd79iOIzvQCEIO6tevH2vWrAEK2bZtHSUlJXu8j4hIezG7uTpAvZ7jU0hTDNJi2qBE\nRNKJeU4yb1Rcd5vdMFlEUwzSCmbf8+LifMs5RETg2WefBWqAOC+//LLtOJKgEYQcM2bMGBYsWAAE\nWL++gp49e1pOJCICjhMEwng8HmKxmO04WUFTDNIiml4QkXTU/NxUgOvWWc2SLTTFIC0UAIqARttB\nRES2O/HEExOfOdTXa1fFdKCCkEO+/PLLxGcNiWOeRUTSw3PPPYd5AxPitNNOsx1H0BRDTvH5fMRi\nHsCL6zbYjiMishNt4JZabX3N9KUwi6Q5s/DHAUK2o4iI7EIM0NVV6UJTDDmlFK09EJF0NWXKFMzL\nUiMPPvig7Tg5TwUhR7z77ruY7Uz93HzzzbbjiIgkeeCBBzAD235+8pOf2I6T87QGIUd069aNTZs2\nASXaqUxE0pbj5GFKQkjP922kfRBkr+hIVRHJBOa5ygGKcN1q23EymvZBkL1UbDuAiMge/eEPf8As\nVHSpra21HSenqSDkAHO8cxQo5sorr7QdR0TkG02dOhVzNUOUMWPG2I6T0zTFkAN+/vOf87vf/Q7w\n4Lra41xE0pvjFAJmN0U957ee1iDIHnk8nsTjlo/ragtTEUlvZh1CEIjiuro0u7W0BkH2yHUdzJye\ndk8UkfRXXFwM+NGmSXapIOSEINBAaWmp7SAiInt0++2307Tlcl2dTna0RQUhy0WjUSAOwLXXXms3\njIjIXvh//+//YS51DHPxxRfbjpOztAYhy91yyy388pc3AeiMdRHJGI4TAKJ07dqFr7/+2nacjKRF\nirJb3bt3T/zjKtamIyKSMZo3d9PVV62lgiC71bwaOIzrxm3HERHZK+Z4egAfrqsTaFtDVzHIHgQT\nv6pYiUjmGDRoEOYlSi9TtuiRz3p+IERBQYHtICIie+2JJ57AHNrkY8WKFbbj5CQVhKxnjni+5JJL\nbAcREdlrBxxwAOAFapk8ebLtODlJaxCy2Lp16+jTpw+g7UpFJPM4ThFQR15eHpFIxHacjKM1CPKN\n7r333sRnml4QkUxUBxQl9nORjqYRhCzWfAZDQKuARSTjNF/qqFHQ1tAIgnyj5r8YOuxERDJZoe0A\nOUkFISdokxERyTw9e/bETJFqisEGFYSsFkDNW0Qy1dSpUzFvcJw9fau0A61ByGKOUwzEyMtr1Apg\nEck4W7dupXPnAUAjNTUbKSoqsh0po2gNguyGB2jgwAMPtB1ERKTFysrKMKMHDn/+859tx8k5KghZ\nLQr4+MUvfmE7iIhIK8WAOI899pjtIDlHUwxZatWqVQwYcDAA0ehmfD6f5UQiIi3nOCVADV6vl8ZG\nXZHVEppikF266667MJuMxFUORCSDRYECYjFdjdXRVBCy1KpVq4A8tPpXRDJbCF3maIfeWmapr7/+\nGnNQk1qsYzdhAAAgAElEQVS3iGQ6FQQbtAYhSzVvs1yM61bbjiMi0irabrn1tAZBdsn8pXCABttR\nRETaSFOlNqggZCnTul00iyQimU/PZTakvCDcc889DBw4kPz8fI444ghee+213X7/Rx99xJgxYygo\nKKBv377ceOONqY6Uk8wIgg+zDkFEJNOpIHS0lBaEv/3tb/zsZz/juuuuY8mSJRx77LGcdNJJrF69\nepffX11dzfjx4+nVqxfvvPMOd955J7fddht33HFHKmPlsEY0NCcimaywsBBzRZYWXHe0lC5SPOqo\nozj00EO59957t9+2//7784Mf/ICZM2cmff+sWbOYPn06GzZsIBAIAHDzzTcza9Ys1qxZkxxWixT3\nmpliKAGq9ZiJSMY68MAD+fTTT4EiXLfGdpyMkjaLFCORCO+99x4TJkzY6fYJEyawaNGiXd7njTfe\n4LjjjtteDpq+f926daxcuTJV0XKYLg0SkcxWXFyMOZlWB851tJQVhE2bNhGLxejRo8dOt3fv3p3K\nyspd3qeysjLp+5u+/qb7SEvoCgYRyWzvvfce4EVrEDqe1Ud8x+tb99YNN9yw/fOxY8cyduzY1AXK\nOh4gbjuEiEirmS2Wo0Cx7Shpr7y8nPLy8pT9vJQVhK5du+L1etmwYcNOt2/YsIFevXrt8j49e/ZM\nGiloun/Pnj13eZ8dC4LsicqBiGS2kpISqqqi6Plsz/7vm+YZM2a06eelbIrB7/dz+OGHM3/+/J1u\nf/nllzn22GN3eZ9jjjmGhQsXEg6Hd/r+Pn36MGDAgFRFExGRDBWNRjEvVVpT1dFSepnjVVddxdy5\nc/nrX//K0qVLmTp1KpWVlVx66aUATJ8+nRNOOGH790+ePJmCggIuuOACPvnkEx5//HF++9vfctVV\nV6UyVg4L7PlbRETSWF1dHVCDWYcgHSmlaxB++MMfsnnzZm666SbWr1/PwQcfzPPPP0+/fv0As/Bw\nxYoV27+/pKSEl19+mcsvv5wjjjiCzp0789///d9MmzYtlbFyWBjtgyAi2UEb/3Y0HdaUpcwCUC/g\n4LoamhORzBQIBIhEfEAU19Wlji2RNvsgSDpSmRKRzGb2QYihRYodTyMIWcq07kYgD9cN2Y4jItIq\nOu659TSCILsUDAYxjdtvO4qIiGQgFYQs1b17d0w50PoDERFpORWELDVkyBBUDkQk8/mBIrxeXebY\n0VQQstSIESMw/3t1maOIZLI8oDZx7LN0JBWELHXRRRdh1iA0JnYiExHJRHWAg9+v9VQdTQUhS/Xp\n0wfIB+L85S9/sR1HRKSVigGXbt262Q6Sc1QQsloeEODRRx+1HUREpJXiQJCZM2faDpJzdMB2VqsD\n4jttby0ikineffddzIZvDqeccortODlHIwhZzQvEWbdune0gIiIt9thjj2EWWvvw+fR+tqOpIGS1\nCGaaQUQk87z++utAA3qpskNbLWcxs0VpPhDBdRttxxERaRHzHJaHOahJz/0tpa2WZQ8aMCVBRCQT\n6TJtW1QQstjhhx9O0yVCIiKZR3sf2KSCkMXuuusuzJUMEInoHHURyTQ+oJAePXrYDpKTVBCy2NFH\nH42Zv6vjxhtvtB1HRKSFIkCE7373u7aD5CQVhJyQx+OPP247hIjIXps3bx4QA+B3v/ud3TA5Slcx\nZDnH8dC0BkGPnYhkimOOOYY333wbcHQVViu19TVTBSHLeb1e4vECoFaPnYhkDL/fTzTaCOThumHb\ncTKSLnOU3SotLQXqgVLbUURE9po5hdbFrEMQG1QQstzUqVOBABAlHo/bjiMispd8gIfi4mLbQXKW\nCkKWu/766zFnMtTzox/9yHYcEZG9ZI6rP/nkk20HyVkqCDnBrAR+6qmnLOcQEdmzG264AagFAtx/\n//2W0+QuLVLMAWY/c9MFXTdmN4yIyB6YBYpxzALFBttxMlZbXzN1fmbOaC4JIiLpzCxQFNv0ipED\nLrzwQsw0Qx6xmEYQRCTdFSV+1fOVTSoIOeCvf/0rZsFPTAsVRSStmautXKCQYcOG2Y6T07QGIUc4\nThAzzRDSYygiaeu0007j6adfAAI0Nlbh9XptR8pY2klR9opZqAhm0Y82HhGR9NT8XOXDdbUWoS20\nk6LslZKSksRnQas5RER2rxAzJarzF2xTQcgRH3zwAebo5zi///3vbccREfkGjUADPXv2tB0k52mK\nIYc4ToCmfc31OIpIupk6dSp//OMfgWLq6iopKCiwHSmjaQ2C7DUzt+cF/Lhuve04IiI7KSgooKEB\nIKYTHFNAaxBkr5mTHeOAV/shiEjaaWiIYp6jtJA6Hagg5JD169cDxUAjxxxzjO04IiLb1dbWYhYn\nhhkxYoTtOIKmGHKO4/gBc+mQHksRSRcDBgxg1apVmB1fQ3g8ev/aVlqDIC1i1iEUA/W4ri4jEpH0\n0Lz/gd68pIrWIEiLzJ8/HzO/5+Ohhx6yHUdEJCEIFFFWVmY7iCSoIOSY8ePHY/ZDCHPZZZfZjiMi\nwrXXXos5XDjKW2+9ZTuOJGiKIQeZobx8wIvr1tiOIyI5zjwnFQF1uG7cdpysoSkGabFRo0YlPguz\ncuVKq1lERCAA1GJOcZR0oYKQgxYuXIj5Xx/gsMMOsx1HRHLYL3/5SyAMBJg9e7btOLIDTTHkqOYV\nwwW4bp3VLCKSu8xzUSHQiOuGbMfJKppikFY54YQTMHN+IRYtWmQ7jojkLLP2wIwiSDpRQchRL7/8\nMmZL0wJGjhxpO46I5KBzzz0XiAGFuuw6DWmKIYc1TzN4tWmSiHS45uegIK7bYDVLNtIUg7TaNddc\ng1k9HODSSy+1HUdEcog5MK448ZXWHqQjjSDkOMcpxpzNENZjKyIdZvjw4bz//vtAgA8+eJtDDjnE\ndqSso7MYpE28Xi/xeBwo0qZJItJhHMeL2bBNVy+0F00xSJusW7cOs8VpLT179rQdR0RywIoVKzDT\nm3X07KmzF9KVRhAEx8kDzCJFPb4i0t7M4sQSzNbKWiDdXjSCIG12wgljE58FEpc/ioi0pyBQjbnE\nUdKVRhAEaFqsWIvX66WxUY1eRNpHYWEh9fVxIM7zzz/JSSedZDtS1tIiRUmJsrIyqqpqAaivryY/\nP99yIhHJRtp/peNoikFS4osvvsAsGmqkoKDAdhwRyUJXXnnl9s8vueRHFpPI3tAIgmznOB7MoSlh\nXDdiO46IZBlzaWMccHDduO04WU8jCJIyDzxwP2bRkMvJJ59sO46IZJFXXnkFM0qZD+iNXibQCILs\nRPODItIednxuCYfr8fv9VvPkAo0gSEqZ09UA8hJHQouItM3XX3+NOdYZPB5X5SBDaARBkjiOH3M+\ng0YRRKTtmkcPAmzZsp6yMu2e2BE0giAp161bJ8ALwP/8z//YDSMiGa22thZzamMACKscZBCNIMgu\nmSsa8oGormgQkVbLz88nFAoDeXz22YcMHjzYdqScoREEaRf77jsQqAcC3H///bbjiEgGqq+vJxTK\nw1y1EFE5yDAaQZBv1DxvqGuWRaTl+vbty9q1m4E469ZV0KtXL9uRckrajCCEw2F++tOf0q1bN4qK\nijjttNNYu3btbu8zd+5cPB7PTh9er5dIREPa6WD06NGYvyJ+Bg0aZDuOiGSQVatWsXZtNeChU6cC\nlYMMlLKC8LOf/YzHH3+cxx57jIULF1JdXc0pp5xCPL77d54FBQVs2LCByspKKisrWb9+vS6BSRP/\n+c9/MIsVwyxfvt52HBHJIPvuux9QA3ioqKiwnEZaw5eKH7Jt2zbuu+8+5s6dy7hx4wB48MEHGTBg\nAK+88goTJkz4xvs6jkO3bt1SEUPawf77D+SLL9YBUQYNGsSXX35pO5KIpLlFixYRi+UD9RQVQWlp\nqe1I0gopGUF49913iUajOxWBvn37csABB7Bo0aLd3rehoYF99tmHfv36MWnSJJYsWZKKSJIin3/+\nOdCAGUWosJxGRDLByJGjMKMHhaxfr9HHTJWSglBZWYnX66VLly473d6jRw82bNjwjfcbMmQIc+bM\n4emnn+bRRx8lGAwycuRIvUtNMxdffAEQBGIEg0GrWUQkvf3pT38CHKCAgQO7UFRUZDuStNJur2K4\n7rrrmDlz5m5/QHl5OWvWrOH8888nGo3u9Hvjxo1j//33Z9asWXsVJh6Pc9hhhzF27FjuvPPO5LC6\nisEac0VDHhBgzZrP6NOnj+1IIpKGmndiDeC6IdtxclpbXzN3uwZh2rRpnHfeebv9Af369aOxsZFY\nLMbmzZt3GkWorKxMrITfOx6Ph+HDh7Ns2bJv/J4bbrhh++djx45l7Nixe/3zpfU+/fRThg49CnDo\n27efLnsUkST77LMPphwUcsEFZ1hOk3vKy8spLy9P2c9LyT4I27Zto3v37sydO5ezzz4bgDVr1jBg\nwABefPFFxo8fv1c/x3VdDj/8cIYPH87//u//JofVCIJVpaWlVFfXAH4mT/4+Dz/8sO1IIpImIpEI\ngUBXmq5ccN2Y7Ug5Ly32QSgtLeXiiy/m6quv5l//+hfvv/8+U6ZMYdiwYTudCDhu3Dh++ctfbv96\nxowZzJ8/nxUrVrBkyRIuvvhiPvnkEy699NJUxJIUq6qqwuyn7uORRx63HUdE0kgwWIApB0X8+tfX\n2Y4jKZCSyxwB/vCHP+Dz+TjzzDNpaGjghBNO4KGHHtphNz5YsWIFAwYM2P71tm3b+PGPf0xlZSWl\npaUMHz6cBQsWcMQRR6QqlqSQ4zgceeTBLF68GCimrKyMrVu32o4lIpY98sgjuG4AqCcQiDJjxgzb\nkSQFtNWytFjzIqQ8Vq9eQd++fW1HEhGLHKcYiAEujY21eL1e25GENJlikNyyYsXnNJ302K9fP9tx\nRMQiM0pcC8ARRxykcpBFVBCkxQYOHEhpqR/z18fLGWdotbJILnrhhReAQqAYiCSmHyVbaIpBWs1x\nfJiSkMfatcvo3bu37Ugi0oHMc0AMKOLDDxdx8MEH244kO9AUg1jzwANzAD9Qr42TRHJM165dMS8h\nQfr376xykIU0giBt4vV6icf9QCOdOhXpqgaRHLBw4UJGj/420Aj4cN3onu4iFmgEQayKxWJAHGik\nqirKihUrbEcSkXY2evTExGde1qypsBlF2pEKgrTZ9ddPx0w1uOy332DbcUSkHZnDl8wuiRMnnqDp\nxSymKQZJibKyMqqq6oBGiouLqK6uth1JRFLst7/9Lddc07QJUljbKac5TTFIWjBrD1zAR01NDddd\np61WRbJJJBLhmmtuBrxAlNpavQnIdhpBkJSpqamhpKQUUxS8xGIRPB51UJFs0LxtfpArr/wxd955\np9U8smdtfc1UQZCUmjRpEs8++yxQAIQ0BCmSBcwUYhUQpLg4T1OIGUIFQdJOIBAgEokARXTu7Gfz\n5s22I4lIK333u9/lqaf+hTl/xcV1w7YjyV7SGgRJO+FwGDD7s2/ZEuKkk06yHUlEWmHt2rU89dTT\nmLMWfGzYsNp2JOlAKgjSLtavX5f4rJ4XXyynoqLCZhwRaQVzUqtZU3TqqePo3r277UjSgTTFIO1m\n4sSJzJ+/EIgA4LqNdgOJyF7z+wNEoy4QoFu3fDZu3Gg7krSQ1iBIWuvTpw/r1tUANQD6/yeSAYLB\nIOFwAKgGvCr3GUprECStrV27Fo+nDjAnP3br1s12JBHZjauvvppwOIopBx5isYjtSGKJCoK0O3Ne\nA4DDpk1bOeOMM6zmEZFdW7BgAbfddjeQD+SzZMl72sskh+n/vHSILVs2Ys5r8PHPf77ILbfcYjuS\niOwgFosxZsxEoB6IMXv27xk2bJjtWGKR1iBIh3n55ZeZMOF0zCVTQV555VnGjRtnO5ZIzovFYvh8\nJZjLk2HcuKN55ZVX7IaSNtMaBMkY48eP5+c/vwwIAiFOOOE7bNq0yXYskZzn8wUwx7bXMWHCSJUD\nATSCIBYccsghfPTRciCMWQQV0jyniCWO4wUKgQjBoIeGhnrbkSRFdJmjZKTevXuzfn01UAf4tX2r\niAV+v59oNICZ9vPo7JQsoykGyUjr1q2jqKjpdLgAjpNnNY9IrgkGg0SjQUw5QOVAkmgEQazy+fKI\nxRppWpeg/78i7c/r9RKPBzFXLEA0GsXn89kNJSmnEQTJaI2NUfz+AOakuJLEfKiItJcePXoQj7uY\np/88QqGQyoHskv5WiHXhcAjHcTA7txXjOB5cN247lkjW6dKlC1u21GNG7ODrr9cRCATshpK0pSkG\nSRumJIDZUCmm/d9FUsjn8xGLeTCjdR6qqrZQWlpqO5a0I00xSNZo/ovsA7w7FAYRaQuPx0MsFseM\nHOTz9dcbVA5kjzTFIGnFdd1EMSgAGnEcH7FYRPskiLSSx+PBdc2/J4izdes6OnXqZDuWZAA960ra\nMSWheV8Er9dHKBSymEgkM5kh5gBmv5EwW7euUTmQvaaCIGkpHm+kpKQIsx6hgPz8zrz++uu2Y4lk\nDMfxYU5ljAAO4XBY5UBaRAVB0ta2bdvw+2OYa7XzGDXqOzz33HO2Y4mktXA4jOMUYKbpGoA4rhvH\n7/dbTiaZRgVB0lo4HKa4uAhzbkOYU075HhMmTLAdSyQtvfXWWwSDnTGnMtYA6MovaTUVBEl71dXV\nTJ78fSAGRHn55XK6detmO5ZIWrnwwgs5+uhxNO2OWFpaqnIgbaKCIBnh4Ycf5qqrfor5Kxtl06Ya\n+vfvbzuWSFo488wzmTv3Icx6gyIOPHAgVVVVtmNJhtNGSZJRFi1axMiRI3e4xUsk0kBeng57ktxk\nLgv2Y8pznNNPP4V58+ZZTiXpQMc9S84Jh8MUFhYTizmYd0xennnmSU455RTb0UQ6TDgcJhgsBAox\n25QHWbXqC/r162c5maQLFQTJWWbr2KYNYDwcdtj+vPfee7ZjibS7l156iRNP/B5m1KAOMIVBVyrI\njlQQJKcddthhLFnyJWbVdiPQoL8jktWGDh3K0qWraVq0W1JSwLZt22zHkjSksxgkp73//vs8/vgD\nmCfLBiAfx3EoLy+3G0wkxaLRKI7jYenSZUAt4GHKlLNVDqTdaARBskJDQwMFBUWYztsIFDN0aD8+\n+eQTy8lE2u6ZZ57h1FNPTXxVCDTwzjtvc/jhh9uMJWlOUwwiOygqKqKuzovZWMkBQsTjcZ0MKRlr\n8ODBfPHFeszGR/k4Tph4PGY7lmQATTGI7KC2tpbbb78e8AIhoBCPJ4+XX37ZcjKRllm+fDmO4/DF\nF19gykGAE04YqXIgHUYjCJKVGhsbycsL0DzlUEJhYYza2lrLyUT2bN999+Wrr75KfFUCNLBmzVf0\n6dPHZizJMBpBENkFn8+H68b461/vBYqAGurq6nAcL4899pjteCK7tGXLFhzHw1dfrU/c4iM/P4rr\nRlQOpMNpBEGy3tatW+ncuRsQB1zMKXf1+rskaWXkyJEsWvQ+ZnosCsS4554/ctlll1lOJplKIwgi\ne1BWVobrNrLffvtiRhMigPnH89///d9Ws4m8+eabOI6PRYvewFyq20CnTvm4blTlQKzSCILklM2b\nN9O1a1cgH/NkXATUsnr1avr27Ws3nOScsrIyqqrqMaU1AER59NGHOeussywnk2ygEQSRFujSpQuu\n63LzzddhFn+Zfzz9+u3HkCFDrGaT3HH22WfjOH6qquow5SBIWVkBrhtTOZC0oREEyWlmf4Sm0QSz\nNuG73/0uTzzxhN1gkpXKy8s5/vhv01RMzahBhM2bN9G5c2eLySQbaQRBpA1c1+WDD97EbKoUBuDJ\nJ5/CcRwef/xxq9kkezQ0NJCfn8/xx4/HlAM/4HDmmd/FdeMqB5KWNIIgkjBq1Chef30xzZsseYFG\n1q5dS+/eve2Gk4zl8Xhw3QLMqYtBwMOgQb1ZtmyZ5WSS7bTVskiKBYNBwmEPZtohCOTRtI9CQUGB\n3XCSMcwCxG00X1rrkJcXpapqq/4eSYfQFINIioVCIcLhqsRXDmabW4fCwlICgQANDQ0W00m6GzJk\nCI5TQFVVCFMOioEQixa9TCQSVjmQjKGCILILfr8f13VZvHhB4pYA4CcScSko6I7H4yEcDtuMKGmm\ntLQUxynm888rMKNPjYDDzJnTcd0YxxxzjN2AIi2kKQaRvbBmzRr69etP07bNRgFeb5jq6mq9K8xR\nsViMrl27UlXVgFnkmk/T+pXf/ObX/OpXv7IbUHKaphhEOkDfvn1x3Th/+MONiVuCQD2xmJfCwm44\njqNLI3NITU0NjuPg8xVQVVWFKQdBIMqpp07CdaMqB5LxVBBEWmDq1Km4rsvnn3+QuMWDWadQwOmn\n/xDHcfjRj35kMaG0pzlz5uA4DiUlvTAbbUWAQgBuvvlXuG6Up556ymZEkZTRFINIG7iuS3FxMXV1\nTZdGupiDdrw4TpytW7dSWlpqN6S0WZ8+fVi3bgNmLUoDphgGgTqee+45Tj75ZKv5RHZFUwwiFjmO\nQ21tLa67jbPP/j7mxMggEMN1i+jUqTeO4/C9733PclJpqQcffBDHcXCcAOvWbQFimPJXQFlZCeHw\nFlzXVTmQrKWCIJIijzzyCK7byI03Xpu4pZamf2JPPvkCjlOK4zg8++yz1jLK7q1evTpRChzOO+/i\nxK1N00gORx01HNetZcuWLfj9fotJRdqfphhE2kldXR29evWipiaMGZquwbzY5AP1gMuyZcsYNGiQ\nzZg5LxqNUlhYSDTatKFRmOZFhy6OE+HVV19lzJgxVnOKtJSmGETSVGFhIdXV1bhumF/96meJWx3M\nlrvmxehb3zoMxylkwIABLFmyxF7YHFNbW8vgwYNxHC9+fynRaABTCKox0wgevvOdcbhuiHg8rnIg\nOUkjCCIdbOjQoSxduhTzbtVL874KRZh3rlHuuusuLr/8clsRs1J5eTmnnHIKdXVNIzphzGZGZhMs\nqMHn81FXV6fpA8kKGkEQyTCffvopruviunXcfvsNBAJBTDmox7xYwRVXXIHjFOM4Dl26dOGrr76y\nGTkjua7LCSeckFhTUMLxxzeVg0bMKE4eAP3796C6ei2u6xKNRlUORBJSVhD+/Oc/c/zxx9OpUyc8\nHg+rVq3aq/vNmzePoUOHEgwGOfDAA3nyySdTFUkk7V111VWEQg24bg0TJ47HLGwswKxTqAV8bNmy\nhX33PQjHKaVTpzJGjBjBunXrrOZOR7FYjJ/+9Kd4PB4cpxCPp5h//etfmH0KzHka5nGFQw89lFBo\nK67rsnLlSoqLiy0mF0lPKSsIDQ0NnHjiicyYMWOv7/PGG29w1llnMWXKFD744APOOecczjjjDN5+\n++1UxRLJGC+++OL2kYWPP16cuLURsyGPA1SzbVstixcvpk+fwThOAMfxkZ+fz7Rp03Ju+u2xxx4j\nLy8vMUJQiM9XzF133ZV4HBox+1IUYy5PhGuuuQLXrcZ1Xd5//30CgYDF9CLpL+VrEN555x1GjBhB\nRUUF/fv33+33nnnmmVRVVfHSSy9tv238+PF069aNRx55JDms1iBIjpo6dSqvvfYa7733CabXN2De\nGddh3hU3JH6NAxE8HoeRI0fy1FNPUVZWZi13KsTjcS666CKeeuopQqEQoVAjzX/mRsCX+BXM6Esc\nCDFx4kSee+45vF6vneAilmX0GoQ333yTCRMm7HTbhAkTWLRokaVEIunpzjvv5N1338V1Q7huPZdd\ndhk+X9Npkl7MO+WmA4PyiMfjLFz4Jp07907sv+DBcXyJ8wN8/OQnP9nracCOUFVVxVlnnUXXrl23\n70PgOH4cpxSvt4j7759HVVUVoVAIUwCapgxKgEaCwSDXXXcdsVgNrtuA67q8+OKLKgcibeCz+R+v\nrKykR48eO93Wo0cPKisrLSUSyQz33HMP99xzz/avw+Ewp556KvPnz6f5fIBo4iOEWZDnADFisRiz\nZs1m1qx7Me+4PZgFkp7EfY3evXsTjUbxeDyMGzeOI488knHjxlFUVET37ubI6/x8M6ff2NhIOBxm\n3bp1fPLJJ1RVVfHoo49SVFTExx9/TGVlJdXV1Tv8CZo2HyqieUOpfJqv6IDmyw6DiT9HCVCH3++l\nU6dOfPbZZxk/OiKSznZbEK677jpmzpy52x9QXl7O6NGjUxpKRFomEAjsNFXX5N133+WSSy7h888/\np7bWLHo0pcBs1GQ+ajEv1pHEr4VALevWVdK0wO+RRx7hkUeeBqZjCkcAM7fvYC4RDCW+zsO8oIcT\nPy+Y+L2mwco8zIt90/23JX7PR/MiQhdowOutZ9iw4Vx22WU6AEvEgt0WhGnTpnHeeeft9gf069ev\n1f/xnj17Jo0WbNiwgZ49e37jfW644Ybtn48dO5axY8e2+r8vku0OP/xw3nnnnV3+3pYtW7j66qtZ\nv349r7/+Otu2baO5MAQxQ/ku5p27S9N2w02XYprP2eFzL6YAeDClwY8pCSU0X53RYL7bCeHxeCgs\nLGHy5MnccsstdOrUKXV/cJEcVF5eTnl5ecp+ntVFimeddRZbt27d6Z3PhAkT6NatGw8//HByWC1S\nFBER2Sttfc1M2RqEyspKKisr+eKLLwD45JNP2LJlCwMGDNg+Tzhu3DiOOuqo7dMWU6dOZfTo0fz2\nt7/ltNNO44knnqC8vJzXX389VbFERESkFVJ2FcPs2bMZPnw45557Lo7j8J3vfIfDDz+cZ555Zvv3\nrFixYqcphWOOOYbHHnuMuXPnMmzYMB566CH+/ve/c+SRR6YqloiIiLSCzmIQERHJQhm9D4KIiIik\nJxUEERERSaKCICIiIklUEERERCSJCoKIiIgkUUEQERGRJCoIIiIikkQFQURERJKoIIiIiEgSFQQR\nERFJooIgIiIiSVQQREREJIkKgoiIiCRRQRAREZEkKggiIiKSRAVBREREkqggiIiISBIVBBEREUmi\ngiAiIiJJVBBEREQkiQqCiIiIJFFBEBERkSQqCCIiIpJEBUFERESSqCCIiIhIEhUEERERSaKCICIi\nIrY5+qAAAAl+SURBVElUEERERCSJCoKIiIgkUUEQERGRJCoIIiIikkQFQURERJKoIIiIiEgSFQQR\nERFJooIgIiIiSVQQREREJIkKgoiIiCRRQRAREZEkKggiIiKSRAVBREREkqggiIiISBIVBBEREUmi\ngiAiIiJJVBBEREQkiQqCiIiIJFFBEBERkSQqCCIiIpJEBUFERESSqCCIiIhIEhUEERERSaKCICIi\nIklUEERERCSJCoKIiIgkUUEQERGRJCoIIiIikkQFQURERJKoIIiIiEgSFQQRERFJooIgIiIiSVQQ\nREREJIkKgoiIiCRRQRAREZEkKggiIiKSRAVBREREkqggiIiISBIVBBEREUmSsoLw5z//meOPP55O\nnTrh8XhYtWrVHu8zd+5cPB7PTh9er5dIJJKqWCIiItIKKSsIDQ0NnHjiicyYMaNF9ysoKGDDhg1U\nVlZSWVnJ+vXr8fv9qYolIiIirZCygjB16lR+8YtfMHLkyBbdz3EcunXrRvfu3bd/SOqUl5fbjpBR\n9Hi1nB6zltNj1nJ6zDqe9TUIDQ0N7LPPPvTr149JkyaxZMkS25Gyiv5RtYwer5bTY9ZyesxaTo9Z\nx7NaEIYMGcKcOXN4+umnefTRRwkGg4wcOZIvv/zSZiwREZGct9uCcN111yUtIvy/HwsWLGj1f/zo\no49mypQpHHLIIYwaNYq//e1vDBo0iD/96U+t/pkiIiLSdo7ruu43/ebmzZvZvHnzbn9Av379yM/P\n3/71O++8w4gRI6ioqKB///4tDnThhReyYcMGnn/++aTfGzRoEMuXL2/xzxQREck1++23X5tG5H27\n+80uXbrQpcv/b+duQ5oK2ziA/8/M+dI0wnK+LFQoqAVB8yW0speVFYVBYSm9kF+kPhgKQQVBaSGR\nJNmHERHVPkySsoKg0AWBhoqaU0yTohfTcoOSlIKs9Ho+uecZm2vnPNs5y12/L+K9++xc5+Kv3Gzn\n3HGS31wsIkJvby8MBoPH1/mrB8YYY0weXhcIYsw8pvj69WsAQH9/P8bGxpCSkoKFCxcCAIxGI9as\nWYOqqioAQEVFBbKzs7F06VJMTEzg6tWr6O/vx/Xr1/1VFmOMMcYk8NtNiteuXYPBYMDBgwchCAJ2\n7tyJ9PR0PHr0yDnn3bt3sNvtzt/Hx8dRUlICvV6Pbdu2YXR0FM3NzcjIyPBXWYwxxhiTwOs9CIwx\nxhgLTYrvgzAb3rpZPCk9A4CGhgbo9XpERkZi5cqVePjwYYArDR6Tk5MoLS3F4sWLodFosHv3bnz6\n9MnrMaGWM5PJhLS0NERFRSEjIwPPnz/3Or+vrw8bNmxAdHQ0dDodzp8/L1OlwUNMzz58+ODxCbGm\npiYZK1ZOc3Mz8vPzodPpoFKpYDab/3pMqGdMbM+kZixoFwi8dbN4UnrW1taGwsJCHDp0CL29vThw\n4AAKCgrQ0dERwEqDR1lZGe7fv487d+6gpaUFExMT2LVrF6anp70eFyo5q6+vR1lZGc6cOYOenh7k\n5ORgx44dGB4e9jh/YmICW7duRWJiIrq6ulBbW4vq6mrU1NTIXLlyxPZsRmNjozNPdrsdmzZtkqli\nZf348QOrVq1CbW0toqKiIAiC1/mcMfE9myE6YxTkOjs7SRAEGhoa+uvcW7dukUajkaGq4CamZ/v2\n7aO8vDyXsS1btlBRUVGgygsa3759I7VaTXV1dc6x4eFhUqlU1NjYOOtxoZSzrKwsKikpcRlbtmwZ\nnT592uN8k8lECxYsoJ8/fzrHLly4QMnJyQGtM5iI7dn79+9JEATq6uqSo7ygptFoyGw2e53DGXPl\nS8+kZixoP0GQirduFqe9vR15eXkuY3l5eWhtbVWoIvm8ePECv3//drl+nU6HFStW/PX6QyFnv379\nQnd3t6h8tLW1Yf369YiIiHCZ//nzZwwNDQW03mAgpWcz9uzZA61Wi3Xr1qGhoSGQZf7TQj1j/w+x\nGZtTCwTeulk8u90OrVbrMqbVal2eNpmr7HY7wsLC3Pb60Gq1cDgcsx4XKjn78uULpqam3PIRHx8/\naz5my9PMa3OdlJ7FxMTg8uXLuHv3Lp48eQKj0Yj9+/fDYrHIUfI/J9QzJoXUjMm6QOCtm8ULdM/m\nIs6Zcnz9LpT9V1xcHMrLy5GVlQWDwYCKigocPXoUly5dUrq0oMQZE09qxvy2UZIvysvLcfjwYa9z\nlixZ4rfzqVQqGAwGvHnzxm/vKbdA9ywhIcFt1e1wOJCQkCD5PZXma8/+/PmDqakpfP361eVTBLvd\njtzcXJ/PNxdy5smiRYsQFhbm9mmKw+FAYmKix2Nmy9PMa3OdlJ55kpmZiZs3b/q7vDkh1DPmL75k\nTNYFQrBt3fwvCHTPsrOzYbVaceLECeeY1WrF2rVrA3bOQPO1Z+np6QgPD0dTUxOKiooAACMjIxgc\nHEROTo7P55sLOfNErVYjPT0dTU1N2Lt3r3PcarWioKDA4zHZ2dk4efIkJicnnd8RW61WJCcnIyUl\nRZa6lSSlZ5709PQgKSkpECX+80I9Y/7iU8ak3zsZWKOjo2Sz2chisZAgCPT48WOy2Ww0NjbmnLN5\n82aXO4PPnTtHjY2N9PbtW7LZbFRcXExqtZo6OzuVuATZSelZa2srzZs3jy5evEivXr2iqqoqCg8P\np46ODiUuQXbHjh0jnU5HT58+pe7ubtq4cSOtXr2apqennXNCOWf19fWkVqvpxo0bNDAwQMePH6eY\nmBj6+PEjERGdOnWKjEajc/74+DglJCRQYWEhvXz5khoaGig2NpZqamqUugTZie3Z7du3qa6ujgYG\nBmhwcJCqq6tJrVbTlStXlLoEWX3//p1sNhvZbDaKjo6myspKstlsnDEvxPZMasaCdoFw9uxZEgSB\nBEEglUrl/Pm/j3OkpqZScXGx8/fy8nJKSUmhiIgIio+Pp+3bt1N7e7sS5StCSs+IiO7du0fLly8n\ntVpNer2eHjx4IHfpipmcnKTS0lKKi4uj6Ohoys/Pp5GREZc5oZ4zk8lEqampFBERQRkZGdTS0uJ8\n7ciRI5SWluYyv6+vj3JzcykyMpKSkpKosrJS7pIVJ6ZnZrOZ9Ho9zZ8/n2JjYykzM5MsFosSZSvi\n2bNnbv+3BEFw/s1xxtyJ7ZnUjPFWy4wxxhhzM6cec2SMMcaYf/ACgTHGGGNueIHAGGOMMTe8QGCM\nMcaYG14gMMYYY8wNLxAYY4wx5oYXCIwxxhhzwwsExhhjjLnhBQJjjDHG3PwHfe2l5b7GofwAAAAA\nSUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x9ffd92c>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compute the covariance matrix."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.cov(np.array([x, y]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "array([[ 0.50309341,  0.00297389],\n",
        "       [ 0.00297389,  0.49683412]])"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Problem 2"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$X$, $Y$, and $Z$ are distributed uniformly in the ranges $[0, \\alpha], [0, \\beta], [0, \\gamma]$.\n",
      "This means they have variance $\\alpha^2 / 12, \\beta^2 / 12, \\gamma^2 / 12$ which give us the diagonal terms in the covariance matrix.\n",
      "The other terms are calculated as follows (shown for $X$ and $Y$, but similar for other pairs of random variables):\n",
      "\n",
      "$$\n",
      "cov(X, Y) = \\langle (X - \\overline{X}) (Y - \\overline{Y}) \\rangle = \n",
      "\\langle (\\alpha \\lambda - \\alpha / 2) (\\beta \\lambda - \\beta / 2) \\rangle =\n",
      "\\langle \\alpha (\\lambda - 1 / 2) \\beta (\\lambda - 1 / 2) \\rangle =\n",
      "\\alpha \\beta \\langle (\\lambda - 1 / 2) (\\lambda - 1 / 2) \\rangle =\n",
      "\\alpha \\beta * var(\\lambda) =\n",
      "\\frac{\\alpha \\beta}{12}\n",
      "$$\n",
      "\n",
      "Thus the covariance matrix is:\n",
      "\n",
      "$$\n",
      "C = \\: \\begin{bmatrix}\n",
      "\\alpha^2 & \\alpha \\beta & \\alpha \\gamma\\\\ \n",
      "\\beta \\alpha & \\beta^2 & \\beta \\gamma \\\\ \n",
      "\\gamma \\alpha & \\gamma \\beta & \\gamma^2\n",
      "\\end{bmatrix}\n",
      "$$\n",
      "\n",
      "The correlation matrix is calculated as:\n",
      "\n",
      "$$\n",
      "r_{ij} = \\frac{C_{i j}}{\\sqrt{C_{i i} C_{j j}}}\n",
      "$$\n",
      "\n",
      "which results in:\n",
      "\n",
      "$$\n",
      "r = \\: \\begin{bmatrix}\n",
      "1 & 1 & 1 \\\\\n",
      "1 & 1 & 1 \\\\\n",
      "1 & 1 & 1 \\\\\n",
      "\\end{bmatrix}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For $\\alpha = 1, \\beta = 2, \\gamma = 3$, compute the expected covariance and correlation matrices as described. Then generate points with these parameters, compute the actual matrices, and verify they match expectations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def cov_and_corr(a, b, c, num_samples=1000000):\n",
      "    \"\"\"Generate points as the homework specifies. Return covariance and correlation.\"\"\"\n",
      "    abc = np.tile(np.array([[a, b, c]]), (num_samples, 1))\n",
      "    lambdas = np.tile(np.array([np.random.random(num_samples)]).T, (1, 3))\n",
      "    xyz = abc * lambdas\n",
      "    cov = np.cov(xyz, rowvar=False)\n",
      "    corr = np.corrcoef(xyz, rowvar=False)\n",
      "    return cov, corr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a, b, c = 1, 2, 3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "abc = np.array([[a, b, c]])\n",
      "cov_expected = np.dot(abc.T, abc) / 12.0\n",
      "corr_expected = np.ones((3, 3), dtype=float)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cov, corr = cov_and_corr(a, b, c)\n",
      "print 'cov:\\n', cov\n",
      "print 'corr:\\n', corr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "cov:\n",
        "[[ 0.08329172  0.16658343  0.24987515]\n",
        " [ 0.16658343  0.33316686  0.49975029]\n",
        " [ 0.24987515  0.49975029  0.74962544]]\n",
        "corr:\n",
        "[[ 1.  1.  1.]\n",
        " [ 1.  1.  1.]\n",
        " [ 1.  1.  1.]]\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.allclose(cov, cov_expected, rtol=0.01)\n",
      "print np.allclose(corr, corr_expected, rtol=0.01)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n",
        "True\n"
       ]
      }
     ],
     "prompt_number": 8
    }
   ],
   "metadata": {}
  }
 ]
}