tobydriscoll/TB-Lecture-04.ipynb

## TB-Lecture-04.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 4: Singular value decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometric interpretation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SVD has many possible interpretations and uses. We can visualize the geometry of the SVD using the same 2D picture used to illustrate the induced 2-norm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAIhCAYAAADn6rozAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd8VfXh//HXuSN7MkLYSxEBGYqi4sA90Vq3VovaqtVaW61a/ba/apet2oGz1tGqdVF3Fa0DRREFGSoOhoKMQMII2eOu8/vj5K7IJvd+7ng/H488cs8h0TcRwzuf8xmWbds2IiIiIiIJ5jIdQERERESyg4qniIiIiCSFiqeIiIiIJIWKp4iIiIgkhYqniIiIiCSFiqeIiIiIJIWKp4iIiIgkhYqniIiIiCSFiqeIiIiIJIWKp4iIiIgkhYqniIiIcOmll3L77bebjiEZTsVTREQky7377rs8+OCDrFmzxnQUyXAqniIiIlksGAxy4403YlkW69atMx1HMpyKp4iISBabOnUql1xyCbZtq3hKwql4ioiIZKmqqireeustLr74YkpLS6murjYdSTKcx3QAka4wd+5cnnrqKSzLYuXKlTzwwAPcf//91NXVUVVVxW9+8xsGDx5sOqaISEq5/vrrufXWWwHo1avXFud46vurdCUVT0l7y5Yt4/HHH2fq1KkAXHTRRRx44IE8+uijBAIBDjvsMPbdd19+9rOfGU4qIpI63nzzTcrLyxk9ejTgFM9ly5bR2NhIcXExoO+v0vVUPCXtTZ06NW4LkKamJrp168aECRNYs2YN1157LVOmTDEXUEQkxfj9fm655RZefvnlyL3KykoA1q1bFyme+v4qXU1zPCXt3XDDDeTn50euP/jgA44++mgA+vXrx2233UZ5ebmpeCIiKeeOO+7gggsuoLS0NHKvV69eAHHzPPX9VbqaRjwl7fXv3z/y+ssvv2Tt2rUcccQRBhOJiKSulStXctdddzF69Gief/55LMvCtm2WL18OELeyXd9fpaupeEpGefPNN8nNzWXixImReytWrNDEdxGRDtdeey0vvPACBxxwQNz9Rx99lClTpmx1SyV9f5WuoEftktba2tq44YYb+PzzzwHnG+OYMWPIzc0FwLZt7rjjDpMRRURSxvTp0xk6dOi3Sic4j84hOuKp76+SCBrxlLQ2ffp07rjjDvbbbz/cbjcrVqygrKws8uu///3vufDCCw0mFBFJDZ9//jlXXnklixcv3uKvV1RUANHi+corr+j7q3Q5y7Zt23QIkV21adMmbrjhBrp37w7AzTffzBVXXEFeXh45OTmccsopHHXUUYZTioiYs3jxYqZMmcKCBQsIBoOMGTOGt956K7IoKBQKMXnyZBYtWkRVVRX5+fkccMABXHHFFbz22mv6/ipdSsVTRERERJJCczxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKVQ8RURERCQpVDxFREREJClUPEVEREQkKTymA4iIiIgANFNPHkV8xUJaaSSfUhYzl74MYwVf4MFLT/pRzUpGMoGNVNObgfRmIEGClFBm+rcg26HiKSIiIknVwEYW8ApuvMzlZTaymgYaWMGnuMmnlVZCQBCLEDYBwI75fBsIYGFjE8IiCISwGcIo3HgZy8EMZE8GsgdHcpKR36NsmWXbtr39DxMRERHZeS3U8zGv4Kedd3iAdSylFR+tNOLHKZHBjjcb8Hd8XgAIdbwFOv0zt/XxdszHBwE3ueSTz2lcyACGcBKn0Y8BCfrdyvaoeIqIiEiXCeLHTzvP8AvW8BmrWUI91TG/7rwB+Drd21qR9BNvSx8fW2I7f174nwWQSx4DGMShHMWPuIYKKimgYHd/27KDVDxFRERkt9RRRS0rmc8zzOBOQrjw49/iaGXsvfhRSufRObiw8eKnnUqG0cBmutOPJppw4cZDDnWsx0Mha1iOlwKaaI6U0VDM+/C/xxcfITIqagMllPITrmd/DmQSR3b1l0Y6UfEUERGRndZGI4t5lfV8xWv8lnbaIkUvPNK4peIJYJGLmzzKqGQ0JwAuDuAUgoToSX96MZggfnLI22YGPz685LCURfjxsZwlLONzvmE5H/EeDbRQx+bIaGhY7IhpuIDawKEcwamcxgVcSCmlu/7Fka1S8RQREZEdtpyZVPMZ73EX61kSKXWxJTM84ui8tvBQQAk9OZaryaWAURxDBYMTntWPn7d4hRqqeYpHWc4yatkYNzoaLZ5uwMIG+tCXC5nCFL7PkCTkzCYqniIiIrJNbdTzFa+xglnM5u64khkunp1HD/fiaCoYxjhOYRTHmQm+BW8wndnM5Bmm8Q3fxGR2NvoJ4cJ53A/FlPAjLuMCzmNvhhtMnTlUPEVERGSL6ljFRpYwnavYtJXRzfC2RzY2AzmQo7ieEioZykHGcu+INtr4gs94hmncw1TaCXXMN3VGPsMFFCCfAn7FjZzDdxmoFfG7RcVTRERE4gTx8SUv8hznE4xZUx47VzKIM1ezkB6czT/pzWiKqcDCMhF5t9jYzGQm1/NzvmAZ7fgIRR69W9i4AIsSSniM+ziGI8gn33TstKTiKSIiIoBTwF7mEj7lXwRwEYpsfBT+dbDJwUsRozid7/B3XBl2+vZSlvFTrmE2H9FIc6SAxtqDIbzHdCrpZSZkGlPxFBERyXItbOBtrqWahVTzGcC3VoK7yaWIXpzO4wzkECM5k2k1a7iMn/Aec2iipdOveiigkDM5hQf4M168RjKmIxVPERGRLNXEWlbyGvO5lxrmx23uHjt3cwhHcz4v4yHXYFozVrCSs7iIhXxKsOMsJOcNwOJQDuKX/JRjOdxgyvSh4ikiIpKFGlnDU4ynhZrIBuuxK9MBjuQ2RjOFQnqaCZlC5rGQ4zmDTTRAZOFRtIA+wt84n+/gjtyTLVHxFBERySLrmcebnE89q/DRBsQfKekml0lMpQ8TqGSssZypqJlm7uNfXM9vOxYcOYuOwANYHMAY3uQxiikyGzSFqXiKiIhkgQa+YjXTmc9tNFMVtyWSDRTQj6F8l705jz5MMJg09c1iDjdyK7OYh1M6wwus3OzPGG7jOibpa7hFKp4iIiIZrp5lvMT++KiPO7c8ALgppJQ9OJ4n6M4IUxHTTjMt/Jzf8m9eoolWYud+evAymyfYn32MZkxFKp4iIiIZqp7FzOFyallMMzVA/GP1YgZxOh+RTw9jGdPd67zHaVxOCwGcx+7Oo3cLi58zhdu4xnDC1KLiKSIikmGC+GhkGTM5gwYWx61WByhnHKO4hgGcQB7dTcXMGBvYxOn8hPeYD3iJ7vvp4vdcxfV8H0/HkZzZTsVTREQkg/io5w2OoJaFcfcDgJcelDOSw3iEYgaaCZih6mjgF/yF+3mu40507ueJTORl/paWpzp1NRVPERGRDGBjs5r/8A3TWMmz3/r1PCo4iY8poLeBdNnjRqbyRx6GTpvKX8Z3uZ2rKabQTLAUoeIpIiKSAeZyKct5gCDROZwAhQxgGFcwiHMo0ihnUrzCe5zKtQTj/kvAZA7lJf5qKFVqUPEUERFJYzW8ThXPsIyHsQli48zntAEvxRzOC/TmSMMps89jvMKPuY0GmjvuuHHh4Xscz738nELyjeYzRcVTREQkTW1iNu9yGDbBuG2SbOBgnqY3x5FDqal4We8TljKeCwhgQ8ziop9yNn/lanPBDHJt/0NEREQklbSznk+4mIVcjN2xXt2Dc+pQAQMYx58ZyFkqnYaNYRj/5a8Mom/HHWerpbt5iauYik32jf1pxFNERCTNzOZwank37vQhgH6cywE8YSqWbMV0ZjOZ6wlFjtl0PM2vOYsjzAUzQCOeIiIiaWItjzKPo9jMbMD5S9wNlDKaIVzJOO43mk+27EQO5m3uoozijjsuwMv3+TOP85bJaEmnEU8REZE0UMtM5jEJAD9EHtK6yOcwPqGIPU1Fkx10G09yA/cTu8l8Dl5qeJoyioxmSxaNeIqIiKQwP3Us5UqWcmXkngfIoZgBXMpBzFDpTBPXcy4Pcj3Rk43c+IBJ3Mgq1htMljwa8RQREUlhi/gOm3jxW/M5+3ARo3jYVCzZRUGCnMgveZ2FxK50P42DeI7/MxcsSTTiKSIikoKaWMhyrqSOGYDzF7YH6MZhDON2Rmg+Z1py4+YVfsfpHBpz18OrfMpfeMlYrmTRiKeIiEiKaeMbPmYfQjQRIHoSkYWXffmQYvY1GU+6wOes5GCuowEfseOAb3ELR7KPuWAJphFPERGRFFLPTNbwW0I0Ac6qdTduKrmE0bym0pkhRjKQT7gr5gQjZ4/P3/Icm2g0GS2hNOIpIiKSImp4kK/54bfuF3EAo5ljIJEk2o95gHt4ldiV7pMYwdv8P6O5EkUjniIiIimgnRVU8/e4e7kMoSffYzjPG0oliXYXP+ByTiC60h3eZwlLWGsuVAKpeIqIiBi2kitZxBBamB93vw8/Y08eI4c+hpJJollY/JSTyMMbuecnyHhu4nNWG0yWGCqeIiIiBrXyBRu4F3BWrTsz/XpSwSVUcrnRbJIce9GHN/g/yiiI3GuijSd432CqxFDxFBERMcDGTxVXsILJkXsWzky/cXzBHjyIFbPPo2S2QxjOwezVceUGvNzGqzzJByZjdTkVTxEREQM2cAebuI8Ay+P+Mu7DLXjpYSyXmHM/P2AMgwjXswAhruQRo5m6moqniIhIEtkEqecRGngxcs8DlHEEo1lFnwxdzSzb14/u3MwZMXfcNODnMWYby9TVVDxFRESSaB0XsI4p+Dptj9SdKeTQ31AqSRUnMZajGInz44ibIPB9HuY9lhpO1jVUPEVERJIkRIAGngaiR2CWcT6DeY1yLjSaTVKDFw//43p6UhK5ZwP/4wtzobqQiqeIiEgSbORKviEfK2a/RhcWFfycYo4zmExSjRsXJzK648oLePkjr/MkH5mM1SVUPEVERBKsmf/SwL1AADdBXOTiZQ8q+Tt5jDUdT1LQP/g+xzGa8MbyQWxu5w2zobqAiqeIiEgCtTKd5piVyS4gh2KGsowyLjUXTFJaDh5OZkzHlbO763LqeI+vTcbabSqeIiIiCdLE/WzkJNp4Nu5+KdcYSiTp5IdM5ITIQiMX9bRzMv9gMy2mo+0yFU8REZEEaWEaEB6vgkJOoA+zKedGo7kkPeTiZSpnEXuOewPtfM1Gc6F2k4qniIhIF/Mxk42MJBCzGMQCCjiWPA4yF0zSziC6M5o+RM+18nIBT7GJZsPJdo2Kp4iISBey8VHHdwjyBS4acWHhZghFXEURV5mOJ2nGi5u3+TG96UZ45HMxG3g4TVe46xBYERGRLmLTRjtvYVMHhB+x25TzKDlMNBtO0lY3CulPGetoxPlT5WIB60zH2iUa8RQREekCIRqp5WDqOTnuvodReNnXUCrJFH9hMkXk4Txud/MUn/FHZpqOtdNUPEVERLpAO88SYGHk2qKMYu6jnFlY5BtMJplgIoO5gPFx915jmaE0u07FU0REZDf5eBYfz8fdc9GdAi7HRamhVJJp9qNv3PUSNvA5NYbS7BoVTxERkd3QzgM0cwZBXopsemNRQjH3Gc0lmecSxnM5B0Suq2niDJ40mGjnqXiKiIjsBj//BaKb3RRyLj3ZQC7HGM0lmekwBsVdr2CzmSC7SMVTRERkF4RYTRtnEuLjuPsexmGRYyiVZLojGUpvigE3kEM7Lv7ELNOxdphl27ZtOoSIiEi6aWV/QszDBoJYwEi8nEA+t2LhNh1PMthsVjGRf0auLWAt11JJkblQO0j7eIqIiOyCEJ8C0b06c7kJD+eaDSVZoQeFMVcWNtCCz1ScnaJH7SIiIjshaD9Me6gSy44d1SzBxcHGMkl2GUZ3LmEczuN25xjNH/AaAUKGk22fHrWLiIjsINtegc/eEwhiAzZeXNYVeK0f4GKU6XiSRVrwU8gdcffe4XwOZ4ChRDtGj9pFRER2UMheBQQB5xG7hZ8cbsSil9Fckn28uMjDQxsBwn8a02HEU4/aRUREtsO2AwT8pxMMTAI7Ombj4lQsS6VTks+Lm39xEu6OR+3g4Ye8SS2tpqNtk4qniIjIdtihx7Dt55xxpWAAKzQQjzUNj/Ws6WiSxc5mBJ6Yh9crqOctVhlMtH0qniIiIttg20Fse03k2gKsUA5u60wsS9smiVl9IlsouQAXrhSvdqmdTkRExCA7tJpQ2whs3/+LecTuwu2+yWgukbBpnEw5hTjLdjz8iJlU0WQ61lapeIqIiGyFHfg92Esjj9hd9nF4vJ/hck8xnEzEMZ5KAjHXG2hlBmu2+vGmqXiKiIhsgW23Y9t1kWtn3XB3LGtvc6FEtmAwJR2vLMBFfgqfnKXiKSIi0okdmANNfbHan455xN4Ny3u90VwiWzKNE+hFMc7jdjc/ZjbVtJiOtUUqniIiIp21XwP2pugqdtdPcOUvxXKNMZ1M5Fv2opx6/JHrGlp5h7UGE22diqeIiEgM2w6B3Ra5tgDL6o1ldTcXSmQ7hkQetzur27uRZzLOVql4ioiIdLAD/4PGbhBaAOF5cq49wXuJ0Vwi2/MMx9KTIsKr23/CHBpjRkFThYqniIhIWOsUoL7jIgi5f4HChViungZDiWzfIIrZgC9yvYR6ZlNjMNGWqXiKiIiE2Y3x167+WFahmSwiOyEPNz07PV7vTYGhNFun4ikiIuJ7FTYPgqAVvefaFzwnGosksjMsLJ7jSIrxRu79mgXY2AZTfZuKp4iIZDe7DRrPhNBKrFAT+C3IexQKZ2FZqTdiJLI1lRTEzet8gVUso8Fgom9T8RQRkexmNwDNkUsLG8vaE8vKN5dJZBcU48VNdNTehUVJzAhoKlDxFBGR7OWbDZvPAyqi99z7g2ecsUgiu6oX+dzLQbjwAjmAl1dSbD9PFU8REclOoc1QeyL43gLfeggWQ8EDUPo2WLmm04nskgGUEOoY9QwBVzPPbKBOVDxFRCQ7BVeBXR9z3QjeI0Cr2CWNpdpios5UPEVEJPsEVkLjo2B1i97zjAL3AHOZRLrAsfRmMn1xDkDwUkERq1Po3HYVTxERyS6hOqieCI1/gfZaoDcU/hK6vw1Wai3EENlZblwcTh+c4mmxghauZqHpWBEqniIikl18n0KwKuZ6HRT+FFw9zGUS6UIbaI+7rqHNUJJvU/EUEZHsYfvAtxpiT3hx9wNXmbFIIl3tAgZ2bCTvBjyMpNx0pAgVTxERyQ62D6qOhvXfg0AbuPpD/klQ8T+w3KbTiXSZkZRyIv1wiqeLB1jFf1lnOhag4ikiItmi9X1oe895bQPtq6HHfyBnhNFYIomwmMa464+oM5QknoqniIhkB1dJ/LVVoMVEkrEmEZ6z7AJcTEiRx+0qniIikvk23gorDoFQPuACqxh6PQaWx3QykYS4nZEMowzwAB5uZQUBQqZjYdm2ndo7jYqIiOyOtk9g+djotasQhm0Gl0Y7JXOtpJVBvB137xMOYTQlW/mM5NCIp4iIZLZgbfx1qBnwG4kikixleMiPqXkeLHqQYzCRQ8VTREQyV9sXsPFRcPeK3iu/HFwF5jKJJEEpXh5hDHnkADmUkM+GFPiBS5NbREQkMwU2wbLDIbjRuc7pBwMehqJjzOYSSZK1+CNbx9cS4GqW8A7jjWbSiKeIiGSmti+jpRPAtwZyx5jLI5JkTQS3eW2CiqeIiGSmUBu4iqPXOUPB091cHpEku4g+9CUHcGPh4QR6mo6k4ikiIhlo1SWw/BgINULOYCj/Pgx9UycUSVbpQy77UA64sXHxJ6r4mCajmVQ8RUQks7R/DbUPR699K6Dy/0HuIGORREyZG3OCkR877toEFU8REcksVi5gxd4AV56pNCJGHRyzb6cbOIDirX9wEmhVu6Q8G5smNhAkQDWfEsJHLcsJ0kIrm/DRRIB2grTgxgKC5JCLCze5FJBLMbkUUkx/3ORQzkhcuCllOG5yTf/2RKQrtSyCVdeAuz8EVzv3ev8RvH3M5hIx5B/syTgWUIMfG/iEJsZSZCyPTi6SlOGjleXMo44qlvMh6/iczXxDC5vwsRkXLtyEcOH8xOTC+ektdsaWeytv4Y8Pv/cCbnIAi2IGUcJQyhlBMXvQi+PJoRyv4dMdRGQn2SH4eAD4q5xrKwdGLYD8kWZziRj0CNVMYWnkug85VHGgsTwa8RQj/Pioo5r3eZLlfMQ3fEQTG/DRCkQLY7hcWoCXEMGOe+H3IeKLp72V96GYjw9/voUPgHqW0MQSqpkeV0yLGUY5Y+jBcRQxnBIOxop7fCciKSXYGC2dALYP/DUqnpLVcjvNqux8nWwqnpIUIUJspppXuI8vmMUS3sUmFCmYTrGMCnbc71wsO+tcA62dfB+2paLaxFLaWEoN/8EL5NKdYkbQnTMo5gCKDf7EKCJb0LoU8kZB22fOdU4/KNjXbCYRw86gJ/+ihv/RgBsXZ1FhNI8etUvC2Ni8zr+ZzxvM4ll8tEQed4NTNK0tvA/L6bj2EH5c7jxqL6CYAorIo4giyvHgJZcSLCxyycemHQ8egjRiEcLHZqCVAPUEqMMiiAt/3L/LE/Pvcm3jfWzGPAZQysFUcBn5jMFDeVd/CUVkR635M3zzc+d1bg/odRFU/sQpnyJZ7ii+YAYNgPN32GxGcqChRUYa8ZQuVc8mPuUDpjGVL5hDK41xI5ohoo+7bb498gjgJZ9+DKOCgQxgNL0YRk8G053BFFBODrt3xrKfJtrYQCtVNLKcBhbRzAoa+Rg/67FpxNWRbWsjpABtrCLIKup4ihzyKWMypZxFCafvVj4R2QVVt0dft2+EnOEqnSIdPqEl8toGFtGi4inpbT4zeZGHeYdnaKel80YmW3xsbgEllNOXoYzmKAYxllEcSz7FeOIevHctL0V4KaKYwVRwyLd+vYHPaOVranmDRhbQzBysjmkBW3q0bwM2rdQxjWamsYEyijmVUq4mh32w9L+ZSOJ5yp35nLHXIgLAsZTyJJsAizwsDjW4eFaP2mWXtdDMC/yLx/grVXwdKWaxj9OJufYApZSxF+OYyGnswb7szUQDyXdOCB8+1rGBR6njDRp5HwhFHr+Hf3/hUd3wtRvIYRglXE4+l+DSKnmRrhdshqU/hYb3ne2TQq3OY/Y9/gGWFgOKAKymnXEsZhNBcrF4jiGcSKmRLCqestO+4GMe426e5SGAby0Qil0o5MJFXwZwKKcwkcnsw6HkpPnemTYBmpnPJh6hgecJUR03P7TzinynhJaRy+kU8Xvc9DIXXiTTLLkSqu6NXu95N/S/0lwekRR0K9XcxNrI9XgK+IjhRrLoGaDssHd5g4f4G+8wPTIH0k10a6MwC+hJX07mIo7iLPZgHxNxE8bCQxETKGICNlNp5yvquJ9mXiLIisjXIvZrYlNHGw/Rzr/J5ygKeAA32tBaZLc1fxF/3brYTA6RFJbbaaJYnsGtATXiKds1g/9xF3/kA96JFM7Y0b3wau8e9GQSJ3EJN9GPobiy8ETWFt6gmado5TGsjpXz4cfwYc7m9ZBjH0Yuf8dl7W0kq0jaCzTC17+Eqjs7blgw5lXofpzRWCKpppkgh7GMBbSTi4tH6MfZhnZiUfGUrZrPHH7HTcxiRuRe5y2G3MAY9uMcLuU0LsajQXQAQtTTxnSauYUQSwnvEPqt7aNs8NjH4OEuXK69zAUWSTdtVTDnEGj9BtxeqPwu9L0Uuh1pOplIymknxDCWsQo/AL3wsJQ9KdnmTtmJkX1DUrJdq1nFeZzGcUzk3ZjSCdEN1i3gOE7jeT7iWeZxJpeqdMZwUUoB59KTxXTjA3KZjAt35CsU/xj+DQLB4QT8F2OH1piIK5J+Vt/vlE6AoB9aqlQ6RbaiikCkdALUEODrjtP7kk3FUyLaaecGrmVfRvASL+AjSIBo2QQooojvcRmLaeLvPMc+jDcVN214mUAZL9GNheRxCR48zs+YNrhivri2/U+CbUMJ+W7HthtNxRVJD668bV+LSERfPAyK2aawEg9DyTGSRcVTAHiapxjNSP7G32ikmSDOvpvgFM9CiriIK1hENX/i7xRQaDBtenKzDwU8SCFf4LV/hNvuOPk9RHgzUMCH7b8emkZg+2eajCuSuuo/hvovIafSuc6thL1u3/bniGSxXFzcSR/yyQG89CcPj6EFRprjmeVWspIfczX/41Ug1PHHMISLEC6ggBzO4lxu5jZ6Gj7fNdPY9iqCwRsIhZ7CssEKOXM+scEV7Pggz1lY+Q+CZeaECZGU07YWZoyAQL1zXToGDv0IXIk7dEIkExzOKt6lNXL9Z3pyDd2SnkMjnlnKxuZ+HmIcB/Iqr2PjwsYdswQGjuV43mEe9/Avlc4EsKwBeDxP4vEsw7IPiJROKxjzQYFp2Jsrwfc/UzFFUkv9J9HSGb62A+byiKSJFuxtXieLimcWWskqTuJMfsTPqKOpo3BGTyYfxjCe5nme51VGZtgenKnI5doDd+6HWDlvYoW6f/vhh90C9cdD7algm5kMLpIySkaBO2aqT/EocOebyyOSJn5F947H6y564eUHhk4uUvHMMi/zGqM5lFd5q6NwOn8EbCzyyONGbmIeCzmF7xhOml0sy8LlOQqrqBpyfklk589gzB6g7S/BugGNgZKqAAAgAElEQVTQOmMr/xSRDFfzJry5H7S3QcEw6H8RHPia6VQiaaESLy68gIcaLF6gxUgOFc8sESTItfyKyZxLAw0xv+LUmoOZwIe8zy3cQm6aH2mZzizLg5X3Wyj8CkIjsMILj8KP34M1sP442Ph/erwo2Wfu+dC+AUJBqFsK/aZAfl/TqUTSwjSa4jZQepwmIzlUPLNALZsZz7H8hX8Q3b7cUUQhf+Q3vMcMRjLCWEaJZ7mHYJUugLw/gN8VXfUeAghA3R9gzXHgX7vtf5BIprBt8NfF3/NvNpNFJA3FbqcEMNDQ3tsqnhnuLWYxjMP4mMUQmcfpnFQwir2Zxztcx89MRpStsXKh8EboMR/swUT2/g3vc9U6A5bvCy2fmEookjzNX0PlidHr0jFQcbS5PCJp5kcUcwwF5OKlFzn8THM8pavdy2MczxQ2UU9s4QSLq7iURcxmL/Y0mFB2iHcsVH4MBd+HANHH7iGcR+9fjYe6Fw0GFEmwDe/C6/vA2hfAlQOjboUj3geP9hMW2VGfEWAGAdqxqAEuoX67n5MIKp4ZKEiQq7iFK7mZACGip6tbFFDADJ7jTv5kOKXsFFcJ9PwX9HwE7HyngIbnf4YC8M1psO6XZjOKJMry+yDU5rwO+aDhM5VOkZ20GD+xu/V9EXOEZjKpeGag07iKu3kSOlavOVwMYwiLmcURHGIwneyWkguh/2xwDXGKZ2R9kQ01v4flF0T/ghbJFDmdNrnO6W4mh0gam0gu5TFrPE7CzDGzKp4ZpI4GxnMO/2UWRCYROyOdp3IM85hOf7QCNO3ljYWhn0LeJL61/++mf8OSUyBoZpsMkS5XuwBa6iC3D1hu6HEojPi16VQiaacfbi6kAHAm3x1o6Kx2HZmZIapYz4lcxacsw5kE6OzBYxHieqbwB67DpZ8zMs+aq2Hjnc7r8AhoCMgbD6NngUtbY0kaa1kDr4wEf8cWcN0mwPEfms0kkqa+IcBgaiLXFrCOSnpF1n8kh5pIBqilnoO5mE/5Guc/qfN43YOXf/Ab/sgNKp2Zqt9UqLzNWfEefuxuA03zYN448K03GE5kN9UujJZOgNo5ENTpXSK7or3TIzJ7C/eSQW0kzVWzkbGcxyqqY+5a5JDLa/ydH3CmsWySJJXXwZBpEPJEV73bQNOXsOAYaFf5lDRVtk/8cZjlY8Ft5vGgSLrbCy8XkAcdx2b+gAIGGNjLU8Uzja2imolczmpqIWaovJA83uNhjuJAc+EkubqfCcNeBLsg/qSjxk9hzqEQbDeZTmTnNSyDmecCBVC0Fwy5GCZNN51KJK3ZeAkvPF6EmRFPzfFMU3U0chCXs5iVRHcUD9CDEt7gHsayl8l4Ykr9h7DwCGdlu40zAmoDeWPgkPfAU2w4oMgOmj4R1s+OXh/xPAz8jrk8ImmujhDlbIq79x5lHNLpRKNE04hnGmrHx+FczWKqcOZzOkPlZRTzOnerdGaz0gNh7EwI5jnzPsNHbTZ8Ah9MhoBWu0uaaF4df92yxkwOkQxRgEXs7rcW0CNme6VkUfFMM+34OJ6b+JTVxG6ZVEQBH/AQ4xhuMp6kgrIDYPwMcBU7j9zDewRvmAlzv2cymciOCQWh/6nR69zu0H+yuTwiGSAHi4cophAPObi5ngKGa46nbM8PmMo7fE54f07wkEcOs7iP4QwyG05SR/lBMPYZCHX8cBLqeFv7PHz8E5PJRLYt0AovT4JFd4OdC3tdCZPnQ9FA08lE0t7dBGjGhQ839xNgbWSqXvKoeKaR3/Ik/2YmzkIiD84RmHm8yp8YozPXpbOKY2G/Z8BHzAlHwFd3wVcPmEolsm3LHoOaWc7rQDusmaHSKdIFfNjMijk0sw5YEHeIZnKoeKaJ5/mQW3gG5/G6M9pp4eYhrmUS4wynk5TV5xQYd1/8PRuYexlUv20kksg22Z1GYOzk/8UokolysBgVU/vygBFJ3jweVDzTwmo2cjZ/I4hF+PE6WNzKRZzDEYbTScobejmMuDm60CiI8+LtU6F+idFoIt9S0B+Khzqv3Xkw4XazeUQyyF/JpxdeSvByDfkMMVADtZ1SiqujmbHcwEo2dtxxzkW8mEk8hObqyU5473RY/Zzz2sZZdFQ+Dk6aDZ48k8lEHHN+BfN+57wu6genvgVlw8xmEskgg2njm469O93Ap+QyIsnlUyOeKczGZgr3xZROAIuDGc59/MhYLklTBz8BJWOcEc/wnM+NC2GW/ixJivj0rujrpjVQM8dcFpEM48eOlE5w/ipYriMzJdZUXuVF5sfd604xM7iFnCRv+CoZwJ0Lhz0DrrLoY/cQsPQRWPa04XAiQH7PTtcVZnKIZCAvFpNjal8f4CADNVDFM0XN5Wt+wX+IbhBvkU8Or/N/5Kp0yq4q2QMOf9p5zB4e9QzZ8Nb3ob3BZDLJdjULoMeBzp6dnkIYdx0MOM50KpGMcg5ecvHixsv55NDdwAbymuOZgnwEGM8tLGIN0Q0YQ9zGOVyHNlGWLjD3l7Dg99FRTxsoGwfnzgG3frCRJKtdDI/vFz1Za+TFcOxDZjOJZJh2bMoI0RZzbxEuRiW5fGrEMwVdxzQWsRbnP48z2nkiY1U6pevs92voPiFaOoPAhoWw4E7DwSQrrZ4Rf5zrilfMZRHJUG0db7HqDeRQ8Uwxc/mG+3gfyIGO/bV6Uc7DXGo0l2QYtxdOegnsjjPdwwX07eth02LD4STrdB+17WsR2W2lWFyOFZnjPwk4wEAOFc8Ucx7/jByt7RRPFw9wEb0oNRdKMlNBBZwwLTqbI4CzefdL50LQZzicZA1fEyx+Fsr3gdI9YdhZcMK/TacSyUjukAtCbgi52TfkwmtgjqeKZwr5Fa/wNXU4j9edeXZTmMhkxhrNJRls6GQYeFJH6ey4t/5jWHDftj5LpOu8fgXMvxNqFsGGZTB8ChRWmk4lknGW2XBPzKqev9gWVQZW+ah4pohNNHMbMwnP6QSLCkq5nTMMJ5OMd9pz4CmIvzfjRti83EweyS7V8+Kva+Zv+eNEZLdsqfCZKIEqninidB7Fh010QRHcx7n0oMhoLskC7hw48Z/R6xDga4UZvzIWSbLIwCOjry0X9D/cXBaRDDbUgp/FXN9sQe/kP2lX8UwFH7CKmazGebzuBixOZh++yxjDySRr7H0W9DgAfEQfu3/yBKx4x2wuyWxLX4CV70BRHxh6Epz+MvQ/1HQqkYzUaMMbMdemZvJrH88U0J8/sYbo5t1eQnzGNQyj5zY+S6SL1a+Gu4ZC0B9dcFQxDq78CFxu0+kk0zRWwf1DogvZcsvgyrXgzTebSyRDPWvDGaHodR7QauBbu0Y8DXuYeXGlE+AmjlLplOQr7Q+jL3ZGPMPbK61ZCJ88ZTiYZKSmtfG7J7TXQdtmc3lEMly3TtflRlKoeBr3c17H+c/gTLSooJCfMtFoJsliJ98LrkKnfIb39frvL6C9yWQqyUQ9RkGPkdHr/odDUW9zeUQy3BEWXGiDOwhlIXjcUANU8TToZt5hM35iz2P/PcdQhh41iSEuF5zwZ+ckI3BGPjevgQUa9ZQuFGiD58+HmqVQ0BcO+S2c+SpYBlY6iGSJOSF4yg/BINQFYHFw+5+TCCqehvgJcg8LcAqnM+K5N72Ywr6Gk0nWO/Ay8HZzRjwDHfdevBFaG7b1WSI77qO7Ycnzznziuiqomq+5nSIJ9p9g/IKix1U8s8stzGIj7US3T3Lxa47AgxZxSAo4477oPM8A0LgRPnnBcCjJGC0bt30tIl1ucKcHCoMMPWBQ8TTAR5CH+Rxn+yRnz8596cPZ6HxiSRGjToX8yvhz3J++FkKh7XyiyA4YdR7kdhwDbLlh/x+bzSOSBS5zwSQbSoJwEPA3r5kcKp4GTGU+62jFWVDkAtxcx0GGU4nE8ObCKX905nqGz3Fv2ghfvmU4mKS9NR/BP4+Gxnoo3xumzIKRZ5tOJZLxfuWDd/zQEIS5PliiR+3ZIUCIv/MFkIMz4uliOD04hxGGk4l0MvY0KKiIP8f96Z+aTCSZYPrV0LzBeV3zJayeazaPSJaY7o++DgKvB7b6oQml4plkL7Kc5TQRPo8d3PxSo52SivJL4Jhr4u+tWwwrF5rJI5kh0BZ/7W81k0Mky+zj3vZ1sqh4JtlfWER0+yQ3xeRwPnsbTiWyFYf9EPJKnNc2EAjB8zebTCTpbv8rwN0xuazbHrDvRWbziGSJG3JgrxD0CsIvvHCG5nhmvjU0MZsNRDeMd3E9+xlOJbINhd1g7BnO4/bwI/dFr8Om1YaDSVp65zZ49lJo9cOQ4+GKhVBUYTqVSMZrt+HUelgSgJogPNIKjYbWiqp4JtFf+ZzoaKcHC4sLGG44lch2fPc3hHdfwAZa2+D9J0wmknTka4bXfgF2x4Thxa/B+i/NZhLJElUhWBFTNNeF4CstLsp8T7A85sriFIYwkGJjeUR2SHlf6DPG2Xk4PDn91Xsg4N/WZ4lsQaeNA3VSkUhS9HXBgJjG18uCoZrjmdleZCXVxE+iv5A9DaUR2Umn/jK6n2cQ2LAa1n1lOJSkFU8eHHglkfI5/mLoN95oJJFskWvBD3OgTwj2Av5bCiU6qz2zPckKYs9k708h32WQ2VAiO2rE4eAudEY8gzgF9Pk7DIeStBEKwoOTYeZd4LNhwo/hzIdMpxLJGi+1wa8aYW0Qlvjh3iZzWVQ8k2AjbbzEOpwvt3NE5qkMMJxKZCcUlcOEM51Rz/CG8gv+ZziUpI0V78OXr0avZ90D/ratf7yIdKkFvk7XBmdKqXgmwTw20YqFs2G8G7D4CSMNpxLZSYedH7+6fUMVzHnFcChJC568+GuXF1yGJpiJZKEjcuML31G5xqKoeCbDPSwnuoWSm34UMYQiw6lEdtLYo6Ggh/OoPYDz/qNXt/NJIkDFcBh5qvPa7YWzH4ju5SkiCbenGw6ynD08z86B20vMZVHxTDAbmw+pxxnt9AIW5zIYt770ko72PT56frsNvP7v6PY4IlvSXAu3jYeFLzpzhE/4E+x/oelUIlnl3A3wfruzh+e0Jviw3VwWtZ8Ee5lqNhIgfESmhYcLGGg6lsiuOfy86ON2P9BUD3UbDIeSlLbwGVi/zHltAzP+YjSOSDb6PGZOpw18oTmemetdanFGOnMADyXksA9lhlOJ7KLRh0F+qTPiGfbcncbiSBrI6/RML7/UTA6RLHZyfvR1oQWT8rb+sYmm4plg/6GG6JfZxQlUmowjsnvyC2GvTnsvfvyOkSiSJgZPhCGHABYU94LzHjCdSCTrjPVA7xAMtmFaD9jT4BRrFc8EaiHIOkJE53e6OBadSyxpbtLZzrMacN5//D40N5hMJKmqZhn8bhwsnuVMzTj1dhhykOlUIlnl7Rb42QZYF4AVfvjtJrN5VDwTaBo1OFtnOfM7c/ByFr3NhhLZXcPGOyUi/GYD61aazSSp6cPHoKnjbznbhhl3mc0jkoW+6rSH51eGTztW8UygOTQSHe300oscCvEYTiWym4aNg259oqOeAWD6YyYTSaoq6rHtaxFJuKMLoSym7Z1ueDdHFc8Eeot6IucSY+kxu2SOwWPAh/MWAj6ZbTiQpKThR8Gg8c6enf1Gw7ka8RRJtoANe3ihmwvOKIJ7DVcRFc8EaSPEWuLHtydgcMdWka500InREc8Q8OlcCIW29RmSbVbOhz8cAN/Mc85qP/EmqBhqOpVI1jlvHcxrh9oQPNMEM1rM5lHxTJDZ1NOMTfTEIjgLPWaSDLHnaGd+p4+OfT2DsHmj4VCSUj54DHwdf8PZIXj3H2bziGSpbwLbvk42Fc8E+Zp2wBN5y8dNkeZ3SqbY7zDn/G0b580Xgpk6t11ilHbaOq5EW8mJmHBecfR1dxccV2AuC6h4Jsy7NAJuwiOe+1KKOzLfUyQDlFfGr2xf+pnhQJJSRhwPA/aF3CIYdhic9WfTiUSyTsiGUAC6BWEo8N8+0N/gHp6g4pkwn+AjOuLpZW8M/4gh0tVGToiOeAaBzxYaDiQp4+s58LuD4OsF0NwMEy/99gioiCTcv+rg7lqoDcLXPvhDCpxwrOKZIEvx4Xx5nRHPg7SwSDLN0JHREc8gsHyp4UCSMj58Avxtzmvbhln/NJtHJEut9m/72gQVzwSoJ0gQL7FzPEeQv53PEkkz++zvrGi3cd5XVUFzk+FQkhLK+277WkSS4owSKIppelPKzGUJU/FMgPm0EsAiuqLdxd7kGU4l0sUGDoWQ5axqDwKWC5oaTaeSVDDiWBg4HgrKYOQxcPYdphOJZKWv26A0CGVB+EU3+GkKbK6j4pkAmwjhLCxyAx6K8FCsL7VkmgFDnMeoYXYIFn5oLo+khm8Wws0Hw7J5UFcHY78DJT1NpxLJOvVBOGcFVPmhLgh/roF1etSemRbhI1o8LfqSg0sr2iXTuN1QHvPjsw1s0l6eWe/DaeBrjV6/96i5LCJZrC4ArTFjA34bNqh4ZqZvCOLM7XQDXirIMZxIJEEq+jnzO0M4j9y/WmI4kBjXfcC2r0UkKfrnwDExe3geXAh7p8ByE+1ongDVhIjO7wwxWMVTMlV5T6dwhq1eZSyKpIgRR8KeE2HdYhi8H3z/TtOJRLLSah/Ut0FBAPYrgul7gDcFHr5qxDMB1mIRHfH0UInh3VpFEqVygLOwKDziuXGT4UBiVPXXcOME+Px9qN3kLDAq0/6dIiZcsQLmNkNLCN5rgH+nyEwoFc8E2EzMpAosemlgWTJV9wqneAZwymdNjeFAYtT8l6GlPnr97r/NZRHJctWd5nNW+8zk6EzFMwHa4oon9MNtKIlIgpX3iJbOIFBbv51PkIzWc2D8dcUgIzFEBC6riL4udcM5KbCVEqh4drkgNrXYxK5q92pFu2SqXr2jj9mDQFu74UBi1J4Hwj7HQHEPGHE4/PgR04lEspIvBG/XQn4Ahrhg+jAYngILi0DFMyFy8RK7j6dGPCVjWS7nyMzwiOfGWsOBxJimzXDdgTD/DWdbrYIeGvEUMeSetfDUBmgNwfI2uGON6URRKp5dbCM27ViESye4Oq5FMlBlH2f/zvCIpyfXcCAx5svZsH5l9PqD58CfIpPKRLJMjX/b1yapeHaxYsATGfF0YeGhj77MkqlaWsDnBV8Z+EuhoR2am02nEhMqBoIr5ntd977g1VZyIiZcUAElHQ9bLeBHvY3GiaNG1MVqgAAunNFODzYutNxCMpcLKAfywM4HukF+ikwkkuTqNQQOORdKK2DofvB/L5lOJJK1nqmBYDt0s+G+ofC9XqYTRal4drFugLvjEXt4xLOHHrVLpvpsMfHfRtywOoUmE0ly2Db86kR483HYsB6aW6DvXqZTiWSl9+vg5uXQHIJaP9z8telE8VQ8u1gtEIxsIO+MeKbInq0iXW/sGIjdPswCBuqIxKxTWw2LZkavV30J3ywyl0cki23oNLV6k9/52TBVqHh2sUrCczyjI54lhjOJJEwgCLEj+jbOaJdkl+JyKCqPXntzoUc/c3lEsthR3WBEYfT6in5gpdCDVxXPLtaARcC2IOiCoAvbtqjWo3bJVLl58deFhVBYYCaLmJOTB5Ovgm59oe8wuOlp6NHXdCqRrPR8NVQ1QkEQrukHf0uxWS8qnl2sAsgLuXBGgSwIufCk0BC3SJeq2QgxI/y0hwwHEiMevAke+Q1UV8G6NVA51HQikay0vh0uWQT1AWgJwp0roCbFzvVQ8exiIduiLe6OxTqNeEqm8gUgZk4z3cq38wmSkWb+J/q6rQXmvGoui0gWqw9AIGawK2BDfQrt4Qkqnl3ObfGtVezNGvGUTLW+luiIpxty9Zg9K/XdY9vXIpIUexTAyTFntJ/QE/Yo3PrHm+AxHSAT5do4p7iEABdU66ssmWr9ZqIjnkBZmck0YsrJl8Ha5eBrg5MvhUNPM51IJCu9ugE+2AxeC87sDY+MBleKPXTViGcCdA+XTpz3azXtTTLVhs0Q2bfWDRUV2/kEyTgL34abz4A1S2H9KsgvMp1IJCuFbDj3Y2f7JL8NT6yFTxpNp/o2Fc8E6BM+u9oH+KFGxVMy1coaiJzU5YZePQwHkqR7/0UIBaPX7z1nLotIFvOFoDEQf6/Wt+WPNUnFMwF620RHPG1YkWITe0W6zMYGIqUTD/RNoQOBJTn6D+t0nWJ7t4hkiTw3XDGQyJke+5fCod2MRtoiFc8E6G/hzPH0OW/rA9v5BJF0tWo9kcfsuGDEEMOBJOkO/g6MnAjdesOks+HyO0wnEslKb2+AaSvB5YeTusHMA50ymmq07CUBxrlwimeHKj8EbWfFu0jGCAahoRWndNqAC8p0TldWaWmCKw91FhYBlK/QHE8RQy5cED0u85VqmLUJjknBafca8UyAUnAetfudt9YgtGhLJck0X6xwZrOHRzxdLjhsnOlUkkzLP4uWToAv50Jtjbk8Illsc6f5nJtTdJqfimcCTPCCNzzi2THf86MU/QMgsstWrou/DtmQl2smi5hRORDyYvZuLa+AkhScVCaSBa7fg8j8zn1K4IQUHO0EFc+EKLA65jD4o2/LVDwl03zwBcQeljCoLxTkbfXDJQN1r4STfgjd+8Ce4+D26ZCjHz5Ekm3GevjbErACcEw3eP9QKPaaTrVlKp4JMsLC+cmjY8RzdqvhQCJd7ctVxK1o32Og4UCSdPfeBE9OhXVrYfEi8KXg3i0iWeDij5xH6zbwRg28s8F0oq1T8UyQfbw4e3l2jHh+2badTxBJN3OX4ox4dqxo33+E4UCSdLNejr4OBmDO6+ayiGSxhk5PVTvv55lKVDwTZFIO0b08gYWtzkkCIhnBtmFjI8457R7n/fBBZjNJ8g0dte1rEUmKXwwnMr9zdCmc0sdonG3SdkoJMrjTVzYA1Aehh77ikgneWADtQZwRT8tZ0X7aIaZTSbKdeiks/hjaW+CMK+DI000nEsk6L62B33wCBODY3vDcoVCYwl1DI54JMjEfil1E5nhiw1MpeGaqyC75dAXOI/aOEU9vPhTlGw4lSbVyGVx1Mny1GFavgk/nmk4kkpV+8CE0dzxaf30dzErh+Z2g4pkwbgv6u3GGOjvmen7QbDiUSFd57RMio5244NDRYOmEhKzy8fvQ1hK9/vANc1lEslhLMP66OYXnd4KKZ0Idl090xDMEMzTiKZkgGIRPV+LM1OmY3zlhb8OhJOmGjXamWITtNdZcFpEs9osRROZ3HtQDTuxrNM52qXgm0P55OEdndox6bmyDjSn+k4jIdi1YARuaiIx24oILjzQcSpJuyAg44nTo0Qf2PxJum2Y6kUjWeWAJ3Dwf8MMplfDO0al5PnssFc8EOrMUcjvmdxKEQBBeqDOdSmQ3vb8U51uH13lz50JfnVaTdW65HF79D6xdCx/OhKqVphOJZBVfEK78EIIdo50vrYKFtWYz7QgVzwTyWDDAi7OXZ8fI5wubDIf6/+3dd3hUVeLG8e+dkkoJEHpHmlIEwbKiSFEs2MWCiq6uHetPd1fdteuu69rLqqiou/YO6qoLCipFUAREBEWpUkJNI23K/f1xZzKTgNTkninv53nyZM5NIi8q5M25554jsrfe+QbnwaLIGs+hfSFXJxalndlTY69DIZg73VgUkXQUtmOlM6oqtP3PTSQqnvXsrDyqZzwJwudbDAcS2RtbSmFO/BPtfhiivRvTUp+DYq8tC3oPNJdFJA1l+eDafale33lyBxjU0mikXaLiWc+OaIhTOkNAGEqrYEaR4VAie+qLH6GsMjKIzHheOcJkIjHlpAugfTdo3RFufxoGDjadSCStXDcDHpwHBODCzvDOMPAkweYiKp71bHhjaO3D+Ykk8pDRa+sNhxLZUy9Oqzlu2Rga55jJIuYs/xnGngq/LIEVK+D18aYTiaSVn4vg4QWx8fgfYU3Zb39+IlHxdMGghjjrPCPbKr22FqrCO/kikUQ05UdqrO8co1mutLR4PlRWxMbfzYaw/lITccv2JjaTYLITUPF0xdg2xNZ5BmBDGfyozeQl2bw9BwrLqbG+8/j+hkOJEb0OgJzc2HjAoJp7eopIverQAEZ3oXp951/6Q5vcHX5JwtDfFC7o3zByfGaI6g3lH1xqOJTI7npjDs6G8R7Ago4t4YiehkOJES3bwpAToVlL6Pc7+Ne7phOJpI3KEAx7F15dDFYA7ugPdx+0869LFCqeLmjsg1Ob17z20QZnKwSRpLC+BN79HudmTuS0opMPMBxKjHnwVpjwKqwrgK9nwsTXTCcSSRsfLINpa53XNvDAXKNxdpuKp0t+H3+ElQ0FlfDyGmNxRHbPJ4sgECY24+mFK4caDiXG/DCv5njRfDM5RNJQ7ZOJMhP8pKLaVDxdMqQZtI8eoRlw3l7QQR+SLP7+GbEjMn3QvTV0zjccSowZfHTN8eHaUkvELf2aw+FtnNfZPhiXZHMAKp4uurQ9NdZ5flYAS0oNhxLZmcXr4ceNOLOdXudt7BDw6q+PtHXwEdC5GzTNh0v/BCNHmU4kkhYWbYa+r8KXa8DvgZdGwMn7mE61e/Sdw0UXtie2n2fkJKMnfzabSWSn7pjsbAWGBXihQQ6cM8BwKDEmHIYLRsKyJbB5I4x/CFbqaUkRNzyzEDZHdjILhOFf35nNsydUPF3UOguOaxEZhIEQPLnEZCKRnVhVBBMW48x0RmY8j+8NzZJk3w6pe6UlsGFdbBwIwK9aNyTihrzMHY+TgYqny67sglM6I7OeFVXwyGLDoUR+y7ivoTy6MbgFXh/cfYzRSGJYo8Zw2JGxcbtO0Ecz4CJuGNAMejYCbNivKfxzkOlEu0/F02XHtoa+DYndcg/AwwvB1tZKkmg2lcGjX+HcYo9sGH9AR9inmeFgYly/QyGzAeS1gHuegoaNTCcSSXmPfAvHvweLN0KeB949Fjo3Np1q96l4GnBZV5w1nmHAhuXF8PSPhkOJ1PbKAiiujAwiR2Q+NtJkIkkEM6+qeY0AACAASURBVD6HB++EklJYvx5uvtJ0IpG08Ezces7CSngrSZfqqXgacGk36JyLM+sZAKrgoe+0obwkkLIA3DK15rWOeXBQOyNxJIGsX1tzvE4bEou4oU2DWuMkXWqv4mmAx4LLeuDcaredt582wTOLDAcTifr3AiiqwnmgKLJ/5yMjwbIMBxPjBh8F7TvFxudcbCyKSLoorIB+TaF9DjTLhMv3h/N6mU61Z1Q8Dfljb2iTQ2xD+SDcMxvKg4aDidg2/N8UYhvG+6FDMxiRZJvFSf3weqFRvvODc34buOga04lEUlowDENfgX/OglWF0DITHhriTGIlIxVPQywLbu6HUzwjVpXAvxYYiyTiuGtW5CeguA3j7zsSsv2Gg0lCeOyfMPcbZ436mjVw542mE4mktOWFMK8gNv5hIyzZYi7P3lLxNGhsb2hda43GXbNh7VYzeUQoD8LfZuMUzsiG8fu2gFN6GA4mCaO0ZMdjEalTLXMhLys2bpix7XrPZKLiaVj1GauR7ZWKtsK9X5tMJGnt/MlQGSK2YbwP/jEUMryGg0nCuOByaNLUeZ2VDWNvMJtHJMX97xc4sDl0yIVD2sDEUdA023SqPafiadjxnZ1NYAni3LoKw+NzYP56w8Ek/Xy/Gd5ZilM4ATwwrBOcoLWdEmfLFigPOGvTm7WGfXubTiSSsiYvhdPfhEm/wMot0D8fhnQ0nWrvqHgmgLeOBcvGKZ4BCFfBFZ9oeyVx2ajJkTXHHiADLD/ce5jhUJJw7r4FikucuzTLlsLTj5tOJJKypq10/qjFj5OdimcC2LcpnLgPzqwnQBhmrIDx80ymkrTy6CL4sQhnttMDWHBZXziwpeFgknBqb6mlLbZE6s0htbZO/l0KbKVs2bYOa0wE5QHIfwTKKqh+0r1hBiy43Nm3W6TeFFZBy9ehqgpn2h1okgkLR0HrHKPRJAF9MxtOPApKiqF7T/jkS8jPN51KJOWsKYaxE2H+OmicA8f2gFuPgCzfzr82kWnGM0Fk++G+oTilM/KgUclWuPpDw8Ek9R07GarCVD9MhBceOVSlU7bvyy9hUzFUASUVzr6vIlLnzn0T3lsEy7bAvNUwvFPyl05Q8UwoYwdAv1bEznEPw8SF8Pwcw8EkdT37M3y1Ie6CB4a3hXP1QJH8hkcfjL1esRwmvGssikgqW7xhx+NkpeKZYF4fBbk+Yue4B+Da9+CXjYaDSepZWgqX1dq7K8sD4wdp3Z78tuhWSlHNmpnJIZLiTto39jo3A45MkfkAFc8E070Z3HAYsRONwlBcCmNe1h0tqUNbg3DstMj/Z36czeIteOpQ6JDEOxNL/fvXOKd8ejxwxllw8qmmE4mknJv+Cy98BXleGNMXZl4KPZqbTlU3VDwT0G3D4NAOOLfbI/t7zlwKN75vOJikjusWwE+lVBdO/DCiPYzpYjiYJLwHHoSCzVARho8nw7p1phOJpJQpP8O9U6AiCIXlMGUJ9GllOlXdUfFMQJYF40+HBhk1rz84BSYtNpNJUsi/V8EzK3AKZ+RhokZZ8MER4NEtdtkB24b33ouNN26E6dPN5RFJQRu37nic7FQ8E1SPFvDUaXEXbAgG4eRxsDnF/icUF31bBJd9j/NHP7Jfp88H7w4Cv/46kJ2wLOjePTb2eKBrV3N5RFLQge2gcxOqd46/9nCjceqcvtMksHMGwkW/iwxCzltZBQx6ALZWmkwmSWlLAM6YD+XRxcJewA/X9oBhLUwmk2Ry5z2Q1Qp8zWHU2dCvn+lEIilj0mLofS8sK4CODeH9C+Dvx5lOVbdUPBPcI6dCn9Y4P/lEnnRfvBJGj4NQ2HA4SS7Hfwe/VOD8sY+cUDSkOfxjP8PBJKmcexVUAEEvvDEJPvvCdCKRlHHzB7C1ynm9YjMsWms2T31Q8UxwORnw5dXQIBPnQaNIAX3/W/i/1wyHk+Rx0RKYUYQzywnggda58HJ/reuUXVdR4TzxUM2CTz83Fkck1XhqtTJvCra0FPwtpZ7GOfDRWMjy4zzpHgBC8OjH8Oj/DIeTxHfbCniugBpPsGdkwMcDoU2W4XCSVLKyoHF23AUbTjneWByRVFJQDIO7OCcZAgxsH7fcLoWoeCaJw7rCg6dTfZQ2IaAKrn8R3phpMJgktifWwp2rIoPI7XW88E5f6NvQYDBJWs8/BpkZkX08T4SB/U0nEkl6G0rhwAfg/ilQHoDzD4KZ1zkbjqQaFc8kcvlQuOEYnNvtkTPdg1Vw9sMweb7hcJJ4XtkIVy2PuxDZPunerjAy31AoSXpX3wqVVRAOwxvvw/TZphOJJL1Ji2FVYWz81jzweX/785OZimeS+ccZcNbBkUHkSfdQFRx3F3y50GQySSgTi+DC5WDXWr95eSv4czsjkSQF2Dasr3V+b0GKHCAtYlDbvB2PU4mKZ5LxeOClK2BYL2JPugchUAFH/RU+08ynTCiGU5dBpU31rXU8cFo+PJoih/2KGZYFl58XG3frAsNTbJNBEZcVl8OEedCnBTTJgr5t4PXzTaeqPyqeScjrgY9ugoO74az5jGzLWBmAk+6CSXNNphOjJpbAqSsjZ7BH79N44cim8Ep38OkJdtlLBw8ETybghQaNnQMIRGSPnfM8PPQZLFgNpWXw73OgXwrfmFLxTFIZPvjwLzCw1qEhpeVw3K3wnk6xSz8vFMMpqyEcLZcewA+HNoZJPSFDf9ylDtz4DwjbgBfmLoTX3zedSCSpTV8aex0IwezlxqK4Qt+JklizhvDpXdCnY+RCZKulYCWcegf8W1stpY8HiuHC9ZFdDyJbJuGFIxrBpG5ms0lq8dea4fT7zeQQSRGDusRe+71wUCdjUVyh4pnkGuXA1w9D/y7U2GrJroLz/wZ3vGAwnNS/kA1XbYEbtkSWXETWc2LBkIbwfifI0R9zqUOP3QFZmc7rQ/rDmdrHU2RPVAXh7HEwdSG0zoET+8CHV8D+KXybHVQ8U0KmH6beC0P6UmOrJYJw+zNwzu1QobPdU0+VDScXwuNbcWY5I4UTH4zMg087QsMU3Y9DzFm0DCrCgB++XgyzvzOdSCQpPfEZvDobSithbSF4w3DUvqZT1T8VzxTRKAc+uhtOiJ5yEI69vfIRDLsC1mjXk9SxKQz9N8IHlcT+GPuctzGNYUJrHYUp9eOtj6k+BSsUggmTTScSSUoFxTsepyoVzxSSlQET74SxJ1NjqyVCMHMe9DsLvtFen8nv0wB02ww/hCIXPFSfSvSnPHghH7wqnVJPunWqOe7e2UgMkWRm29CrNeRElkxbFlw2xGgk11i2bdumQ0jdu/9V+OOjOOe625H3OPuAPvVXuOg05390SR62bcM9lXBHGVYwcnpAVCbwUCO4PMdUPEkXm7bAoNGw4lfo2glmvAYNG5hOJZI0bBvOeBze+toZD+4J/zwLDuqy469LFZrxTFE3jIaJ/4TcbGL9JAzhSrjkFjjtaucumSQHe3UY++hyuKUSKwix9ZxAAws+barSKe74z/vw4yqosOH7ZXDTw6YTiSSVhatjpRPgi8WpfVJRbSqeKeyEw2Das9C1Q+RCXNF8dzK0OQK+mmckmuyG0Pshwn0rYFIIsLCj6+vwQY9MWJIPgzIMp5S08dPyWuMVRmKIJKvczJpjrwey0mhXMhXPFNevB3z3Bhx7WK0P2LC+AAafBVfdrtnPRGRvtgleHCR8YhA2EymcAB5sPNijM+HrPGilJ9fFRScNc9bsRJ0y3FwWkSRTUQX/+QIGtgfLBp8XHhvj7MudLlQ800B2Fvz3X/DU7XHrOkNAGAKV8Ph42HcYTJ1pMKTUEJxsU7WfTfhZZ3NWZyG2hzAe7GwPPJWN9UoDaKiFuuKyYQfDgL6ABzwZkKMlHiK76pJn4LY34ZtfnO2TvrgJLk+zn91UPNPIpWfC/PegQ2tqPvVuw5IlcOQoOP8qKCwymzOd2euh4kyoOgoogOi2NTYeZ8aztxdrZi7WpZk7/geJ1JcPvoSvfwB8zpZt1z1gOpFI0pgSt7NMMASzfzaXxRQVzzTTpzvMexcuOjPuYhgIQSgA/34Feh4Mj40zlTA92eVQ9SSUdYTQG3HX8Tg/I/g9WFf48czPwtpff2zFoNrbYWh7DJFdNiDuyXXLggPScDcyfQdLQ00awzP3wAfPQptWRO/jVt9+L1gH19wABwyCzz43GDQN2GGo+gSK9oPKK4CKyHWcnwdsLOzWPnyz/Xif8GFpU3gx7fjD4ejISRUeD/z1IrN5RJJA0VYYeitMnAnNs+HI3vDqVXB4GpxUVJuKZxobOQy+nwQXnxt5ViBSQK2Qs+h5/jw45jgYcRwsScPbAfWtYgZsPgyKjgF7ee2PWoTxYN1ikbnGg6efCqckCJ8PiitxbrV74J+vwGatzxHZkb+9DVO/d/bw3FAE7RrDmYeaTmWGimeaa5IH4+6DyW/CgH6AHdkd0gYr7MzITZkC+/eHU0+HDTp2c6+VfwkFx8HGQRCMe6ArRGzy2eoP2Ssh804VTkkwhSUwcwHVx2au3QjzfjKdSiShbdm643E6UfEUAIYeBrMnw8P3QvPmxApoRCAAH7wH7dvCyGPh+wXOT26y67Z+Cr8eA2uOgMqPYtfDRG+rQ7AJZLwD2V+Dp72hoCI70rgBdGgVG2dlQlf9zyryW9ZsgmY5kBU5HjPTD1ceazaTSToyU7ZRWAi33Q3Pj4fycuda9PY7YWcm1AIGDoAbboKRJzl332T7it6H9bdBcG71ier4cf4deiNveKDBXZBzLVjanUYS3Rdz4aTroagUDugJU56EhrmmU4kknLWbYcC1znuAEw+B+y+Ebm3M5jJJM56yjbw8eOR+mPsNnDPa2eC2ev1nuPrMHL6bA+eNggO6wN9ugdJSg6ETTMVPsPYOWNAYlp8IlXOd69Gf8iK7WBHygP8yaFEMuTerdEqSeG4CFJY6/xPPWQz3/cd0IpGE9PGcWOkE+PL79C6doOIpO7BPF3hhPMycDuecW/OwkuhLH7BmFTxyN/RpDqOGwtRPoKjQRGKzbBsKxsOio+D7HlBwO4SLIx+LfE5k4wDsbMg6H9qWQtMnwdJkkSSTTUU7HosIAO3ya47b52//89KJbrXLLlv4PYx7CsY/5dx69xC7dRy90x69hZzbAEacAqedB4cdaSpx/QsUQtEMWPs0lHwOFMX+HdQ+etcD+JtA4zOg5f3gaeB6XJG68dEMOOkGCAQhKwO+GAcH7mc6lUhCeW8m3PgibCmBUBg6tYQXr4NeHU0nM0vFU3ZbQQG88gI8+DcoL46tU6xerxgRvSWflQWDhsPRp8HwU6BRnoHQdahsFWyeBmtejpTNUuf3vb31m9GZYX8HaPoHaH49eDW7KcluUxH0GQNrNwAeOO9YePGvplOJJIx1W6DThVAZcMYNsmHNi9BQy6lUPGXv/O9DeOxemPcVEKxZPKNFLFpILSA7A1q3h8Ejoe+hMPRU8NWeGkwwW1dAyRJY/TpsnAplPzu/Nx/blu7o79MXed9oMLS4BhqfAFaC/z5Fdtk7n8Npf4mNvV6omlJzPY5IGpu3FPpfXfPa8vHQsYWZPIlExVPqxJJFMOFVeP0Z2LDOuRZ/yzn+fbSQenC2lchrBgccAV37wQFHQn5baNbayG+DUADWTIbCH2D9VChdCOXLYtlh2yUG8UXTC+T2gObnQuvLwdfMwG9CpL7NWQwHXhzbU61Ta1j2ptlMIgnCtuHd6XD9s7B8A2DB4b1g6t/1sxmoeEodCwRg2RJ47kH4ajKsXRErbRlsW0ShZqkD8HugSXPIbw1d+kCrLtCmG7TuCrl50LQtZObu/hHR4RBUFkNVERSvhKJlzlvhQihdCcWLIFzilMftZY6KL5rRAprRHNr9HpqfAo1/t3u5RJLSrc/CP16GYBhGDYPXbtW57SLARQ/Ccx87r9s0gxvPhotGQHam2VyJQsVT6tWsqTB7Cnz4Emxa5cwo+omthYx/MAm2nRmNfx+9nukHgtCoCfgzISMHMrLA53HeWwHnn2tXgl0FdjkESsAbhFA5+C2qHzOvPXMZXziJe1+7HFsWNBsIzYdDp7GQ1QYs/SQr6eTQK2HmD7Hxq3+Fs4aZyyOSACqrIPuEmgesfPw3OHqguUyJRtt+S706eIjzdtUdUFkBH/4H5k+Drz+B8s0QDjgdMFrq4m9nW3HvvXHv7YDzvnwzRPa33+5ay/g3H86v44PY3kbULJPxwpFfOxT5+iCQ1x2aDIQO50KjXpDTYc//vYgkvYItOx6LpKEMPzRpAJtLYtdaNjGXJxFpxlOMKSuFaRNg6Xz4YTqsWQKFG2re3o5/OGl7b1HxM5e78j4q+s+J/joWzl8cvixoOQCa9YeOp0JuB+dNRCLufx3++LTzulEOzH/GWespkqYWrYAx98CytRCywOOHv4yG60eZTpZYVDwlYYTDsH4lFG+C76fCmh+h4BfYtBwKfwWfBeHKbW/Fw45LZrRYVq8t9QBhyGoCmU2gaUdo2sMpmg3aQJsjwK8tj0R2bNUG6HUxlGwFPHDh0fDc/5lOJWLMwEtgzk+x8Qd/h5Fa878N3WqXhOHxQKtOzlv3Adt+vKoC1v8CdhgKfgJsp5CGg1Be5BznGQ6A13LWYGblOLu8NGjhvG+6j3PefPM+zsezm7r8GxRJJZ/Ng5IKqu89vDVNxVPS2q8bao5XbzSTI9GpeErSyMiCdr2c1+37mM0ikva6ttnxWCSNvDYJerSFgk2Ax1nXOfIQ06kSk57DFRGR3TeoF/z5DPBFVkdnZEF5pelUIq67czyMvg2++NZZ0nXz2TBnHLRtbjpZYlLxFBGRPTP9JwhGdr796id4+mPTiURc9+ZnsdeBIGT5VTp3RMVTRET2TFnljsciKa6iMlIy4x7T7trWWJykoOIpIiJ75pYzwR95uCg3C4b2NZtHxEVzF0PHE+CTL6FxFvTsALddCKNHmE6W2FQ8RURkzxy2H2TnAj7YGoLRjzj3GkXSwJ8fg/WbnddFJc7BXbdfZDZTMlDxFBGRPfPzOiiuoPqMsRUbYEOx6VQirqj9M1ZVwEyOZKPiKSIie6ZnW2iVFxt3bgEtG5vLI+KCrWVwyR1QUAB+C7Chc1sYe4bpZMlBxVNERPZMXi48eSn4IgfcLi+EN2eZTiVSr254AJ55GxYthUAl3HQ+fPcqtNGT7LtExVNERPbcF4shGDnE1gae/mxnXyGS1H5YWnO8tQwa5JjJkoxUPEVEZM/VvrXeSrfaJXW9/T/I8ABhZ2xZMHKw0UhJR8VTRET23LXHwtD9YuOVG6GiylwekXry0Asw6hqYPAO8Now+BiY9DSMONZ0suah4iojInsv0Q2FZbDxjCbzwhbk8IvXk7Umx16EwtG4Gw3Ue+25T8RQRkb1TWWsfmQrtKyOpZemqbddxdu9kJErSU/EUEZG9c9fpsROMMvyQ18BsHpE69N5k6DkSPpkGOdnQszP8+SK4RNsn7REVTxER2TvH9oOsHMAHVcAlL8CqTYZDidSNv42DQGQSv6wczj4e7r3eebBIdp+Kp4iI7J3CMiippPoEo0AI1hSaTiWy1wIByPDVvKatk/aOiqeIiOydVo1heNyT7Xm5kN/QXB6ROvD5LGhxEEyfBZkewIZhh8ClZ5pOltxUPEVEZO9YFjwxBjIyAS8UVsFpT5lOJbJXxt4GhcXO68oKePJW+PR5Z52n7DkVTxER2XtLNkCVDXgBC+b/uu3T7iJJpLyy5ti2zeRINSqeIiKy9/q3h0ZZsXHHZuDVtxhJPjO+gW6HQ8EasCInFPXqBmefaDZXqtDfCiIisvfaNoFx54HHC3hgRTGMfcN0KpHddsbl8PNy5wx2OwhP3A6z34XGWrZcJ3w7/xQREZFdsLYEwt7Y+OMfzGUR2QO2Detr7QTWOl/rOuuSZjxFRKRu9GlTc9xcU0SSPKZMg/7DoWEGELnF3q0zDBtkNFbKUfEUEZG6MbwH3DQCLA/ggTlr4f4pplOJ7FTpVjj59zB/IWzeAhkWPHYnzJoIjRuZTpdaVDxFRKTu2B6wfTgruSx4c77pRCI7tWkzFJfExlUBOOxAaJJnLlOqUvEUEZG6s09+zXFOhpkcIrvotbfh5HMimzJEbrHv3wv27WY0VsqybFs7U4mISB2xbTj3ZXhlbuSCF8adBhcfbDSWyPYs/gl6HwqhkDNukge33wznn6lb7PVFM54iIlJ3LAuaNcS51R653f6abrdLYlq2IlY6AbYUwiVjVDrrk7ZTEhGRutW5ac2xR3McknhuvQtefROyfVAeACw4/mjIytrpl8pe0K12ERGpW8EQnPYSTFwUueCDx0+AsbrdLonh9bfh7N/Hxt27wZWXw0XnQWamsVhpQT+GiohI3fJ5oVM+4I+8WfDWQsOhRGKWLas59npg7MUqnW7QrXYREal7XZvVHFeFIBzWbXcxavNmOOc8+GqWU4CCNmDBGaeZTpY+9DeAiIjUvbEHwem9YuMZK+HmyebyiAC33gGTP4XSUggHYehh8MZLcOtNppOlDxVPERGpex4PdKs16znpFzNZRICqKli1sua1jh3gtJPM5ElXKp4iIlI/+reuOS6pgk1lZrJIWps0GVq1g/9+7Oz4Bc56zvPHmM2VjlQ8RUSkfozqDdcdCliAB5YUwZgJplNJGrriSiiJHIlp23DJxTBrBgw+3GyudKSHi0REpP7sE326PWLuOmNRJP2EwzBtGhQX1bx+yMHQaz8zmdKdZjxFRKT+DO4Amd7YuCQAny43FkfSh23D6LPgqOFQuBE8kXPY+/WDU7Su0xhtIC8iIvXrpQUwZiLVt9zzsmDj1c7miSL1ZOH3cED/mtdeewOOPVanE5mkW+0iIlK/2jSixrebwkqoCEJuhrFIktpmfQX/eR4sG2wAC7xeGDxYpdM0/bgpIiL163dtYP8WcRc8cOEk516oSB2b8w0cNRieGwdWyCk6fj888ig0a7bTL5d6puIpIiL1K9sPH58OeHFmPr3wxk8wWw8aSd379H8QCMTG3bvCliK4+BJzmSRGxVNEROpfo0zw+nC+7UQ2Uly0xWQiSTGbNsLok+CZx5wfcaL26+3MeEpiUPEUEZH6l+OHR4dWd07wwCWfw9cFJlNJCvnT1fDfibB+nVM8u3SGM8+GJ8aZTibxVDxFRMQdV+wP3Zvj7Ovpg0AYPlllOpWkgO/mwg/f1bx2+pnwwsvQtKmZTLJ9Kp4iIuKevvnETXvCK7/Amq3G4kjyu/kaGHYALFkYu8WekQEnnGI0lvwG7eMpIiLu2VwBI96HORtw5j68MKozvHmk6WSShDash31b1rx28TUw+vfQt5+RSLITmvEUERH3NM2C4R1wbrdH5qcWbIFg2GQqSUITXofrLgCfVfP6H65Q6UxkKp4iIuKu0V0hN25D+R9LYNRUY3Ek+UyfAleMhs/+Cx47Vj5vuBW6djebTXZMxVNERNzVrxn8eyg19vWcsArWlhkOJslgQwF8/G7N8wc6dYZlxXDjHeZyya5R8RQREfcNyI/b1xPneaO7vtNpRrJDUz6CwzvDfx5zFmtEHToEGjY0lUp2hx4uEhERM55fApd9BVVx6ztfPAzO28dcJkloJwyEBXNi414DYNgJcOWNkJlpLpfsOs14ioiIGRd0g26Nal6bsEqznrKNGVPgzCGwbEnN66ecA9ffptKZTFQ8RUTEnAu61hy/swaunrP9z5W0tKEALjoBZn0OpcVgRR4k6n8IjL7YbDbZfSqeIiJizvW94Mr9cB408gMWvL5Ss54CwMJv4Zn7oDzujAHbhgmz4e3pkNvAXDbZMyqeIiJi1ontcIpnZCprQxWMnK69PdPcV1PgjEPghQerfyQBoM8A582jBpOU9J9NRETMOqo1PNgffB6ceuGDjwrgvTWmk4khmzfAa09BIBC7tk83+L874aXJKp3JTE+1i4hIYmjxAWyojI175cGHv4OOOeYyieveeBruGQuhEIRw3gDOvhxu/5fJZFIXVDxFRCQxvLQSfv8NhGycmU8vDMmHKYNMJxOXBINwYC4EqmLXGjSD/oPgHy9Cozxz2aRuaLJaREQSw7kd4Nb9qD7NCAu+LoQP1xsOJvXNtuG+q2FEa6Cq5scefROenKDSmSpUPEVEJHGc3Q4axz1KstWGk+bC7EKjsaR+ffQKvPYYFG6sXuULwJGnwsAjTCaTuqZb7SIikliWlMLA6VAconp+5OSW8GpfyPIajSZ1a9M6ePAaWDgLfl0B0X0M2naBh9+HLvvG9u2U1KAZTxERSSzdGsDRLanxLeq9TXDa98YiSf247Rz49A1Yt8KZ5Yx2zNMvh332U+lMRZrxFBGRxFMahPMXwDsbcApopIROPQCO0GK/ZDfzvzDxaZj1PyiriF0fcQ6MvAAOHm4um9QvFU8REUlMi7dCn9kQjPs25cmA9/aFE5qZyyV75ZcFcNEBEAo64zAQBLJy4IU50KmnyXRS33SrXUREElPPXHilFzSNPmzkd1rKbSthTeVOvlgSTUUZPHU93Hc+2MHY9YwMOPdP8PQ0lc50oBlPERFJbH9eBvetrnmtXQbM7Q/5fjOZZLfd93v434uxcQCwgUOOg/s+NBRKXKfiKSIiiW1rCC79GV7eUPP6mJbwVBfI0ZPuiWz6BHj7IVgyB8pLY9f3OQAGHA3n3gQ5Dc3lE3epeIqISHLY5xtYGvckCj4YngeTdX82Ua1YBJf0ja3njPJ44J+fwf7aozPtaI2niIgkh3f3hR7ZkYEX8MCnxfCn1VAZ3tFXiss2roY/D4EbDsV5cijKghHnwZ0TVDrTlWY8RUQkefxcAfstgED8t64MuCofHm1nLJbEbFoNj1wM33wUuxbEeS5s4NFw78emkkkiUPEUEZHkMnELXLECVgeo3nbcZ8FZ+fBcG8jQruMm2DY8cj5M/Y8zjpZNgO6/g4NOgFOvcbZNkvSl4ikiIsnnwyI4BYAnRQAACH9JREFUflncBQ/gg6uawv0tVT5dtuAz+HoiTHwkds3GeXI9uwHcPwM69zGVThKJiqeIiCSnNwvh6tWwLoxTPC3n/UG5MKUt5OgxBjd88i94dqzz2ia2pNPjgwsehv5HQbvuptJJolHxFBGR5PVyEYxZ7TQewLn17oEjc+GJfOiufT7ry9fvwIcPwKqFUFwUux4CbAsueABOvM5YPElQKp4iIpLcppfBiaths01ssxYvNPfDvFbQxmcyXcrZ/Cv8+AU8eR6EQs616G11gNNugaEXQstOhgJKQlPxFBGR5PffrXDGOthq49xyjzx05PHBQ43h6lzDAVPDvPfhidMhWFmzbAJkt4SDToGLnnD26RTZHhVPERFJDaUh6LUOVobiLvrBsuCqbLitITRVI9oTc9+FKY/BinlQvCV2Pfrker/j4I869lJ2gYqniIikjl8CcNUW+KiS6k3mo/b1wodNoLNuve+qX2bCks/hvb+AHdkbKUzsAaKeR0Lvo+GosZCR/Vv/FJEYFU8REUk9ZxfCqxXb+YAf/pEDf1JL2pFAJXw5Dl6/evsfrwJ6DIbrPoJM7cspu0HFU0REUo9tw4QKOLsIyuM/4MfGgj5+GJeNdYhmP+PZNrx8Icx+EbCgKhy3YUBEn5Fw3nho2NxZxSCyO1Q8RUQkdX1ZBX8sgVlBog8cOd/0vNhZFpyfiXWrH6tNeq/9LF4L71wJaxfAuiWxE4fiHyDqdQx0PRyOvE631WXPqXiKiEjq+0MpjK+KlE4P0bWfYbzQwgM3+PFc6sVqlF5TeEsmwbyX4eepsHFF7HoAp3RaXghnwEFnw7nj9LS67D0VTxERSQv2ohAcXwZLnW97NhY2XiBSQDt7sP7oxXuOldIFNFgFsx6DdfPh21cgHNkEIBR5g9jT6qc+BEOvNZNTUpOKp4iIpA17UxjGBQjfFYByi+itd6eAWk4B7WhhXWPhO8vC09pw4Dr060xY/DYsmwq/znGuxd9Kj77OyIUTH4b2B0K7/Y1ElRSm4ikiImnH/jZM+J4A9sQwdtA55z2+gAaxoImF7zzwnQXeQwwH3kOVRTDzbtj0Iyz5CEKRfZC2N7vZan/oNhL6nQFtVDilnqh4iohI2gpPChO6MYQ938YOOaXTBkI4s6EAeME+BLzHQdaNYCXwOkc7DKEq+PZx+PF1KF4FZQWRj7Ht7CYWtBsMLfvCiLshq5GR2JJGVDxFRCTthT+zCY61Cf8M4aDlbLkUJzorGPID3SH7OvB1h8zDTaStqXQ1/PoplK2HWXdBVXHNGc140d9HbhtotA8MuAT2P9fdvJLeVDxFREQiwt9B5R8g/ANQFrseLWzhuNdBwNMdPN2gwbng6wQZPcGbV3/5AiWw4SvAAwsfh9JVsPlHqCx1Ph4ithVS9Mn0KMsHTftCfi848iHIaVZ/OUV+i4qniIhILfZmqLwKQtMhuCJW5qKlM/62dVVkHARoBFZ78HeCnEPAkw0Nj4FgIWT3Bizw5jgl0LadDdij7wPFzvZFVRth61Lw5cHaCYAH1n0OZSuhohDKNm6bJRh5X52DWPFs3g+a7AP7XwqdjqrXf20iO6XiKSIisgOVz0HgE6j4BELFzrXtFdDq2/HUvM1tWxC0wcqGynLw5EAoA4Ll4GsFW1dBRisoXQN4IOSJPATkobrxxv8zo4Vye7929Ndvsi90PB6a9oTeF+iEIUkcKp4iIiK7ILQOqiZDxTQofQ4IQzC8bQHd3vrK6Ixk9GPxnx/YwceifqtkBgGP3ymyjbrCwNvB3wDaHAa+rDr5bYvUKRVPERGR3WQHIVwGhfdB4BconQqBdWB7IBCu9bnEbn/vqHiG4z6vqtavF/2Y5YPKIGTmQ4fRzm37bmOgSa+6/h2K1A8VTxERkb1kByC4AYIboegtCFc4ZTRU5KzZDGwGsqGq3Pn86lv10aIaeW8D5ELVVsjtBr5cyGwJzY+GcBW0PQVy2oE3W7fPJTmpeIqIiNQjOwjBTWA1gK1zwN8CKtdGHjjqAaU/QIPeECwByw+N+kKgCDKamE4uUvdUPEVERETEFQl8/oKIiIiIpBIVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQVKp4iIiIi4goVTxERERFxhYqniIiIiLhCxVNEREREXKHiKSIiIiKuUPEUEREREVeoeIqIiIiIK1Q8RURERMQV/w/gkDXaTQhY5AAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x320625390>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x3203c98d0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "A = [-2 0; -1 3];\n",
    "t = linspace(0,2*pi,300);\n",
    "x1,x2 = (cos(t),sin(t));\n",
    "x = [x1.'; x2.'];  Ax = A*x;\n",
    "subplot(121)\n",
    "scatter(x1,x2,40,t,\".\",edgecolor=\"none\",cmap=\"hsv\")\n",
    "axis(\"equal\"), axis(\"off\")\n",
    "title(L\"x\")\n",
    "subplot(122)\n",
    "scatter(Ax[1,:].',Ax[2,:].',40,t,\".\",edgecolor=\"none\",cmap=\"hsv\")\n",
    "axis(\"equal\"), axis(\"off\")\n",
    "title(L\"Ax\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **left singular vectors** and their associated **singular values** are found from the major and minor axes of the ellipse. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: sigma not defined\nwhile loading In[2], in expression starting on line 4",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: sigma not defined\nwhile loading In[2], in expression starting on line 4",
      ""
     ]
    }
   ],
   "source": [
    "normAx = sqrt( sum(Ax.^2,1) );  # 2-norms of all the vectors\n",
    "σ1,k1 = findmax(normAx);  \n",
    "σ2,k2 = findmin(normAx);\n",
    "σ = [ σ1 σ2 ];  println(\"σ = \",sigma)\n",
    "U = [ Ax[:,k1]/σ1 Ax[:,k2]/σ2 ];  println(\"U = \",U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **right singular vectors** are the pre-images of those vectors on the ellipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2x2 Array{Float64,2}:\n",
       " -0.469368  -0.880524\n",
       "  0.883002  -0.474001"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = [ x[:,k1] x[:,k2] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: U not defined\nwhile loading In[4], in expression starting on line 13",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: U not defined\nwhile loading In[4], in expression starting on line 13",
      ""
     ]
    }
   ],
   "source": [
    "subplot(121)\n",
    "scatter(x1,x2,40,t,\".\",edgecolor=\"none\",cmap=\"hsv\")\n",
    "axis(\"equal\"), axis(\"off\"), title(L\"x\")\n",
    "hold(\"on\")\n",
    "plot([0 0; V[1,:]],[0 0; V[2,:]], \"k\" )\n",
    "text(1.2*V[1,1],1.2*V[2,1],L\"v_1\")\n",
    "text(1.2*V[1,2],1.2*V[2,2],L\"v_2\")\n",
    "\n",
    "subplot(122)\n",
    "scatter(Ax[1,:].',Ax[2,:].',40,t,\".\",edgecolor=\"none\",cmap=\"hsv\")\n",
    "axis(\"equal\"), axis(\"off\"), title(L\"Ax\")\n",
    "hold(\"on\")\n",
    "plot([0 0; U[1,:].*σ],[0 0; U[2,:].*σ], \"k\" )\n",
    "text(1.1*U[1,1].*σ1,1.1*U[2,1].*σ1,L\"\\sigma_1 u_1\") \n",
    "text(1.1*U[1,2].*σ2,1.1*U[2,2].*σ2,L\"\\sigma_2 u_2\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having searched over only a finite number of vectors, we have approximately (up to choices of signs in the singular vectors) found the SVD of this matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(\n",
       "2x2 Array{Float64,2}:\n",
       " 0.289784   0.957092\n",
       " 0.957092  -0.289784,\n",
       "\n",
       "[3.25661653798294,1.8424029756098448],\n",
       "2x2 Array{Float64,2}:\n",
       " -0.471858  -0.881675\n",
       "  0.881675  -0.471858)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U,σ,V = svd(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full versus thin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of a rectangular matrix, some of the SVD is essentially dead weight. We consider (as always in this text) the case with $m>n$, a tall and skinny matrix. The range of $A$ can have at most only $n$ dimensions. The last $m-n$ columns of $U$ therefore describe the remainder of $\\mathbb{C}^m$ in which the range is embedded (i.e., the orthogonal complement). That description can be changed without affecting the matrix at all.\n",
    "\n",
    "Algebraically, we have \n",
    "\n",
    "$$U = \\bigl[ U_1 \\: U_2 \\bigr], \\qquad S = \\begin{bmatrix} S_1 \\\\ 0 \\end{bmatrix},$$\n",
    "\n",
    "where the breaks occur after $n$ columns in $U$ and $n$ rows in $S$. Therefore $US=U_1S_1$, and we change nothing by replacing the SVD with $U_1S_1V^*$. The book calls this the **reduced SVD**, though it is also called the **thin SVD**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "full U = \n",
      "[-0.39581592558291195 0.46223308060818946 0.35763867012477596 -0.4088944319107176 -0.5652314102975107 -0.12281575730860446\n",
      " -0.11161826541211242 -0.13971308537270824 -0.32897240789305343 0.571281283804118 -0.7016074643093034 0.2029370154694903\n",
      " -0.45102437287686287 -0.005646089309049103 -0.5009609551653109 0.03554490968208342 0.12647121755098079 -0.7268595797492602\n",
      " -0.27399964890033257 -0.5053052800683752 0.6868463382960395 0.34605499477896445 0.023397215306391606 -0.27844458122088755\n",
      " -0.5129908730857501 0.49236165992048264 0.06914683936350981 0.437193262597549 0.4111861588071027 0.359759737216207\n",
      " -0.5377849825772198 -0.5186863449311001 -0.19071098566062492 -0.440780659116373 0.05141948162997606 0.45656276184400846]\n",
      "\n",
      "full S = \n",
      "[1.8695611702425259 0.0 0.0\n",
      " 0.0 0.8021585655333503 0.0\n",
      " 0.0 0.0 0.581838619284657]\n",
      "\n",
      "\n",
      "\n",
      "thin U = \n",
      "[-0.39581592558291195 0.46223308060818946 0.35763867012477596\n",
      " -0.11161826541211242 -0.13971308537270824 -0.32897240789305343\n",
      " -0.45102437287686287 -0.005646089309049103 -0.5009609551653109\n",
      " -0.27399964890033257 -0.5053052800683752 0.6868463382960395\n",
      " -0.5129908730857501 0.49236165992048264 0.06914683936350981\n",
      " -0.5377849825772198 -0.5186863449311001 -0.19071098566062492]\n",
      "\n",
      "thin S = \n",
      "[1.8695611702425259 0.0 0.0\n",
      " 0.0 0.8021585655333503 0.0\n",
      " 0.0 0.0 0.581838619284657]\n"
     ]
    }
   ],
   "source": [
    "A = rand(6,3); \n",
    "U,σ,V = svd(A,thin=false);               # full SVD\n",
    "println(\"full U = \\n\",U,\"\\n\\nfull S = \\n\",diagm(σ))\n",
    "U1,σ1,V = svd(A)   # thin SVD\n",
    "println(\"\\n\\n\\nthin U = \\n\",U1,\"\\n\\nthin S = \\n\",diagm(σ1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The price we pay is that while the columns of $U_1$ are orthonormal, we can't call it a unitary matrix because it isn't square. Take note:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.9999999999999999 4.412484725792806e-16 -7.607664998086165e-17\n",
      " 4.412484725792806e-16 1.0000000000000007 1.856185990309045e-16\n",
      " -7.607664998086165e-17 1.856185990309045e-16 1.0]\n",
      "\n",
      "\n",
      "[0.49823508612221257 -0.13805297728366447 -0.0032501894337621643 0.12052741939355656 0.45536538775772584 -0.095095749733363\n",
      " -0.13805297728366447 0.14020122855292397 0.21593372239789466 -0.12477236849451877 -0.034277617444250195 0.19523254867540252\n",
      " -0.0032501894337621643 0.21593372239789466 0.4544167418535938 -0.217649679131227 0.1939516022265859 0.341021441473415\n",
      " 0.12052741939355656 -0.12477236849451877 -0.217649679131227 0.8021671260931621 -0.06074037392490781 0.27845870302925235\n",
      " 0.45536538775772584 -0.034277617444250195 0.1939516022265859 -0.06074037392490781 0.5103609254228962 0.0073104560859966125\n",
      " -0.095095749733363 0.19523254867540252 0.341021441473415 0.27845870302925235 0.0073104560859966125 0.5946188919552119]\n"
     ]
    }
   ],
   "source": [
    "println(U1'*U1)           # n by n identity\n",
    "println(\"\\n\\n\",U1*U1')    # NOT the m by m identity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll discuss the significance of $U_1U_1^*$ in due course."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.3",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}