Trefethen & Bau lecture notes in Julia

TB-Lecture-01-julia.ipynb

TB-Lecture-02.ipynb

TB-Lecture-03.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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/URIParser.ji for module URIParser.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/BinDeps.ji for module BinDeps.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/SHA.ji for module SHA.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/PyPlot.ji for module PyPlot.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/Conda.ji for module Conda.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/PyCall.ji for module PyCall.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/MacroTools.ji for module MacroTools.\n",
      "INFO: Recompiling stale cache file /Users/driscoll/.julia/lib/v0.4/LaTeXStrings.ji for module LaTeXStrings.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAIhCAYAAADn6rozAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XeYVOXh9vHvmbK9UpfeFBGQoigoFuy9xa7RoCZqNMZYYvslbzSmYjSxR40mamzEHkUjiqKAojTBQpHedinL9jI7M+f94+y0pcPOPFPuz3XttXMOu3qz4nLvc55i2bZtIyIiIiISZy7TAUREREQkM6h4ioiIiEhCqHiKiIiISEKoeIqIiIhIQqh4ioiIiEhCqHiKiIiISEKoeIqIiIhIQqh4ioiIiEhCqHiKiIiISEKoeIqIiIhIQqh4ioiICDfeeCN33XWX6RiS5lQ8RUREMtysWbN46KGHWLNmjekokuZUPEVERDKYbdvcfvvt2LZNeXm56TiS5lQ8RUREMtgTTzzBeeedh8fjUfGUuLNs27ZNhxAREZHE27RpExdddBGTJ0+mV69e2Latx+0SVx7TAUTaw5dffsnzzz+P2+1m5cqVPPnkkzz++ONUVVWxdu1afvvb39KvXz/TMUVEksptt93GPffcA0DXrl2ZP3/+Vh+j76/SnlQ8JeUtXbqUZ599loceegiAyy+/nDFjxvDss8/i9/s58sgjOfDAA7nxxhsNJxURSR7Tp08nGAwyZswYwCmegUCAjRs30rlzZ0DfX6X9qXhKyrv//vuZMGFC+Lquro4OHTowevRo1qxZw80338z48ePNBRQRSTKBQIDbb7+dV199NXyvrKwMgPLy8nDx1PdXaW9aXCQp77bbbiM/Pz98/dlnn3HccccB0LNnTyZMmEBpaampeCIiSefBBx/k7LPPpkuXLuF7Xbt2BYhZYKTvr9LeNOIpKa93797h14sWLWLdunUcffTRBhOJiCSvdevWcd999zF48GAmT56MZVnYts3KlSuB2OKp76/S3lQ8Ja18+OGHZGdnM3bs2PC95cuXa+K7iEirW265heeff56jjjoq5v5rr73Gueeeu90tlfT9VdqDHrVLSmtqauK2227jm2++AeCDDz5g+PDhZGdnA87GyH/5y19MRhQRSRpTpkyhU6dOW5VOcB6dQ2TEU99fJR404ikpbdKkSdx7770cdNBBeDwevv/+e0pKSsK//oc//EET30VEcB6VX3HFFeEi2VZovmeoeL7zzjv6/irtThvIS0rbvHkzt956K506dcKyLH7zm99w7bXXkpOTQ3Z2NmeeeabmI4lIRluxYgUXX3wxs2fPxu/3M2TIEN5///3wKnaAs846i3nz5rF69WpycnIYPXo0V155ZXiEVN9fpb2oeIqIiIhIQmiOp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEiqeIiIiIJISKp4iIiIgkhIqniIiIiCSEx3QAEREREYB6qsmhgO+ZSyO15FLMQr6gBwNZzrd48NKZnpSzkiGMZhPldKMP3ehDgABFlJj+LchOqHiKiIhIQtWwiTm8gxsvX/A2m1hNDTUsZz5ucmmkkSAQwCKIjR+woz7fBvxY2NgEsQgAQWz6MxQ3XkZwGH3Ylz7swzGcauT3KNtm2bZt7/zDRERERHZfA9XM4x1aaOZjnmQ9i2nERyO1tOCUyEDrmw20tH6eHwi2vvnb/DN39PF21McHADfZ5JLL2VxGb/pzKmfTk95x+t3Kzqh4ioiISLsJ0EILzbzC7azha1aziGrKo37deQPwtbm3vSLZQqxtfXx0iW37eaF/FkA2OfSmL0dwLD/lJrpQRh55e/vbll2k4ikiIiJ7pYq1VLKS2bzCFB4kiIsWWrY5Whl9L3aU0nl0Di5svLTQTBkDqWELHelJHXW4cOMhiyo24CGfNSzDSx511IfLaDDqfejf44uNEB4VtYEiivk5t3IwYxjHMe39pZE2VDxFRERktzVRy0LeZQPf8x730ExTuOiFRhq3VTwBLLJxk0MJZQzjZMDFIZxBgCCd6UVX+hGghSxydpihBR9esljMAlrwsYxFLOEbVrCML/mUGhqoYkt4NDQkesQ0VEBt4AiO5kzO5lIuo5jiPf/iyHapeIqIiMguW8ZUyvmaT3mIDSwKl7rokhkacXReW3jIo4jOnMANZJPHUI6nC/3inrWFFj7kHSoo5yWeZRlLqGRTzOhopHi6AQsb6E4PLmM84/kR/ROQM5OoeIqIiMgONVHN97zHcqYxg4djSmaoeLYdPdyP4+jCQEZyBkM50UzwbZjMJGYwlVeYyApWRGV2NvoJ4sJ53A+FFPFTruZSLmZ/BhlMnT5UPEVERGSbqljFJhYxievZvJ3RzdC2RzY2fRjDsdxKEWUM4FBjuXdFE018y9e8wkQe4QGaCbbON3VGPkMFFCCXPH7NHVzID+ijFfF7RcVTREREYgTw8R1v8hqXEIhaUx49VzKAM1czn05cwD/pxjAK6YKFZSLyXrGxmcpUbuUWvmUJzfgIhh+9W9i4AIsiiniOxzieo8kl13TslKTiKSIiIoBTwN7mSubzL/y4CIY3Pgr9Othk4aWAoZzDWfwdV5qdvr2YJfyCm5jBl9RSHy6g0fahP58yiTK6mgmZwlQ8RUREMlwDG/mImylnLuV8DbDVSnA32RTQlXN4nj4cbiRnIq1mDVfzcz5lJnU0tPlVD3nkcx5n8CT34cVrJGMqUvEUERHJUHWsYyXvMZtHqWB2zObu0XM3+3Mcl/A2HrINpjVjOSs5n8uZy3wCrWchOW8AFkdwKL/iF5zAUQZTpg4VTxERkQxUyxpeYhQNVIQ3WI9emQ5wDBMYxnjy6WwmZBKZxVxO4lw2UwPhhUeRAvoMf+MSzsIdvifbouIpIiKSQTYwiw+4hGpW4aMJiD1S0k0243iA7oymjBHGciajeup5jH9xK/e0LjhyFh2BB7A4hOF8wHMUUmA2aBJT8RQREckANXzPaiYxmwnUszZmSyQbyKMnA/gB+3Mx3RltMGnym8ZM7uCPTGMWTukMLbByczDDmcAvGaev4TapeIqIiKS5apbwFgfjozrm3HI/4CafYvbhJF6gI4NNRUw59TRwC/fwb96ijkai53568DKDFziYA4xmTEYqniIiImmqmoXM5BoqWUg9FUDsY/VC+nIOX5JLJ2MZU937fMrZXEMDfpzH7s6jdwuLWxjPBG4ynDC5qHiKiIikmQA+alnCVM6lhoUxq9UBShnJUG6iNyeTQ0dTMdPGRjZzDj/nU2YDXiL7frr4PddzKz/C03okZ6ZT8RQREUkjPqqZzNFUMjfmvh/w0olShnAkz1BIHzMB01QVNdzO/TzOa613InM/T2Esb/O3lDzVqb2peIqIiKQBG5vV/IcVTGQlr2716zl04VTmkUc3A+kyxx08wJ94GtpsKn81P+BebqCQfDPBkoSKp4iISBr4gqtYxpMEiMzhBMinNwO5lr5cSIFGORPiHT7lTG4mEPNfAk7nCN7ir4ZSJQcVTxERkRRWwfus5RWW8DQ2AWyc+Zw24KWQo3iDbhxjOGXmeY53+BkTqKG+9Y4bFx5+yEk8yi3kk2s0nykqniIiIilqMzP4hCOxCcRsk2QDh/Ey3TiRLIpNxct4X7GYUVyKHxuiFhf9ggv4KzeYC2aQa+cfIiIiIsmkmQ18xRXM5Qrs1vXqHpxTh/LozUjuow/nq3QaNpyB/Je/0pcerXecrZYe5i2u5wFsMm/sTyOeIiIiKWYGR1HJJzGnDwH05CIO4QVTsWQ7JjGD07mVYPiYTcfL/IbzOdpcMAM04ikiIpIi1vEssziWLcwAnL/E3UAxw+jPdYzkcaP5ZNtO4TA+4iFKKGy94wK8/Ij7eJ4PTUZLOI14ioiIpIBKpjKLcQC0QPghrYtcjuQrCtjXVDTZRRN4kdt4nOhN5rPwUsHLlFBgNFuiaMRTREQkibVQxWKuYzHXhe95gCwK6c1VHMoUlc4UcSsX8Q9uJXKykRsfMI47WMUGg8kSRyOeIiIiSWwBZ7GZN7eaz9mdyxnK06ZiyR4KEOAUfsX7zCV6pfvZHMpr/J+5YAmiEU8REZEkVMdclnEdVUwBnL+wPUAHjmQg9zJY8zlTkhs37/A7zuGIqLse3mU+9/OWsVyJohFPERGRJNPECuZxAEHq8BM5icjCy4F8TiEHmown7eAbVnIYv6QGH9HjgB9yN8dwgLlgcaYRTxERkSRSzVTWcA9B6gBn1bobN2VcyTDeU+lME0Pow1c8FHWCkbPH5z28xmZqTUaLK414ioiIJIkK/sFSfrLV/QIOYRgzDSSSePsZT/II7xK90n0cg/mI/2c0V7xoxFNERCQJNLOccv4ecy+b/nTmhwzidUOpJN4e4sdcw8lEVrrDdBaxiHXmQsWRiqeIiIhhK7mOBfSngdkx97tzI/vyHFl0N5RM4s3C4hecSg7e8L0WAoziTr5htcFk8aHiKSIiYlAj37KRRwFn1boz068zXbiSMq4xmk0SYz+6M5n/o4S88L06mniB6QZTxYeKp4iIiAE2LazlWpZzeviehTPTbyTfsg//wIra51HS2+EM4jD2a71yA14m8C4v8pnJWO1OxVNERMSAjfyFzTyGn2Uxfxl35268dDKWS8x5nB8znL6E6pmfINfxjNFM7U3FU0REJIFsAlTzDDW8Gb7nAUo4mmGsonuarmaWnetJR+7i3Kg7bmpo4TlmGMvU3lQ8RUREEmg9l7Ke8fjabI/UkfFk0ctQKkkWpzKCYxmC8+OImwDwI57mUxYbTtY+VDxFREQSJIifGl4GIkdglnAJ/XiPUi4zmk2SgxcP/+NWOlMUvmcD/+Nbc6HakYqniIhIAmziOlaQixW1X6MLiy7cQiEnGkwmycaNi1MY1nrlBbz8ifd5kS9NxmoXKp4iIiJxVs9/qeFRwI+bAC6y8bIPZfydHEaYjidJ6Al+xIkMI7SxfACbe5lsNlQ7UPEUERGJo0YmUR+1MtkFZFHIAJZQwlXmgklSy8LDaQxvvXJ2d11GFZ+y1GSsvabiKSIiEid1PM4mTqWJV2PuF3OToUSSSn7CWE4OLzRyUU0zp/EEW2gwHW2PqXiKiIjESQMTgdB4FeRzMt2ZQSl3GM0lqSEbLw9wPtHnuNfQzFI2mQu1l1Q8RURE2pmPqWxiCP6oxSAWkMcJ5HCouWCScvrSkWF0J3KulZdLeYnN1BtOtmdUPEVERNqRjY8qziLAt7ioxYWFm/4UcD0FXG86nqQYL24+4md0owOhkc+FbOTpFF3hrkNgRURE2olNE818iE0VEHrEblPKs2Qx1mw4SVkdyKcXJaynFudPlYs5rDcda49oxFNERKQdBKmlksOo5rSY+x6G4uVAQ6kkXdzP6RSQg/O43c1LfM2fmGo61m5T8RQREWkHzbyKn7nha4sSCnmMUqZhkWswmaSDsfTjUkbF3HuPJYbS7DkVTxERkb3k41V8vB5zz0VH8rgGF8WGUkm6OYgeMdeL2Mg3VBhKs2dUPEVERPZCM09Sz7kEeCu86Y1FEYU8ZjSXpJ8rGcU1HBK+LqeOc3nRYKLdp+IpIiKyF1r4LxDZ7Cafi+jMRrI53mguSU9H0jfmejlbzATZQyqeIiIieyDIapo4jyDzYu57GIlFlqFUku6OYQDdKATcQBbNuPgz00zH2mWWbdu26RAiIiKpppGDCTILGwhgAUPwcjK5/BELt+l4ksZmsIqx/DN8bQHruJkyCsyF2kXax1NERGQPBJkPRPbqzOZOPFxkNpRkhE7kR11Z2EADPlNxdosetYuIiOyGgP00zcEyLDt6VLMIF4cZyySZZSAduZKROI/bnWM0f8x7+AkaTrZzetQuIiKyi2x7OT57XyCADdh4cVnX4rV+jIuhpuNJBmmghXz+EnPvYy7hKHobSrRr9KhdRERkFwXtVUAAcB6xW7SQxR1YdDWaSzKPFxc5eGjCT+hPYyqMeOpRu4iIyE7Yth9/yzkE/OPAjozZuDgTy1LplMTz4uZfnIq79VE7ePgJH1BJo+loO6TiKSIishN28Dls+zVnXCngxwr2wWNNxGO9ajqaZLALGIwn6uH1cqr5kFUGE+2ciqeIiMgO2HYA214TvrYAK5iF2zoPy9K2SWJW9/AWSi7AhSvJq11ypxMRETHIDq4m2DQY2/f/oh6xu3C77zSaSyRkIqdRSj7Osh0PP2Uqa6kzHWu7VDxFRES2w/b/HuzF4UfsLvtEPN6vcbnHG04m4hhFGf6o6400MoU12/1401Q8RUREtsG2m7HtqvC1s264I5a1v7lQItvQj6LWVxbgIjeJT85S8RQREWnD9s+Euh5YzS9HPWLvgOW91WgukW2ZyMl0pRDncbubnzGDchpMx9omFU8REZG2mm8Ce3NkFbvr57hyF2O5hptOJrKV/SilmpbwdQWNfMw6g4m2T8VTREQkim0HwW4KX1uAZXXDsjqaCyWyE/3Dj9ud1e0dyDEZZ7tUPEVERFrZ/v9BbQcIzoHQPDnXvuC90mgukZ15hRPoTAGh1e0/Zya1UaOgyULFU0REJKRxPFDdehGA7Pshfy6Wq7PBUCI715dCNuILXy+imhlUGEy0bSqeIiIiIXZt7LWrF5aVbyaLyG7IwU3nNo/Xu5FnKM32qXiKiIj43oUtfSFgRe65DgTPKcYiiewOC4vXOIZCvOF7v2EONrbBVFtT8RQRkcxmN0HteRBciRWsgxYLcp6F/GlYVvKNGIlsTxl5MfM632AVS6gxmGhrKp4iIpLZ7BqgPnxpYWNZ+2JZueYyieyBQry4iYzau7AoihoBTQYqniIikrl8M2DLxUCXyD33weAZaSySyJ7qSi6PciguvEAW4OWdJNvPU8VTREQyU3ALVJ4Cvg/BtwEChZD3JBR/BFa26XQie6Q3RQRbRz2DwA3MMhuoDRVPERHJTIFVYFdHXdeC92jQKnZJYcm2mKgtFU8REck8/pVQ+yxYHSL3PEPB3dtcJpF2cALdOJ0eOAcgeOlCAauT6Nx2FU8REckswSooHwu190NzJdAN8n8FHT8CK7kWYojsLjcujqI7TvG0WE4DNzDXdKwwFU8REcksvvkQWBt1vR7yfwGuTuYyibSjjTTHXFfQZCjJ1lQ8RUQkc9g+8K2G6BNe3D3BVWIskkh7u5Q+rRvJuwEPQyg1HSlMxVNERDKD7YO1x8GGH4K/CVy9IPdU6PI/sNym04m0myEUcwo9cYqniydZxX9ZbzoWoOIpIiKZonE6NH3qvLaB5tXQ6T+QNdhoLJF4WEhtzPWXVBlKEkvFU0REMoOrKPbaytNiIklb4wjNWXYBLkYnyeN2FU8REUl/m/4Iyw+HYC7gAqsQuj4Hlsd0MpG4uJchDKQE8AAe/shy/ARNx8KybTu5dxoVERHZG01fwbIRkWtXPgzcAi6Ndkr6Wkkjffko5t5XHM4wirbzGYmhEU8REUlvgcrY62A90GIkikiilOAhN6rmebDoRJbBRA4VTxERSV9N38KmZ8HdNXKv9Bpw5ZnLJJIAxXh5huHkkAVkUUQuG5PgBy5NbhERkfTk3wxLjoLAJuc6qyf0fhoKjjebSyRB1tES3jq+Ej83sIiPGWU0k0Y8RUQkPTV9FymdAL41kD3cXB6RBKsjsMNrE1Q8RUQkPQWbwFUYuc4aAJ6O5vKIJNjldKcHWYAbCw8n09l0JBVPERFJQ6uuhGXHQ7AWsvpB6Y9gwAc6oUgySneyOYBSwI2Niz+zlnnUGc2k4ikiIumleSlUPh259i2Hsv8H2X2NRRIx5YuoE4xasGOuTVDxFBGR9GJlA1b0DXDlmEojYtRhUft2uoFDKNz+ByeAVrVL0rOxqWMjAfyUM58gPipZRoAGGtmMjzr8NBOgATcWECCLbFy4ySaPbArJJp9CeuEmi1KG4MJNMYNwk236tyci7alhAay6Cdy9ILDaudftT+DtbjaXiCFPsC8jmUMFLdjAV9QxggJjeXRykSQNH40sYxZVrGUZn7Oeb9jCChrYjI8tuHDhJogL5ycmF85Pb9EzttzbeQt9fOi9F3CTBVgU0pciBlDKYArZh66cRBaleA2f7iAiu8kOwrze0LLWubayYOgcyB1iNpeIQc9QzngWh6+7k8VaxhjLoxFPMaIFH1WUM50XWcaXrOBL6tiIj0YgUhhD5dICvAQJtN4LvQ8SWzzt7bwPRn186PMtfABUs4g6FlHOpJhiWshAShlOJ06kgEEUcRhWzOM7EUkqgdpI6QSwfdBSoeIpGS27zazKtteJpuIpCREkyBbKeYfH+JZpLOITbILhgukUy4hA6/22xbKttjXQ2s33IdsqqnUsponFVPAfvEA2HSlkMB05l0IOodDgT4wisg2NiyFnKDR97Vxn9YS8A81mEjHsXDrzLyr4HzW4cXE+XYzm0aN2iRsbm/f5N7OZzDRexUdD+HE3OEXT2sb7kKzWaw+hx+XOo/Y8CsmjgBwKKKAUD16yKcLCIptcbJrx4CFALRZBfGwBGvFTjZ8qLAK4aIn5d3mi/l2uHbyPzphDb4o5jC5cTS7D8VDa3l9CEdlVa+6DFbc4r7M7QdfLoeznTvkUyXDH8i1TqAGcv8NmMIQxhhYZacRT2lU1m5nPZ0zkAb5lJo3UxoxoBok87rbZeuQRwEsuPRlIF/rQm2F0ZSCd6UdH+pFHKVns3RnLLdTRxEYaWUsty6hhAfUsp5Z5tLABm1pcrdm2N0IK0MQqAqyiipfIIpcSTqeY8yninL3KJyJ7YO29kdfNmyBrkEqnSKuvaAi/toEFNKh4SmqbzVTe5Gk+5hWaaWi7kck2H5tbQBGl9GAAwziWvoxgKCeQSyGemAfv7ctLAV4KKKQfXTh8q1+v4WsaWUolk6llDvXMxGqdFrCtR/s2YNNIFROpZyIbKaGQMynmBrI4AEv/m4nEn6fUmc8ZfS0iAJxAMS+yGbDIweIIg4tn9ahd9lgD9bzBv3iOv7KWpeFiFv04nahrD1BMCfsxkrGczT4cyP6MNZB89wTx4WM9G3mWKiZTy3QgGH78Hvr9hUZ1Q9duIIuBFHENuVyJS6vkRdpfoB4W/wJqpjvbJwUbncfs+zwBlhYDigCsppmRLGQzAbKxeI3+nEKxkSwqnrLbvmUez/Ewr/IUwFYLhKIXCrlw0YPeHMEZjOV0DuAIslJ870wbP/XMZjPPUMPrBCmPmR/adkW+U0JLyOYcCvg9brqaCy+SbhZdB2sfjVzv+zD0us5cHpEk9EfKuZN14etR5PElg4xk0TNA2WWfMJmn+BsfMyk8B9JNZGujEAvoTA9O43KO5Xz24QATcePGwkMBoylgNDYP0Mz3VPE49bxFgOXhr0X018Smiiaeopl/k8ux5PEkbrShtcheq/829rpxoZkcIkksu81EsRyDWwNqxFN2agr/4yH+xGd8HC6c0aN7odXenejMOE7lSu6kJwNwZeCJrA1Mpp6XaOQ5rNaV86HH8CHO5vWQZR9JNn/HZe1vJKtIyvPXwtJfwdoHW29YMPxd6Hii0VgiyaaeAEeyhDk0k42LZ+jJBYZ2YlHxlO2azUx+x51MY0r4XtsthtzAcA7iQq7ibK7Ao0F0AIJU08Qk6rmbIIsJ7RC61fZRNnjs4/HwEC7XfuYCi6SaprUw83BoXAFuL5T9AHpcBR2OMZ1MJOk0E2QgS1hFCwBd8bCYfSna4U7Z8ZF5Q1KyU6tZxcWczYmM5ZOo0gmRDdYt4ETO5nW+5FVmcR5XqXRGcVFMHhfRmYV04DOyOR0X7vBXKPYx/GT8gUH4W67ADq4xEVck9ax+3CmdAIEWaFir0imyHWvxh0snQAV+lrae3pdoKp4S1kwzt3EzBzKYt3gDHwH8RMomQAEF/JCrWUgdf+c1DmCUqbgpw8toSniLDswlhyvx4HF+xrTBFfXFte1/EmgaQNB3L7ZdayquSGpw5ez4WkTCeuChb9Q2hWV4GECWkSwqngLAy7zEMIbwN/5GLfUEcPbdBKd45lPA5VzLAsr5M38nj3yDaVOTmwPI4x/k8y1e+6e47daT34OENgMFfNgtt0LdYOyWqSbjiiSv6nlQ/R1klTnX2WWw3707/hyRDJaNiwfpTi5ZgJde5OAxtMBIczwz3EpW8jNu4H+8CwRb/xgGcRHEBeSRxflcxF1MoLPh813TjW2vIhC4jWDwJSwbrKAz5xMbXIHWD/Kcj5X7D7DMnDAhknSa1sGUweCvdq6Lh8MRX4IrfodOiKSDo1jFJzSGr++jMzfRIeE5NOKZoWxsHucpRjKGd3kfGxc27qglMHACJ/Exs3iEf6l0xoFl9cbjeRGPZwmWfUi4dFqBqA/yT8TeUga+/5mKKZJcqr+KlM7Qte03l0ckRTRg7/A6UVQ8M9BKVnEq5/FTbqSKutbCGTmZfCADeZnXeZ13GZJme3AmI5drH9zZn2NlfYAV7Lj1ww+7AapPgsozwTYzGVwkaRQNBXfUVJ/CoeDONZdHJEX8mo6tj9dddMXLjw2dXKTimWHe5j2GcQTv8mFr4XT+CNhY5JDDHdzJLOZyBmcZTppZLMvC5TkWq6Acsn5FeOfPQNQeoM1vwfre0DhlO/8UkTRX8QHshOxbAAAgAElEQVR8cBA0N0HeQOh1OYx5z3QqkZRQhhcXXsBDBRZv0GAkh4pnhggQ4GZ+zelcRA01Ub/i1JrDGM3nTOdu7iY7xY+0TGWW5cHKuQfyv4fgYKzQwqPQ4/dABWw4ETb9nx4vSub54hJo3gjBAFQthp7jIbeH6VQiKWEidTEbKD1PnZEcKp4ZoJItjOIE7ucJItuXOwrI50/8lk+ZwhAGG8sosSx3f6ziOZDzB2hxRVa9BwH8UPUHWHMitKzb8T9IJF3YNrRUxd5r2WImi0gKit5OCaCPob23VTzT3IdMYyBHMo+FEJ7H6ZxUMJT9mcXH/JIbTUaU7bGyIf8O6DQb7H6E9/4N7XPVOAWWHQgNX5lKKJI49Uuh7JTIdfFw6HKcuTwiKeanFHI8eWTjpStZ3Kg5ntLeHuU5TmI8m6kmunCCxfVcxQJmsB/7Gkwou8Q7AsrmQd6PwE/ksXsQ59H796Og6k2DAUXibOMn8P4BsO4NcGXB0D/C0dPBo/2ERXbV1/iZgp9mLCqAK6ne6efEg4pnGgoQ4Hru5jruwk+QyOnqFnnkMYXXeJA/G04pu8VVBJ3/BZ2fATvXKaCh+Z9BP6w4G9b/ymxGkXhZ9hgEm5zXQR/UfK3SKbKbFtJC9G5930YdoZlIKp5p6Gyu52FehNbVaw4XA+nPQqZxNIcbTCd7pegy6DUDXP2d4hleX2RDxe9h2aWRv6BF0kVWm02uszqaySGSwsaSTWnUGo9TMXPMrIpnGqmihlFcyH+ZBuFJxM5I55kczywm0QutAE15OSNgwHzIGcdW+/9u/jcsOgMCZrbJEGl3lXOgoQqyu4Plhk5HwODfmE4lknJ64uYy8gBn8t0YQ2e168jMNLGWDZzC9cxnCc4kQGcPHosgtzKeP/BLXPo5I/2suQE2Pei8Do2ABoGcUTBsGri0NZaksIY18M4QaGndAq7DaDjpc7OZRFLUCvz0oyJ8bQHrKaNreP1HYqiJpIFKqjmMK5jPUpz/pM7jdQ9enuC3/InbVDrTVc8HoGyCs+I99NjdBupmwayR4NtgMJzIXqqcGymdAJUzIaDTu0T2RHObR2T2Nu4lgtpIiitnEyO4mFWUR921yCKb9/g7P+Y8Y9kkQcp+Cf0nQtATWfVuA3XfwZzjoVnlU1JUyQGxx2GWjgC3mceDIqluP7xcSg60Hpv5Y/LobWAvTxXPFLaKcsZyDauphKih8nxy+JSnOZYx5sJJYnU8Dwa+CXZe7ElHtfNh5hEQaDaZTmT31SyBqRcBeVCwH/S/AsZNMp1KJKXZeAktPF6AmRFPzfFMUVXUcijXsJCVRHYU99OJIibzCCPYz2Q8MaX6c5h7tLOy3cYZAbWBnOFw+KfgKTQcUGQXTRoLG2ZEro9+HfqcZS6PSIqrIkgpm2PufUoJh7c50SjeNOKZgprxcRQ3sJC1OPM5naHyEgp5n4dVOjNZ8RgYMRUCOc68z9BRmzVfwWeng1+r3SVF1K+OvW5YYyaHSJrIwyJ691sL6BS1vVKiqHimmGZ8nMSdzGc10VsmFZDHZzzFSAaZjCfJoOQQGDUFXIXOI/fQHsEbp8IXPzSZTGTXBAPQ68zIdXZH6HW6uTwiaSALi6coJB8PWbi5lTwGaY6n7MyPeYCP+YbQ/pzgIYcspvEYg+hrNpwkj9JDYcQrEGz94STY+rbudZj3c5PJRHbM3whvj4MFD4OdDftdB6fPhoI+ppOJpLyH8VOPCx9uHsfPuvBUvcRR8Uwh9/Ai/2YqzkIiD84RmDm8y58ZrjPXpa0uJ8BBr4CPqBOOgO8fgu+fNJVKZMeWPAcV05zX/mZYM0WlU6Qd+LCZFnVoZhUwJ+YQzcRQ8UwRr/M5d/MKzuN1Z7TTws1T3Mw4RhpOJ0mr+xkw8rHYezbwxdVQ/pGRSCI7ZLcZgbET/xejSDrKwmJoVO3LAQYnePN4UPFMCavZxAX8jQAWocfrYPFHLudCjjacTpLegGtg8F2RhUYBnBcfnQnVi4xGE9lKXi8oHOC8dufA6HvN5hFJI38ll654KcLLTeTS30AN1HZKSa6KekZwGyvZ1HrHORfxCsbxFJqrJ7vh03Ng9WvOaxtn0VHpSDh1BnhyTCYTccz8Ncz6nfO6oCec+SGUDDSbSSSN9KOJFa17d7qB+WQzOMHlUyOeSczGZjyPRZVOAIvDGMRj/NRYLklRh70ARcOdEc/QnM9Nc2Ga/ixJkpj/UOR13RqomGkui0iaacEOl05w/ipYpiMzJdoDvMubzI6515FCpnA3WQne8FXSgDsbjnwFXCWRx+5BYPEzsORlw+FEgNzOba67mMkhkoa8WJweVfu6A4caqIEqnknqC5ZyO/8hskG8RS5ZvM//ka3SKXuqaB846mXnMXto1DNow4c/guYak8kk01XMgU5jnD07Pfkw8pfQ+0TTqUTSyoV4ycaLGy+XkEVHAxvIa45nEvLhZxR3s4A1RDZgDDKBC/kl2kRZ2sEXv4I5v4+MetpAyUi4aCa49YONJFjlQnj+oMjJWkOugBOeMptJJM00Y1NCkKaoewtwMTTB5VMjnknol0xkAetw/vM4o52nMEKlU9rPQb+BjqMjpTMAbJwLcx40HEwy0uopsce5Ln/HXBaRNNXU+hat2kAOFc8k8wUreIzpQBa07q/VlVKe5iqjuSTNuL1w6ltgt57pHiqgH90KmxcaDicZp+PQHV+LyF4rxuIarPAc/3HAIQZyqHgmmYv5Z/hobad4uniSy+lKsblQkp7yusDJEyOzOfw4m3e/dREEfIbDScbw1cHCV6H0ACjeFwaeDyf/23QqkbTkDrog6IagmwODLrwG5niqeCaRX/MOS6nCebzuzLMbz1hOZ4TRXJLGBpwOfU5tLZ2t9zbMgzmP7eizRNrP+9fC7AehYgFsXAKDxkN+melUImlniQ2PRK3qud+2WGtglY+KZ5LYTD0TmEpoTidYdKGYeznXcDJJe2e/Bp682HtT7oAty8zkkcxSPiv2umL2tj9ORPbKtgqfiRKo4pkkzuFZfNhEFhTBY1xEJwqM5pIM4M6CU/4ZuQ4CvkaY8mtjkSSD9Dkm8tpyQa+jzGURSWMDLLgx6vouC7ol/km7imcy+IxVTGU1zuN1N2BxGgfwA4YbTiYZY//zodMh4CPy2P2rF2D5x2ZzSXpb/Aas/BgKusOAU+Gct6HXEaZTiaSlWhsmR12bmsmvfTyTQC/+zBoim3d7CfI1NzGQzjv4LJF2Vr0aHhoAgZbIgqMuI+G6L8HlNp1O0k3tWni8f2QhW3YJXLcOvLlmc4mkqVdtODcYuc4BGg18a9eIp2FPMyumdALcybEqnZJ4xb1g2BXOiGdoe6U1c+GrlwwHk7RUty5294TmKmjaYi6PSJrr0Oa61EgKFU/jbuF9nP8MzkSLLuTzC8YazSQZ7LRHwZXvlM/Qvl7/vR2a60ymknTUaSh0GhK57nUUFHQzl0ckzR1twWU2uANQEoTnDTVAFU+D7uJjttBC9Hnsv+d4StCjJjHE5YKT73NOMgJn5HPLGpijUU9pR/4meP0SqFgMeT3g8HvgvHfBMrDSQSRDzAzCSy0QCECVHxYGdv458aDiaUgLAR5hDk7hdEY896cr4znQcDLJeGOuBm8HZ8TT33rvzTugsWZHnyWy6758GBa97swnrloLa2drbqdInP0nELug6HkVz8xyN9PYRDOR7ZNc/Iaj8aBFHJIEzn0sMs/TD9Rugq/eMBxK0kbDph1fi0i769fmgUJfQw8YVDwN8BHgab7B2T7J2bPzQLpzATqfWJLE0DMhtyz2HPeXb4ZgcCefKLILhl4M2a3HAFtuOPhnZvOIZICrXTDOhqIAHAr8zWsmh4qnAQ8wm/U04iwocgFufsmhhlOJRPFmwxl/cuZ6hs5xr9sE331oOJikvDVfwj+Pg9pqKN0fxk+DIReYTiWS9n7tg49boCYAX/hgkR61ZwY/Qf7Ot0AWzoini0F04kIGG04m0saIsyGvS+w57i//wmQiSQeTboD6jc7riu9g9Rdm84hkiEktkdcB4H3/dj80rlQ8E+xNlrGMOkLnsYObX2m0U5JRbhEcf1PsvfULYeVcM3kkPfibYq9bGs3kEMkwB7h3fJ0oKp4Jdj8LiGyf5KaQLC5hf8OpRLbjyJ9ATpHz2gb8QXj9LpOJJNUdfC24WyeXddgHDrzcbB6RDHFbFuwXhK4BuN0L52qOZ/pbQx0z2Ehkw3gXt3KQ4VQiO5DfAUac6zxuDz1yX/A+bF5tOJikpI8nwKtXQWML9D8Jrp0LBV1MpxJJe802nFkNi/xQEYBnGqHW0FpRFc8E+ivfEBnt9GBhcSmDDKcS2Ykf/JbQ7gvYQGMTTH/BZCJJRb56eO92sFsnDC98DzZ8ZzaTSIZYG4TlUUVzfRC+1+Ki9PcCy6KuLM6gP30oNJZHZJeU9oDuw52dh0OT0999BPwtO/oskW1os3GgTioSSYgeLugd1fi6WjBAczzT25uspJzYSfSXsa+hNCK76cxfRfbzDAAbV8P67w2HkpTiyYEx1xEun6OugJ6jjEYSyRTZFvwkC7oHYT/gv8VQpLPa09uLLCf6TPZe5PMD+poNJbKrBh8F7nxnxDOAU0Bf/4vhUJIyggH4x+kw9SHw2TD6Z3DeU6ZTiWSMt5rg17WwLgCLWuDROnNZVDwTYBNNvMV6nC+3c0TmmfQ2nEpkNxSUwujznFHP0Ibyc/5nOJSkjOXT4bt3I9fTHoGWpu1/vIi0qzm+NtcGZ0qpeCbALDbTiIWzYbwbsPg5QwynEtlNR14Su7p941qY+Y7hUJISPDmx1y4vuAxNMBPJQEdnxxa+Y7ONRVHxTIRHWEZkCyU3PSmgPwWGU4nsphHHQV4n51G7H+f9l+/u5JNEgC6DYMiZzmu3Fy54MrKXp4jE3b5uONRy9vC8IAvuLTKXRcUzzmxsPqcaZ7TTC1hcRD/c+tJLKjrwpMj57Tbw/r8j2+OIbEt9JUwYBXPfdOYIn/xnOPgy06lEMspFG2F6s7OH58Q6+LzZXBa1nzh7m3I24Sd0RKaFh0vpYzqWyJ456uLI4/YWoK4aqjYaDiVJbe4rsGGJ89oGptxvNI5IJvomak6nDXyrOZ7p6xMqcUY6swAPRWRxACWGU4nsoWFHQm6xM+IZ8tqDxuJICshp80wvt9hMDpEMdlpu5HW+BeNytv+x8abiGWf/oYLIl9nFyZSZjCOyd3LzYb82ey/O+9hIFEkR/cZC/8MBCwq7wsVPmk4kknFGeKBbEPrZMLET7GtwirWKZxw1EGA9QSLzO12cgM4llhQ37gLnWQ047+dNh/oak4kkWVUsgd+NhIXTnKkZZ94L/Q81nUoko3zUADduhPV+WN4C92w2m0fFM44mUoGzdZYzvzMLL+fTzWwokb01cJRTIkJvNrB+pdlMkpw+fw7qWv+Ws22Y8pDZPCIZ6Ps2e3h+b/i0YxXPOJpJLZHRTi9dySIfj+FUIntp4Ejo0D0y6ukHJj1nMpEkq4JOO74Wkbg7Lh9KotreOYZ3c1TxjKMPqSZ8LjGWHrNL+ug3HHw4b0HgqxmGA0lSGnQs9B3l7NnZcxhcpBFPkUTz27CPFzq44NwCeNRwFVHxjJMmgqwjdnx7NAZ3bBVpT4eeEhnxDALzv4BgcEefIZlm5Wz4wyGwYpZzVvspd0KXAaZTiWSci9fDrGaoDMIrdTClwWweFc84mUE19dhETiyC89FjJkkT+w5z5nf6aN3XMwBbNhkOJUnls+fA1/o3nB2ET54wm0ckQ63w7/g60VQ842QpzYAn/JaLmwLN75R0cdCRzvnbNs6bLwhTdW67RClus3VckbaSEzHh4sLI644uODHPXBZQ8YybT6gF3IRGPA+kGHd4vqdIGigti13Zvvhrw4EkqQw+CXofCNkFMPBIOP8+04lEMk7QhqAfOgRgAPDf7tDL4B6eoOIZN1/hIzLi6WV/DP+IIdLehoyOjHgGgK/nGg4kSWPpTPjdobB0DtTXw9irth4BFZG4+1cVPFwJlQFY6oM/JMEJxyqecbIYH86X1xnxPFQLiyTdDBgSGfEMAMsWGw4kSePzF6ClyXlt2zDtn2bziGSo1S07vjZBxTMOqgkQwEv0HM/B5O7ks0RSzAEHOyvabZz3a9dCfZ3hUJIUSnvs+FpEEuLcIiiIanrjS8xlCVHxjIPZNOLHIrKi3cX+5BhOJdLO+gyAoOWsag8Algvqak2nkmQw+AToMwrySmDI8XDBX0wnEslIS5ugOAAlAbi9A/wiCTbXUfGMg80EcRYWuQEPBXgo1Jda0k3v/s5j1BA7CHM/N5dHksOKuXDXYbBkFlRVwYizoKiz6VQiGac6ABcuh7UtUBWA+ypgvR61p6cF+IgUT4seZOHSinZJN243lEb9+GwDm7WXZ8b7fCL4GiPXnz5rLotIBqvyQ2PU2ECLDRtVPNPTCgI4czvdgJcuZBlOJBInXXo68zuDOI/cv19kOJAY17H3jq9FJCF6ZcHxUXt4HpYP+yfBchPtaB4H5QSJzO8M0k/FU9JVaWencIasXmUsiiSJwcfAvmNh/ULodxD86EHTiUQy0mofVDdBnh8OKoBJ+4A3CR6+asQzDtZhERnx9FCG4d1aReKlrLezsCg04rlps+FAYlT5UrhjNHwzHSo3OwuMSrR/p4gJ1y6HL+qhIQif1sC/k2QmlIpnHGwhalIFFl01sCzpqmMXp3j6ccpnRYXhQGLU7LehoTpy/cm/zWURyXDlbeZzlvvM5GhLxTMOmmKKJ/TEbSiJSJyVdoqUzgBQWb2TT5C01rlP7HWXvkZiiAhc3SXyutgNFybBVkqg4tnuAthUYhO9qt2rFe2Srrp2izxmDwBNzYYDiVH7joEDjofCTjD4KPjZM6YTiWQkXxA+qoRcP/R3waSBMCgJFhaBimdcZOMleh9PjXhK2rJczpGZoRHPTZWGA4kxdVvgl2Ng9mRnW628ThrxFDHkkXXw0kZoDMKyJvjLGtOJIlQ829kmbJqxCJVOcLVei6Shsu7O/p2hEU9PtuFAYsx3M2DDysj1Z69BS5JMKhPJMBUtO742ScWznRUCnvCIpwsLD931ZZZ01dAAPi/4SqClGGqaob7edCoxoUsfcEV9r+vYA7zaSk7EhEu7QFHrw1YL+Gk3o3FiqBG1swrAjwtntNODjQstt5D05QJKgRywc4EOkJskE4kksbr2h8MvguIuMOAg+L+3TCcSyVivVECgGTrY8NgA+GFX04kiVDzbWQfA3fqIPTTi2UmP2iVdfb2Q2G8jblidRJOJJDFsG359CnzwPGzcAPUN0GM/06lEMtL0KrhrGdQHobIF7lpqOlEsFc92VgkEwhvIOyOeSbJnq0j7GzEcorcPs4A+OiIx41SWw4KpketV38GKBebyiGSwjW2mVm9ucX42TBYqnu2sjNAcz8iIZ5HhTCJx4w9A9Ii+jTPaJZmlsBQKSiPX3mzo1NNcHpEMdmwHGJwfub62J1hJ9OBVxbOd1WDhty0IuCDgwrYtyvWoXdJVdk7sdX4+5OeZySLmZOXA6ddDhx7QYyDc+TJ06mE6lUhGer0c1tZCXgBu6gl/S7JZLyqe7awLkBN04YwCWRB04UmiIW6RdlWxCaJG+GkOGg4kRvzjTnjmt1C+FtavgbIBphOJZKQNzXDlAqj2Q0MAHlwOFUl2roeKZzsL2hZNMXcs1mvEU9KVzw9Rc5rpULqTT5C0NPU/kddNDTDzXXNZRDJYtR/8UYNdfhuqk2gPT1DxbHdui61WsddrxFPS1YZKIiOebsjWY/aM1GOfHV+LSELskwenRZ3RfnJn2Cd/+x9vgsd0gHSUbeOc4hIEXFCur7Kkqw1biIx4AiUlJtOIKaddDeuWga8JTrsKjjjbdCKRjPTuRvhsC3gtOK8bPDMMXEn20FUjnnHQMVQ6cd6v07Q3SVcbt0B431o3dOmyk0+QtDP3I7jrXFizGDasgtwC04lEMlLQhovmOdsntdjwwjr4qtZ0qq2peMZB99DZ1T6gBSpUPCVdrayA8EldbujayXAgSbjpb0IwELn+9DVzWUQymC8Itf7Ye5W+bX+sSSqecdDNJjLiacPyJJvYK9JuNtUQLp14oEcSHQgsidFrYJvrJNu7RSRD5Ljh2j6Ez/Q4uBiO6GA00japeMZBLwtnjqfPedvg38kniKSqVRsIP2bHBYP7Gw4kCXfYWTBkLHToBuMugGv+YjqRSEb6aCNMXAmuFji1A0wd45TRZKNlL3Ew0oVTPFutbYGA7ax4F0kbgQDUNOKUThtwQYnO6cooDXVw3RHOwiKA0uWa4yliyGVzIsdlvlMO0zbD8Uk47V4jnnFQDM6j9hbnrTEADdpSSdLNt8ud2eyhEU+XC44caTqVJNKyryOlE+C7L6CywlwekQy2pc18zi1JOs1PxTMORnvBGxrxbJ3v+WWS/gEQ2WMr18deB23IyTaTRcwo6wM5UXu3lnaBoiScVCaSAW7dh/D8zgOK4OQkHO0EFc+4yLNa5zC0RN6WqHhKuvnsW4g+LKFvD8jL2e6HSxrqWAan/gQ6dod9R8K9kyBLP3yIJNqUDfC3RWD54fgOMP0IKPSaTrVtKp5xMtjC+cmjdcRzRqPhQCLt7btVxKxo36eP4UCScI/eCS8+AOvXwcIF4EvCvVtEMsAVXzqP1m1gcgV8vNF0ou1T8YyTA7w4e3m2jnh+17STTxBJNV8sxhnxbF3RfvBgw4Ek4aa9HXkd8MPM981lEclgNW2eqrbdzzOZqHjGybgsInt5AnMbnZMERNKCbcOmWpxz2j3O+0F9zWaSxBswdMfXIpIQtw8iPL9zWDGc0d1onB3Sdkpx0q/NV9YPVAegk77ikg4mz4HmAM6Ip+WsaD/7cNOpJNHOvAoWzoPmBjj3WjjmHNOJRDLOW2vgt18BfjihG7x2BOQncdfQiGecjM2FQhfhOZ7Y8FISnpkqskfmL8d5xN464unNhYJcw6EkoVYugetPg+8XwupVMP8L04lEMtKPP4f61kfr76+HaUk8vxNUPOPGbUEvN85QZ+tcz8/qDYcSaS/vfUV4tBMXHDEMLJ2QkFHmTYemhsj155PNZRHJYA2B2Ov6JJ7fCSqecXViLpERzyBM0YinpINAAOavxJmp0zq/c/T+hkNJwg0c5kyxCNlvhLksIhns9sGE53ce2glO6WE0zk6peMbRwTk4R2e2jnpuaoJNSf6TiMhOzVkOG+sIj3bigsuOMRxKEq7/YDj6HOjUHQ4+BiZMNJ1IJOM8uQjumg20wBll8PFxyXk+ezQVzzg6rxiyW+d3EgB/AN6oMp1KZC9NX4zzrcPrvLmzoYdOq8k4d18D7/4H1q2Dz6fC2pWmE4lkFF8ArvscAq2jnW+tgrmVZjPtChXPOPJY0NuLs5dn68jnG5sNhxLZW6/NwllY1DrH8+hhkK8TizLOF/+/vfsOk6o63Dj+vVO2UhaW3kGaUgTBElGkKBbsYkFFo7Fj/WkSNbFrYoy9REVFTewd1KgBBZUiKAIigqJUKUvdxrYp9/fHndmZXZC6e8+U9/M888ycu7vwrgL7zrnnnjs19joUgrnTjUURSUdhO1Y6o6pC2//cRKLiWc/OyqN6xpMgfL7FcCCRvbGlFObEX9HuhyHauzEt9Tko9tqyoPdAc1lE0lCWD67dl+r1nSd3gEEtjUbaJSqe9eyIhjilMwSEobQKZhQZDiWyp774EcoqI4PIjOeVI0wmElNOugDad4PWHeH2p2HgYNOJRNLKdTPgwXlAAC7sDO8MA08SbC6i4lnPhjeG1j6cdySRi4xeW284lMieenFazXHLxtA4x0wWMWf5zzD2VPhlCaxYAa+PN51IJK38XAQPL4iNx/8Ia8p++/MTiYqnCwY1xFnnGdlW6bW1UBXeyReJJKIpP1JjfecYzXKlpcXzobIiNv5uNoT1j5qIW7Y3sZkEk52AiqcrxrYhts4zABvK4EdtJi/J5u05UFhOjfWdx/c3HEqM6HUA5OTGxgMG1dzTU0TqVYcGMLoL1es7/9If2uTu8EsShv6lcEH/hpHbZ4ao3lD+waWGQ4nsrjfm4GwY7wEs6NgSjuhpOJQY0bItDDkR8ltCv9/Bv941nUgkbVSGYNi78OpisAJwR3+4+6Cdf12iUPF0QWMfnNq85rGPNjhbIYgkhfUl8O73OCdzIncrOvkAw6HEmAdvhQmvwroC+HomTHzNdCKRtPHBMpi21nltAw/MNRpnt6l4uuT38bewsqGgEl5eYyyOyO75ZBEEwsRmPL1w5VDDocSYH+bVHC+abyaHSBqqfWeizAS/U1FtKp4uGZIP7aO30Aw4jxd0ow9JFn//jNgtMn3QvTV0bmY4lBgz+Oia48O1pZaIW/o1h8PbOK+zfTAuyeYAVDxddGl7aqzz/KwAlpQaDiWyM4vXw48bcWY7vc5j7BDw6p+PtHXwEdC5GzRtBpf+CUaOMp1IJC0s2gx9X4Uv14DfAy+NgJP3MZ1q9+gnh4subE9sP8/InYye/NlsJpGdumOysxUYFuCFBjlwzgDDocSYcBguGAnLlsDmjTD+IVipqyVF3PDMQtgc2cksEIZ/fWc2z55Q8XRR6yw4rkVkEAZC8OQSk4lEdmJVEUxYjDPTGZnxPL435CfJvh1S90pLYMO62DgQgF+1bkjEDXmZOx4nAxVPl13ZBad0RmY9K6rgkcWGQ4n8lnFfQ3l0Y3ALvD64+xijkcSwRo3hsCNj43adoI9mwEXcMCAfejYCbNivKfxzkOlEu0/F02XHtoa+DYmdcg/AwwvB1tZKkmg2lcGjX+GcYo9sGH9AR9gn33AwMa7foZDZAPJawD1PQcNGphOJpLxHvoXj34PFGyHPA+8eC50bm061+1Q8DbisK84azzBgw/JiePpHw6FEantlARRXRgaRW2Q+NtJkIkkEMz6HB++EklJYv4JFCfUAACAASURBVB5uvtJ0IpG08Ezces7CSngrSZfqqXgacGk36JyLM+sZAKrgoe+0obwkkLIA3DK15rGOeXBQOyNxJIGsX1tzvE4bEou4oU2DWuMkXWqv4mmAx4LLeuCcaredx0+b4JlFhoOJRP17ARRV4VxQFNm/85GRYFmGg4lxg4+C9p1i43MuNhZFJF0UVkC/ptA+B/Iz4fL94bxeplPtGRVPQ/7YG9rkENtQPgj3zIbyoOFgIrYN/zeF2IbxfuiQDyOSbLM4qR9eLzRq5rxxbtYGLrrGdCKRlBYMw9BX4J+zYFUhtMyEh4Y4k1jJSMXTEMuCm/vhFM+IVSXwrwXGIok47poVeQcUt2H8fUdCtt9wMEkIj/0T5n7jrFFfswbuvNF0IpGUtrwQ5hXExj9shCVbzOXZWyqeBo3tDa1rrdG4azas3WomjwjlQfjbbJzCGdkwft8WcEoPw8EkYZSW7HgsInWqZS7kZcXGDTO2Xe+ZTFQ8Dau+x2pke6WirXDv1yYTSVo7fzJUhohtGO+DfwyFDK/hYJIwLrgcmjR1Xmdlw9gbzOYRSXH/+wUObA4dcuGQNjBxFDTNNp1qz6l4GnZ8Z2cTWII4p67C8PgcmL/ecDBJP99vhneW4hROAA8M6wQnaG2nxNmyBcoDztr0/Nawb2/TiURS1uSlcPqbMOkXWLkF+jeDIR1Np9o7Kp4J4K1jwbJximcAwlVwxSfaXklcNmpyZM2xB8gAyw/3HmY4lCScu2+B4hLnLM2ypfD046YTiaSsaSudv2rx42Sn4pkA9m0KJ+6DM+sJEIYZK2D8PJOpJK08ugh+LMKZ7fQAFlzWFw5saTiYJJzaW2ppiy2RenNIra2Tf5cCWylbtq2bNSaC8gA0ewTKKqi+0r1hBiy43Nm3W6TeFFZBy9ehqgpn2h1okgkLR0HrHKPRJAF9MxtOPApKiqF7T/jkS2jWzHQqkZSzphjGToT566BxDhzbA249ArJ8O//aRKYZzwSR7Yf7huKUzsiFRiVb4eoPDQeT1HfsZKgKU30xEV545FCVTtm+L7+ETcVQBZRUOPu+ikidO/dNeG8RLNsC81bD8E7JXzpBxTOhjB0A/VoRu497GCYuhOfnGA4mqevZn+GrDXEHPDC8LZyrC4rkNzz6YOz1iuUw4V1jUURS2eINOx4nKxXPBPP6KMj1EbuPewCufQ9+2Wg4mKSepaVwWa29u7I8MH6Q1u3Jb4tupRSVn28mh0iKO2nf2OvcDDgyReYDVDwTTPd8uOEwYnc0CkNxKYx5WWe0pA5tDcKx0yJ/zvw4m8Vb8NSh0CGJdyaW+vevcU759HjgjLPg5FNNJxJJOTf9F174CvK8MKYvzLwUejQ3napuqHgmoNuGwaEdcE63R/b3nLkUbnzfcDBJHdctgJ9KqS6c+GFEexjTxXAwSXgPPAgFm6EiDB9PhnXrTCcSSSlTfoZ7p0BFEArLYcoS6NPKdKq6o+KZgCwLxp8ODTJqHn9wCkxabCaTpJB/r4JnVuAUzsjFRI2y4IMjwKNT7LIDtg3vvRcbb9wI06ebyyOSgjZu3fE42al4JqgeLeCp0+IO2BAMwsnjYHOK/SEUF31bBJd9j/NXP7Jfp88H7w4Cv/45kJ2wLOjePTb2eKBrV3N5RFLQge2gcxOqd46/9nCjceqcftIksHMGwkW/iwxCzqOsAgY9AFsrTSaTpLQlAGfMh/LoYmEv4Idre8CwFiaTSTK58x7IagW+5jDqbOjXz3QikZQxaTH0vheWFUDHhvD+BfD340ynqlsqngnukVOhT2ucdz6RK90Xr4TR4yAUNhxOksvx38EvFTh/7SN3KBrSHP6xn+FgklTOvQoqgKAX3pgEn31hOpFIyrj5A9ha5bxesRkWrTWbpz6oeCa4nAz48mpokIlzoVGkgL7/Lfzfa4bDSfK4aAnMKMKZ5QTwQOtceLm/1nXKrquocK54qGbBp58biyOSajy1Wpk3BVtaCn5LqadxDnw0FrL8OFe6B4AQPPoxPPo/w+Ek8d22Ap4roMYV7BkZ8PFAaJNlOJwklawsaJwdd8CGU443FkcklRQUw+Auzp0MAQa2j1tul0JUPJPEYV3hwdOpvpU2IaAKrn8R3phpMJgktifWwp2rIoPI6XW88E5f6NvQYDBJWs8/BpkZkX08T4SB/U0nEkl6G0rhwAfg/ilQHoDzD4KZ1zkbjqQaFc8kcvlQuOEYnNPtkXu6B6vg7Idh8nzD4STxvLIRrloedyCyfdK9XWFkM0OhJOldfStUVkE4DG+8D9Nnm04kkvQmLYZVhbHxW/PA5/3tz09mKp5J5h9nwFkHRwaRK91DVXDcXfDlQpPJJKFMLIILl4Nda/3m5a3gz+2MRJIUYNuwvtb9ewtS5AbSIga1zdvxOJWoeCYZjwdeugKG9SJ2pXsQAhVw1F/hM818yoRiOHUZVNpUn1rHA6c1g0dT5Ga/YoZlweXnxcbdusDwFNtkUMRlxeUwYR70aQFNsqBvG3j9fNOp6o+KZxLyeuCjm+DgbjhrPiPbMlYG4KS7YNJck+nEqIklcOrKyD3Yo+dpvHBkU3ilO/h0BbvspYMHgicT8EKDxs4NCERkj53zPDz0GSxYDaVl8O9zoF8Kn5hS8UxSGT748C8wsNZNQ0rL4bhb4T3dxS79vFAMp6yGcLRcegA/HNoYJvWEDP11lzpw4z8gbANemLsQXn/fdCKRpDZ9aex1IASzlxuL4gr9JEpi+Q3h07ugT8fIgchWS8FKOPUO+Le2WkofDxTDhesjux5EtkzCC0c0gkndzGaT1OKvNcPp95vJIZIiBnWJvfZ74aBOxqK4QsUzyTXKga8fhv5dqLHVkl0F5/8N7njBYDipfyEbrtoCN2yJLLmIrOfEgiEN4f1OkKO/5lKHHrsDsjKd14f0hzO1j6fInqgKwtnjYOpCaJ0DJ/aBD6+A/VP4NDuoeKaETD9MvReG9KXGVksE4fZn4JzboUL3dk89VTacXAiPb8WZ5YwUTnwwMg8+7QgNU3Q/DjFn0TKoCAN++HoxzP7OdCKRpPTEZ/DqbCithLWF4A3DUfuaTlX/VDxTRKMc+OhuOCF6l4Nw7PHKRzDsClijXU9Sx6Yw9N8IH1QS+2vscx5jGsOE1roVptSPtz6m+i5YoRBMmGw6kUhSKije8ThVqXimkKwMmHgnjD2ZGlstEYKZ86DfWfCN9vpMfp8GoNtm+CEUOeCh+q5Ef8qDF5qBV6VT6km3TjXH3TsbiSGSzGwberWGnMiSacuCy4YYjeQay7Zt23QIqXv3vwp/fBTnvu525BlnH9Cn/goXneb8QZfkYds23FMJd5RhBSN3D4jKBB5qBJfnmIon6WLTFhg0Glb8Cl07wYzXoGED06lEkoZtwxmPw1tfO+PBPeGfZ8FBXXb8dalCM54p6obRMPGfkJtNrJ+EIVwJl9wCp13tnCWT5GCvDmMfXQ63VGIFIbaeE2hgwadNVTrFHf95H35cBRU2fL8MbnrYdCKRpLJwdax0AnyxOLXvVFSbimcKO+EwmPYsdO0QORBXNN+dDG2OgK/mGYkmuyH0fohw3wqYFAIs7Oj6OnzQIxOWNINBGYZTStr4aXmt8QojMUSSVW5mzbHXA1lptCuZimeK69cDvnsDjj2s1gdsWF8Ag8+Cq27X7GcisjfbBC8OEj4xCJuJFE4ADzYe7NGZ8HUetNKV6+Kik4Y5a3aiThluLotIkqmogv98AQPbg2WDzwuPjXH25U4XKp5pIDsL/vsveOr2uHWdISAMgUp4fDzsOwymzjQYUmoITrap2s8m/KyzOauzENtDGA92tgeeysZ6pQE01EJdcdmwg2FAX8ADngzI0RIPkV11yTNw25vwzS/O9klf3ASXp9l7NxXPNHLpmTD/PejQmppXvduwZAkcOQrOvwoKi8zmTGf2eqg4E6qOAgogum2NjceZ8eztxZqZi3Vp5o5/IZH68sGX8PUPgM/Zsu26B0wnEkkaU+J2lgmGYPbP5rKYouKZZvp0h3nvwkVnxh0MAyEIBeDfr0DPg+GxcaYSpie7HKqehLKOEHoj7jge5z2C34N1hR/P/Cys/fXXVgyqvR2GtscQ2WUD4q5ctyw4IA13I9NPsDTUpDE8cw988Cy0aUX0PG716feCdXDNDXDAIPjsc4NB04AdhqpPoGg/qLwCqIgcx3k/YGNht/bhm+3H+4QPS5vCi2nHHw5HR+5U4fHAXy8ym0ckCRRthaG3wsSZ0DwbjuwNr14Fh6fBnYpqU/FMYyOHwfeT4OJzI9cKRAqoFXIWPc+fB8ccByOOgyVpeDqgvlXMgM2HQdExYC+v/VGLMB6sWywy13jw9FPhlATh80FxJc6pdg/88xXYrPU5Ijvyt7dh6vfOHp4biqBdYzjzUNOpzFDxTHNN8mDcfTD5TRjQD7Aju0PaYIWdGbkpU2D//nDq6bBBt93ca+VfQsFxsHEQBOMu6AoRm3y2+kP2Ssi8U4VTEkxhCcxcQPVtM9duhHk/mU4lktC2bN3xOJ2oeAoAQw+D2ZPh4XuheXNiBTQiEIAP3oP2bWHksfD9Auedm+y6rZ/Cr8fAmiOg8qPY8TDR0+oQbAIZ70D21+BpbyioyI40bgAdWsXGWZnQVX9YRX7Lmk2QnwNZkdtjZvrhymPNZjJJt8yUbRQWwm13w/PjobzcORY9/U7YmQm1gIED4IabYORJztk32b6i92H9bRCcW31Hdfw4/w29kQceaHAX5FwLlnankUT3xVw46XooKoUDesKUJ6FhrulUIgln7WYYcK3zDHDiIXD/hdCtjdlcJmnGU7aRlweP3A9zv4FzRjsb3Fav/wxX3zOH7+bAeaPggC7wt1ugtNRg6ART8ROsvQMWNIblJ0LlXOd49F1eZBcrQh7wXwYtiiH3ZpVOSRLPTYDCUucP8ZzFcN9/TCcSSUgfz4mVToAvv0/v0gkqnrID+3SBF8bDzOlwzrk1b1YSfekD1qyCR+6GPs1h1FCY+gkUFZpIbJZtQ8F4WHQUfN8DCm6HcHHkY5HPiWwcgJ0NWedD21Jo+iRYmiySZLKpaMdjEQGgXbOa4/bNtv956USn2mWXLfwexj0F459yTr17iJ06jp5pj55Czm0AI06B086Dw440lbj+BQqhaAasfRpKPgeKYv8Nat961wP4m0DjM6Dl/eBp4Hpckbrx0Qw46QYIBCErA74YBwfuZzqVSEJ5bybc+CJsKYFQGDq1hBevg14dTSczS8VTdltBAbzyAjz4Nygvjq1TrF6vGBE9JZ+VBYOGw9GnwfBToFGegdB1qGwVbJ4Ga16OlM1S5/ve3vrN6MywvwM0/QM0vx68mt2UZLepCPqMgbUbAA+cdyy8+FfTqUQSxrot0OlCqAw44wbZsOZFaKjlVCqesnf+9yE8di/M+woI1iye0SIWLaQWkJ0BrdvD4JHQ91AYeir4ak8NJpitK6BkCax+HTZOhbKfne/Nx7alO/p9+iLPjQZDi2ug8QlgJfj3KbLL3vkcTvtLbOz1QtWUmutxRNLYvKXQ/+qax5aPh44tzORJJCqeUieWLIIJr8Lrz8CGdc6x+FPO8c/RQurB2VYiLx8OOAK69oMDjoRmbSG/tZFvg1AA1kyGwh9g/VQoXQjly2LZYdslBvFF0wvk9oDm50Lry8GXb+CbEKlvcxbDgRfH9lTr1BqWvWk2k0iCsG14dzpc/yws3wBYcHgvmPp3vTcDFU+pY4EALFsCzz0IX02GtStipS2DbYso1Cx1AH4PNGkOzVpDlz7Qqgu06Qatu0JuHjRtC5m5u3+L6HAIKouhqgiKV0LRMudRuBBKV0LxIgiXOOVxe5mj4otmtIBmNId2v4fmp0Dj3+1eLpGkdOuz8I+XIRiGUcPgtVt133YR4KIH4bmPnddt8uHGs+GiEZCdaTZXolDxlHo1ayrMngIfvgSbVjkzin5iayHjL0yCbWdG45+jxzP9QBAaNQF/JmTkQEYW+DzOsxVwfl27EuwqsMshUALeIITKwW9RfZl57ZnL+MJJ3HPtcmxZkD8Qmg+HTmMhqw1Yeicr6eTQK2HmD7Hxq3+Fs4aZyyOSACqrIPuEmjdY+fhvcPRAc5kSjbb9lnp18BDncdUdUFkBH/4H5k+Drz+B8s0QDjgdMFrq4k9nW3HP3rhnO+A8l2+GyP72211rGf/w4fw+PojtbUTNMhkvHPm9Q5GvDwJ53aHJQOhwLjTqBTkd9vy/i0jSK9iy47FIGsrwQ5MGsLkkdqxlE3N5EpFmPMWYslKYNgGWzocfpsOaJVC4oebp7fiLk7b3iIqfudyV56jorxP9fSycfzh8WdByAOT3h46nQm4H5yEiEfe/Dn982nndKAfmP+Os9RRJU4tWwJh7YNlaCFng8cNfRsP1o0wnSywqnpIwwmFYvxKKN8H3U2HNj1DwC2xaDoW/gs+CcOW2p+JhxyUzWiyr15Z6gDBkNYHMJtC0IzTt4RTNBm2gzRHg15ZHIju2agP0uhhKtgIeuPBoeO7/TKcSMWbgJTDnp9j4g7/DSK3534ZOtUvC8HigVSfn0X3Ath+vqoD1v4AdhoKfANsppOEglBc5t/MMB8BrOWsws3KcXV4atHCem+7j3G++eR/n49lNXf4GRVLJZ/OgpILqcw9vTVPxlLT264aa49UbzeRIdCqekjQysqBdL+d1+z5ms4ikva5tdjwWSSOvTYIebaFgE+Bx1nWOPMR0qsSk63BFRGT3DeoFfz4DfJHV0RlZUF5pOpWI6+4cD6Nvgy++dZZ03Xw2zBkHbZubTpaYVDxFRGTPTP8JgpGdb7/6CZ7+2HQiEde9+VnsdSAIWX6Vzh1R8RQRkT1TVrnjsUiKq6iMlMy4y7S7tjUWJymoeIqIyJ655UzwRy4uys2CoX3N5hFx0dzF0PEE+ORLaJwFPTvAbRfC6BGmkyU2FU8REdkzh+0H2bmAD7aGYPQjzrlGkTTw58dg/WbndVGJc+Ou2y8ymykZqHiKiMie+XkdFFdQfY+xFRtgQ7HpVCKuqP0eqypgJkeyUfEUEZE907MttMqLjTu3gJaNzeURccHWMrjkDigoAL8F2NC5LYw9w3Sy5KDiKSIieyYvF568FHyRG9wuL4Q3Z5lOJVKvbngAnnkbFi2FQCXcdD589yq00ZXsu0TFU0RE9twXiyEYuYmtDTz92c6+QiSp/bC05nhrGTTIMZMlGal4iojInqt9ar2VTrVL6nr7f5DhAcLO2LJg5GCjkZKOiqeIiOy5a4+FofvFxis3QkWVuTwi9eShF2DUNTB5BnhtGH0MTHoaRhxqOllyUfEUEZE9l+mHwrLYeMYSeOELc3lE6snbk2KvQ2FonQ/DdT/23abiKSIie6ey1j4yFdpXRlLL0lXbruPs3slIlKSn4ikiInvnrtNjdzDK8ENeA7N5ROrQe5Oh50j4ZBrkZEPPzvDni+ASbZ+0R1Q8RURk7xzbD7JyAB9UAZe8AKs2GQ4lUjf+Ng4CkUn8snI4+3i493rnwiLZfSqeIiKydwrLoKSS6jsYBUKwptB0KpG9FghAhq/mMW2dtHdUPEVEZO+0agzD465sz8uFZg3N5RGpA5/PghYHwfRZkOkBbBh2CFx6pulkyU3FU0RE9o5lwRNjICMT8EJhFZz2lOlUIntl7G1QWOy8rqyAJ2+FT5931nnKnlPxFBGRvbdkA1TZgBewYP6v217tLpJEyitrjm3bTI5Uo+IpIiJ7r397aJQVG3fMB69+xEjymfENdDscCtaAFblDUa9ucPaJZnOlCv2rICIie69tExh3Hni8gAdWFMPYN0ynEtltZ1wOPy937sFuB+GJ22H2u9BYy5brhG/nnyIiIrIL1pZA2Bsbf/yDuSwie8C2YX2tncBaN9O6zrqkGU8REakbfdrUHDfXFJEkjynToP9waJgBRE6xd+sMwwYZjZVyVDxFRKRuDO8BN40AywN4YM5auH+K6VQiO1W6FU7+PcxfCJu3QIYFj90JsyZC40am06UWFU8REak7tgdsH85KLgvenG86kchObdoMxSWxcVUADjsQmuSZy5SqVDxFRKTu7NOs5jgnw0wOkV302ttw8jmRTRkip9j37wX7djMaK2VZtq2dqUREpI7YNpz7MrwyN3LAC+NOg4sPNhpLZHsW/wS9D4VQyBk3yYPbb4bzz9Qp9vqiGU8REak7lgX5DXFOtUdOt7+m0+2SmJatiJVOgC2FcMkYlc76pO2URESkbnVuWnPs0RyHJJ5b74JX34RsH5QHAAuOPxqysnb6pbIXdKpdRETqVjAEp70EExdFDvjg8RNgrE63S2J4/W04+/excfducOXlcNF5kJlpLFZa0NtQERGpWz4vdGoG+CMPC95aaDiUSMyyZTXHXg+MvVil0w061S4iInWva37NcVUIwmGddhejNm+Gc86Dr2Y5BShoAxaccZrpZOlD/wKIiEjdG3sQnN4rNp6xEm6ebC6PCHDrHTD5UygthXAQhh4Gb7wEt95kOln6UPEUEZG65/FAt1qznpN+MZNFBKiqglUrax7r2AFOO8lMnnSl4ikiIvWjf+ua45Iq2FRmJouktUmToVU7+O/Hzo5f4KznPH+M2VzpSMVTRETqx6jecN2hgAV4YEkRjJlgOpWkoSuuhJLILTFtGy65GGbNgMGHm82VjnRxkYiI1J99ole3R8xdZyyKpJ9wGKZNg+KimscPORh67WcmU7rTjKeIiNSfwR0g0xsblwTg0+XG4kj6sG0YfRYcNRwKN4Inch/2fv3gFK3rNEYbyIuISP16aQGMmUj1Kfe8LNh4tbN5okg9Wfg9HNC/5rHX3oBjj9XdiUzSqXYREalfbRpR48dNYSVUBCE3w1gkSW2zvoL/PA+WDTaABV4vDB6s0mma3m6KiEj9+l0b2L9F3AEPXDjJORcqUsfmfANHDYbnxoEVcoqO3w+PPAr5+Tv9cqlnKp4iIlK/sv3w8emAF2fm0wtv/ASzdaGR1L1P/weBQGzcvStsKYKLLzGXSWJUPEVEpP41ygSvD+fHTmQjxUVbTCaSFLNpI4w+CZ55zHmLE7Vfb2fGUxKDiqeIiNS/HD88OrS6c4IHLvkcvi4wmUpSyJ+uhv9OhPXrnOLZpTOceTY8Mc50Momn4ikiIu64Yn/o3hxnX08fBMLwySrTqSQFfDcXfviu5rHTz4QXXoamTc1kku1T8RQREff0bUbctCe88gus2WosjiS/m6+BYQfAkoWxU+wZGXDCKUZjyW/QPp4iIuKezRUw4n2YswFn7sMLozrDm0eaTiZJaMN62LdlzWMXXwOjfw99+xmJJDuhGU8REXFP0ywY3gHndHtkfmrBFgiGTaaSJDThdbjuAvBZNY//4QqVzkSm4ikiIu4a3RVy4zaU/7EERk01FkeSz/QpcMVo+Oy/4LFj5fOGW6Frd7PZZMdUPEVExF398uHfQ6mxr+eEVbC2zHAwSQYbCuDjd2vef6BTZ1hWDDfeYS6X7BoVTxERcd+AZnH7euJcb3TXd7qbkezQlI/g8M7wn8ecxRpRhw6Bhg1NpZLdoYuLRETEjOeXwGVfQVXc+s4XD4Pz9jGXSRLaCQNhwZzYuNcAGHYCXHkjZGaayyW7TjOeIiJixgXdoFujmscmrNKsp2xjxhQ4cwgsW1Lz+CnnwPW3qXQmExVPEREx54KuNcfvrIGr52z/cyUtbSiAi06AWZ9DaTFYkQuJ+h8Coy82m012n4qniIiYc30vuHI/nAuN/IAFr6/UrKcAsPBbeOY+KI+7x4Btw4TZ8PZ0yG1gLpvsGRVPEREx68R2OMUzMpW1oQpGTtfenmnuqylwxiHwwoPVb0kA6DPAeXjUYJKS/reJiIhZR7WGB/uDz4NTL3zwUQG8t8Z0MjFk8wZ47SkIBGLH9ukG/3cnvDRZpTOZ6ap2ERFJDC0+gA2VsXGvPPjwd9Axx1wmcd0bT8M9YyEUghDOA+Dsy+H2f5lMJnVBxVNERBLDSyvh999AyMaZ+fTCkGYwZZDpZOKSYBAOzIVAVexYg3zoPwj+8SI0yjOXTeqGJqtFRCQxnNsBbt2P6rsZYcHXhfDhesPBpL7ZNtx3NYxoDVTV/Nijb8KTE1Q6U4WKp4iIJI6z20HjuEtJttpw0lyYXWg0ltSvj16B1x6Dwo3Vq3wBOPJUGHiEyWRS13SqXUREEsuSUhg4HYpDVM+PnNwSXu0LWV6j0aRubVoHD14DC2fBrysguo9B2y7w8PvQZd/Yvp2SGjTjKSIiiaVbAzi6JTV+RL23CU773lgkqR+3nQOfvgHrVjiznNGOefrlsM9+Kp2pSDOeIiKSeEqDcP4CeGcDTgGNlNCpB8ARWuyX7Gb+FyY+DbP+B2UVseMjzoGRF8DBw81lk/ql4ikiIolp8VboMxuCcT+mPBnw3r5wQr65XLJXflkAFx0AoaAzDgNBICsHXpgDnXqaTCf1TafaRUQkMfXMhVd6QdPoxUZ+p6XcthLWVO7kiyXRVJTBU9fDfeeDHYwdz8iAc/8ET09T6UwHmvEUEZHE9udlcN/qmsfaZcDc/tDMbyaT7Lb7fg//ezE2DgA2cMhxcN+HhkKJ61Q8RUQksW0NwaU/w8sbah4f0xKe6gI5utI9kU2fAG8/BEvmQHlp7Pg+B8CAo+HcmyCnobl84i4VTxERSQ77fANL465EwQfD82Cyzs8mqhWL4JK+sfWcUR4P/PMz2F97dKYdrfEUEZHk8O6+0CM7MvACHvi0GP60GirDO/pKcdnG1fDnIXDDoThXDkVZMOI8uHOCSme60oyniIgkj58rYL8FEIj/0ZUBVzWDR9sZiyUxm1bDIxfDNx/FjgVxrgsbeDTc+7GpZJIIVDxFRCS5TNwCV6yA1QGqtx33WXBWM3iuDWRo13ETbBseOR+m/scZR8smQPffwUEnwKnXONsmSfpS8RQRkeTzYREcvyzugAfwwVVNGAIlogAACHVJREFU4f6WKp8uW/AZfD0RJj4SO2bjXLme3QDunwGd+5hKJ4lExVNERJLTm4Vw9WpYF8YpnpbzfFAuTGkLObqMwQ2f/AueHeu8tokt6fT44IKHof9R0K67qXSSaFQ8RUQkeb1cBGNWO40HcE69e+DIXHiiGXTXPp/15et34MMHYNVCKC6KHQ8BtgUXPAAnXmcsniQoFU8REUlu08vgxNWw2Sa2WYsXmvthXito4zOZLuVs/hV+/AKePA9CIedY9LQ6wGm3wNALoWUnQwEloal4iohI8vvvVjhjHWy1cU65Ry468vjgocZwda7hgKlh3vvwxOkQrKxZNgGyW8JBp8BFTzj7dIpsj4qniIikhtIQ9FoHK0NxB/1gWXBVNtzWEJqqEe2Jue/ClMdgxTwo3hI7Hr1yvd9x8Efd9lJ2gYqniIikjl8CcNUW+KiS6k3mo/b1wodNoLNOve+qX2bCks/hvb+AHdkbKUzsAqKeR0Lvo+GosZCR/Vu/ikiMiqeIiKSeswvh1YrtfMAP/8iBP6kl7UigEr4cB69fvf2PVwE9BsN1H0Gm9uWU3aDiKSIiqce2YUIFnF0E5fEf8GNjQR8/jMvGOkSzn/FsG16+EGa/CFhQFY7bMCCiz0g4bzw0bO6sYhDZHSqeIiKSur6sgj+WwKwg0QuOnB96XuwsC87PxLrVj9Umvdd+Fq+Fd66EtQtg3ZLYHYfiLyDqdQx0PRyOvE6n1WXPqXiKiEjq+0MpjK+KlE4P0bWfYbzQwgM3+PFc6sVqlF5TeEsmwbyX4eepsHFF7HgAp3RaXghnwEFnw7njdLW67D0VTxERSQv2ohAcXwZLnR97NhY2XiBSQDt7sP7oxXuOldIFNFgFsx6DdfPh21cgHNkEIBR5QOxq9VMfgqHXmskpqUnFU0RE0oa9KQzjAoTvCkC5RfTUu1NALaeAdrSwrrHwnWXhaW04cB36dSYsfhuWTYVf5zjH4k+lR19n5MKJD0P7A6Hd/kaiSgpT8RQRkbRjfxsmfE8Ae2IYO+jc5z2+gAaxoImF7zzwnQXeQwwH3kOVRTDzbtj0Iyz5CEKRfZC2N7vZan/oNhL6nQFtVDilnqh4iohI2gpPChO6MYQ938YOOaXTBkI4s6EAeME+BLzHQdaNYCXwOkc7DKEq+PZx+PF1KF4FZQWRj7Ht7CYWtBsMLfvCiLshq5GR2JJGVDxFRCTthT+zCY61Cf8M4aDlbLkUJzorGPID3SH7OvB1h8zDTaStqXQ1/PoplK2HWXdBVXHNGc140e8jtw002gcGXAL7n+tuXklvKp4iIiIR4e+g8g8Q/gEoix2PFrZw3Osg4OkOnm7Q4FzwdYKMnuDNq798gRLY8BXggYWPQ+kq2PwjVJY6Hw8R2wopemV6lOWDpn2hWS848iHIya+/nCK/RcVTRESkFnszVF4FoekQXBErc9HSGX/auioyDgI0Aqs9+DtBziHgyYaGx0CwELJ7AxZ4c5wSaNvOBuzR50Cxs31R1UbYuhR8ebB2AuCBdZ9D2UqoKISyjdtmCUaeq3MQK57N+0GTfWD/S6HTUfX6n01kp1Q8RUREdqDyOQh8AhWfQKjYOba9Alp9Op6ap7ltC4I2WNlQWQ6eHAhlQLAcfK1g6yrIaAWlawAPhDyRi4A8VDfe+F8zWii393tHf/8m+0LH46FpT+h9ge4wJIlDxVNERGQXhNZB1WSomAalzwFhCIa3LaDbW18ZnZGMfiz+8wM7+FjUb5XMIODxO0W2UVcYeDv4G0Cbw8CXVSfftkidUvEUERHZTXYQwmVQeB8EfoHSqRBYB7YHAuFan0vs9PeOimc47vOqav1+0Y9ZPqgMQmYz6DDaOW3fbQw06VXX36FI/VDxFBER2Ut2AIIbILgRit6CcIVTRkNFzprNwGYgG6rKnc+vPlUfLaqRZxsgF6q2Qm438OVCZktofjSEq6DtKZDTDrzZOn0uyUnFU0REpB7ZQQhuAqsBbJ0D/hZQuTZywVEPKP0BGvSGYAlYfmjUFwJFkNHEdHKRuqfiKSIiIiKuSOD7L4iIiIhIKlHxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFyh4ikiIiIirlDxFBERERFXqHiKiIiIiCtUPEVERETEFSqeIiIiIuIKFU8RERERcYWKp4iIiIi4QsVTRERERFzx/xWfPSzvKBJtAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31df23cd0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x31e09f210>"
      ]
     },
     "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": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "σ = [3.256607724682357 1.842414551168056]\n",
      "U = [0.28825607506522444 0.9558374294215863\n",
      " 0.9575533589247086 -0.29389591442674756]\n"
     ]
    }
   ],
   "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 ];  @show(σ);\n",
    "U = [ Ax[:,k1]/σ1 Ax[:,k2]/σ2 ];  @show(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.5",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

TB-Lecture-05.ipynb

TB-Lecture-06.ipynb

TB-Lecture-07.ipynb

TB-Lecture-08.ipynb

TB-Lecture-09.ipynb

TB-Lecture-10.ipynb

TB-Lecture-11.ipynb

TB-Lecture-12.ipynb