{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"interact-js-shim\">\n",
       "    <script charset=\"utf-8\">\n",
       "(function (IPython, $, _, MathJax) {\n",
       "    $.event.special.destroyed = {\n",
       "\tremove: function(o) {\n",
       "\t    if (o.handler) {\n",
       "\t\to.handler.apply(this, arguments)\n",
       "\t    }\n",
       "\t}\n",
       "    }\n",
       "\n",
       "    var OutputArea = IPython.version >= \"4.0.0\" ? require(\"notebook/js/outputarea\").OutputArea : IPython.OutputArea;\n",
       "\n",
       "    var redrawValue = function (container, type, val) {\n",
       "\tvar selector = $(\"<div/>\");\n",
       "\tvar oa = new OutputArea(_.extend(selector, {\n",
       "\t    selector: selector,\n",
       "\t    prompt_area: true,\n",
       "\t    events: IPython.events,\n",
       "\t    keyboard_manager: IPython.keyboard_manager\n",
       "\t})); // Hack to work with IPython 2.1.0\n",
       "\n",
       "\tswitch (type) {\n",
       "\tcase \"image/png\":\n",
       "            var _src = 'data:' + type + ';base64,' + val;\n",
       "\t    $(container).find(\"img\").attr('src', _src);\n",
       "\t    break;\n",
       "\tcase \"text/latex\":\n",
       "\t\tif (MathJax){\n",
       "\t\t\tvar math = MathJax.Hub.getAllJax(container)[0];\n",
       "\t\t\tMathJax.Hub.Queue([\"Text\", math, val.replace(/^\\${1,2}|\\${1,2}$/g, '')]);\n",
       "\t\t\tbreak;\n",
       "\t\t}\n",
       "\tdefault:\n",
       "\t    var toinsert = OutputArea.append_map[type].apply(\n",
       "\t\toa, [val, {}, selector]\n",
       "\t    );\n",
       "\t    $(container).empty().append(toinsert.contents());\n",
       "\t    selector.remove();\n",
       "\t}\n",
       "    }\n",
       "\n",
       "\n",
       "    $(document).ready(function() {\n",
       "\tfunction initComm(evt, data) {\n",
       "\t    var comm_manager = data.kernel.comm_manager;\n",
       "        //_.extend(comm_manager.targets, require(\"widgets/js/widget\"))\n",
       "\t    comm_manager.register_target(\"Signal\", function (comm) {\n",
       "            comm.on_msg(function (msg) {\n",
       "                var val = msg.content.data.value;\n",
       "                $(\".signal-\" + comm.comm_id).each(function() {\n",
       "                var type = $(this).data(\"type\");\n",
       "                if (typeof(val[type]) !== \"undefined\" && val[type] !== null) {\n",
       "                    redrawValue(this, type, val[type], type);\n",
       "                }\n",
       "                });\n",
       "                delete val;\n",
       "                delete msg.content.data.value;\n",
       "            });\n",
       "\t    });\n",
       "\n",
       "\t    // coordingate with Comm and redraw Signals\n",
       "\t    // XXX: Test using Reactive here to improve performance\n",
       "\t    $([IPython.events]).on(\n",
       "\t\t'output_appended.OutputArea', function (event, type, value, md, toinsert) {\n",
       "\t\t    if (md && md.reactive) {\n",
       "                // console.log(md.comm_id);\n",
       "                toinsert.addClass(\"signal-\" + md.comm_id);\n",
       "                toinsert.data(\"type\", type);\n",
       "                // Signal back indicating the mimetype required\n",
       "                var comm_manager = IPython.notebook.kernel.comm_manager;\n",
       "                var comm = comm_manager.comms[md.comm_id];\n",
       "                comm.then(function (c) {\n",
       "                    c.send({action: \"subscribe_mime\",\n",
       "                       mime: type});\n",
       "                    toinsert.bind(\"destroyed\", function() {\n",
       "                        c.send({action: \"unsubscribe_mime\",\n",
       "                               mime: type});\n",
       "                    });\n",
       "                })\n",
       "\t\t    }\n",
       "\t    });\n",
       "\t}\n",
       "\n",
       "\ttry {\n",
       "\t    // try to initialize right away. otherwise, wait on the status_started event.\n",
       "\t    initComm(undefined, IPython.notebook);\n",
       "\t} catch (e) {\n",
       "\t    $([IPython.events]).on('kernel_created.Kernel kernel_created.Session', initComm);\n",
       "\t}\n",
       "    });\n",
       "})(IPython, jQuery, _, MathJax);\n",
       "</script>\n",
       "    <script>\n",
       "        window.interactLoadedFlag = true\n",
       "       $(\"#interact-js-shim\").bind(\"destroyed\", function () {\n",
       "           if (window.interactLoadedFlag) {\n",
       "               console.warn(\"JavaScript required by Interact will be removed if you remove this cell or run using Interact more than once.\")\n",
       "           }\n",
       "       })\n",
       "       $([IPython.events]).on(\"kernel_starting.Kernel kernel_restarting.Kernel\", function () { window.interactLoadedFlag = false })\n",
       "   </script>\n",
       "</div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using PyPlot, Interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection onto a line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose $b$ is a vector of data and we want to find $p$, a multiple of $a=(1,1,\\ldots,1)$, say, closest to $b$.\n",
    "Which vector is that?\n",
    "\n",
    "Let us call this vector $p=\\hat{x}a$.\n",
    "\n",
    "Here is an example in 2d:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.405977\n",
       " 0.752382"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = rand(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAGwCAYAAADfbKDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XlcVOX+B/DPLDAssqgoi6K4r4kriNY1jaI0b7fyillqmppmZtG1tLqa93bzV5mtmklulXtldW+Ger2ZuSCK4gLuqKCyCMoOAzPn/P44ipJgzDAzzyyf9+s1r+dxOjPz4Wjz5TnnPM9RybIsg4iIyEmpRQcgIiKyJhY6IiJyaix0RETk1FjoiIjIqbHQERGRU2OhIyIip8ZCR0RETo2FjoiInBoLHREROTUWOiIicmosdERE5NS0ogPUhyRJuHz5Mnx8fKBSqUTHISIiAWRZRnFxMUJCQqBW13+c5hCF7vLlywgNDRUdg4iI7EBmZiZatmxZ7+0dotD5+PgAUH44X19fwWmIiMiWZFnGs18nY1dqJi599nR1TagvlSPcpqeoqAh+fn4oLCxkoSMickHllUas+jUNUx/oYXIt4MUoRERkl/5z5DKuFOsBAJ7uGoyObG3W+7DQERGR3VmblIHn1xzC6PhEFFVUNei9WOiIiMiunMktwWubjgIA7u4QAB9dwy4ncYiLUYiIyHW0b94I/3ikO85dKcXfH+7S4GllHNEREZHFJSQk4O6774a/vz+aNm2Khx9+GGfPnr3jay4VlFf3x/RvjTnDu1pk7jQLHRERWVxpaSni4uJw4MABbN++HWq1Go8++igkSap1+6/2nsfg93bgv2k5Fs/CQ5dERGRxjz/+eI0/L1++HM2aNUNaWhq6d+9e47+t35+Bv/+QCgBIzriG6K6BFs3CER0REVnc6dOn8cQTT6Bt27bw9fVFWFgYACAjI+O2bfuFNUFzHx2eHdQWr8R0sngWjuiIiMjihg8fjtatWyM+Ph4hISGQJAndu3dHZWVl9TayLEOlUqFts0bYPOMeNPV2t8p6xix0RERkUfn5+Th58iTi4+Nxzz33AAB27dpVY5svfkvH8axivDuiBzRqFQIa6ayWh4WOiIgsqnHjxmjatCmWLl2K4OBgZGRkYNasWdX//Yvf0vHWT8cBAPd3bY4HuwdbNQ/P0RERkUWp1WqsW7cOycnJ6N69O1566SW89957AJTDlckXrgEAXhjSHjHdgqyeh4s6ExGRTVUZJSQcy8bDPYJNOidnbi3giI6IiKzuq8QLyLxaBgBw06gxPDzEZjfSZqEjIiKrWvTLGfz9+2N4wgILNJuDhY6IiKzmTG4x3t96EgAQ2zcUvh5uNs/Aqy6JiMhq2jf3wQexPXHxWjmmDW4vJAMLHRERWdyJ7CJ0DlIuGHmkZwuhWXjokoiILOrD/57CQx/9hk2HLoqOAoAjOiIisqDlu87hw/+eBgDkFukFp1FwREdERBbzYPcgtGrihdeGdsazg9qJjgOAIzoiImogWZZhlGRoNWqE+Hvi5xn3wFtnP+WFIzoiIjKbLMtYsPUkpq4+iEqDclNVeypyAEd0RERkJlmW8d6Wk1i84ywAYOepKxa/aaolcERHRERmkWUgv0S5v9zc4V3tssgBHNEREZGZ1GoV5j92F4aHh+DuDgGi49SJIzoiIqo3WZbx8fbTOJ1TDEApdvZc5AAWOiIiqidZlvGvn45j4bZTGP3FPhQLWKDZHCYXup07d2L48OEICVFusfD999//4Wt27NiB3r17Q6fToX379li5cqU5WYmISKAzuSX4cu8FAMCM+zrAR8ACzeYwudCVlpYiPDwcixYtqtf2586dw7BhwzB48GCkpKTgxRdfxMSJE7FlyxaTwxIRkTgdAn3w+Zg+ePvRu/BU/9ai49Rbg+4wrlKpsGnTJvzlL3+pc5tXX30VP/30E44dO1b93KhRo1BQUICEhIR6fQ7vME5EJIYsy0hMv4qodk1FR7HfO4zv3bsX0dHRNZ6LiYnB3r1763yNXq9HUVFRjQcREdmWLMuY+2MqnohPxMrd50THMZvVC112djYCA2vOrQgMDERRURHKy8trfc38+fPh5+dX/QgNDbV2TCIi+p0lv6bjy70XoFIBXu6OOxvNLq+6nD17NgoLC6sfmZmZoiMREbmckX1bokuwL955vAdG9nPcAYfVS3RQUBBycnJqPJeTkwNfX194enrW+hqdTgedTmftaERE9DuSJKPSKMHDTYOmjXT48fmBcNPY5Zio3qyePioqCtu3b6/x3LZt2xAVFWXtjyYiIhNIkozXvz+K8Sv2o7zSCAAOX+QAMwpdSUkJUlJSkJKSAkCZPpCSkoKMjAwAymHHsWPHVm8/ZcoUpKen45VXXsGJEyewePFibNiwAS+99JKFfgQiImooSZLx2qajWJuUiX3n8pF0/qroSBZjcqE7cOAAevXqhV69egEA4uLi0KtXL8yZMwcAkJWVVV30AKBNmzb46aefsG3bNoSHh+P999/HF198gZiYGAv9CERE1FAqFeDr6Qa1Clg4sicGdWwmOpLFNGgena1wHh0RkfXJsozUy0Xo3sJPdJRa2e08OiIisk9GScabP6biyMUCAMoiIPZa5BqChY6IyAUZJRkzvzmMlXvO4+kV+x1mgWZzsNAREbmgc3klSDiWDY1ahX8+0t1hFmg2h+NOdSciIrO1b+6DVRMikFesx0N3BYuOY1UsdERELsJglLAtLae6sPULayI4kW3w0CURkQswGCXEbTiMqasP4pPtp0XHsSkWOiIiF7B4x1n8ePgytGoVOgb5iI5jUzx0SUTkAsYPDMPuM3l45u42eKBbkOg4NsVCR0TkpKqMEsqrjPD1cIOPhxvWTe4PlUolOpbN8dAlEZETqjJKeGHtITz1xT4Ulitz5FyxyAEsdERETkeSZExfcwg/H8vGiaxipF0uEh1JKBY6IiIno1ar0CPUD+4aNT4f0wdR7ZqKjiQUz9ERETmh5+5tj4fvCkGrpl6iowjHER0RkROoNEiYvvYQ9pzNq36ORU7BQkdE5OD0BiOeW52Mfx++jGmrD6JEbxAdya6w0BERObhL18px4MI16LRqfPxELzTS8azUrbg3iIgcXNtmjfD1M5EoLK/CwPYBouPYHY7oiIgcUEWVEav3XYAsywCA7i38WOTqwBEdEZGDqagyYvJXydh56goyrpZh9kNdREeyaxzRERE5mMU7zmLnqSvwdNPg3o7NRcexexzRERE5mOfubYeT2UUYP7AN+rd17cng9cFCR0TkAMorjSitNCCgkQ4ebhp8Pqav6EgOg4cuiYjsXHmlEc+s2o9RSxNxpVgvOo7DYaEjIrJjkiTjmVX7sedsPrIKynHxWpnoSA6HhY6IyI6p1SoM6xEMH50WXz4TgV6tGouO5HB4jo6IyA7Jslx9/7gnI1vjwW5BaNpIJziVY+KIjojIzpTqDRizLAnb0nKqn2ORMx9HdEREdqRUb8D4FfuRdP4qjmcVYUC7wfDm2pUNwhEdEZEduVpaiQtXS+HjocXyp/uxyFkA9yARkR0JbeKFtZP6o0RvQI+W/qLjOAWO6IiIBCuuqMKn/zsNo6Qs0Ny2WSMWOQviiI6ISKCiiiqMW56EQxkFuFKsx7xHuouO5HQ4oiMiEmjZb+dwKKMAfp5u+GvfUNFxnBJHdEREAj0/pD1yiyvwZGRrdG/hJzqOU2KhIyKyscLyKhRXVKFlYy+4adSY/1gP0ZGcGg9dEhHZUGFZFcYs24fYzxOReZXrVtoCCx0RkY0YJRnjViThyMVClFcpt90h62OhIyKyEY1ahWfuboNmPjqsmRSJzkG+oiO5BJ6jIyKyMoNRglajjCuGh4fgvi7N4eXOr19b4YiOiMiKrpZW4i+Ld+Pb5IvVz7HI2RYLHRGRlVwtrcTo+EQcu1SEd7ecQKme5+REYKEjIrKSSoOE8iojmvnosHpify7QLAj3OhGRlQT5eWDtpP4orzKiXbNGouO4LI7oiIgsKK9Ej3/+Jw2VBgkAEOLvySInGEd0REQWcqVYj9HxiTidW4KySiPmP3aX6EgEjuiIiCxmzb4MnM4tQZCvByb/qa3oOHQdR3RERBYyfUh7VBiMiO0birAAb9Fx6DoWOiKiBsgtqkBRRRXaN/eBWq3Cqw92Fh2JfoeHLomIzJRTVIFRSxMxamkiTucUi45DdWChIyIyg1GSMW55EtLzSqHTauDhphEdierAQkdEZAaNWoVZD3VG2wBvrJvcH6FNvERHojrwHB0RkQnKK43wdFdGb/d2ao6B7QPgpuGYwZ7xb4eIqJ4uF5TjoY92Yvmuc9XPscjZP/4NERHVw+WCcoxamojz+WVYseccynjTVIfBQkdEVA8ebhp4uWvQqokX1k2O4q12HAj/poiI6qGJtzvWTOqPiiojQvw9RcchE3BER0RUh8yrZYhbn1J9mLKJtzuLnAPiiI6IqBaZV8swamkiLhWUQ6tR4d0R4aIjkZk4oiMiqsV/jmThUkE52gR4I+7+TqLjUAOYVegWLVqEsLAweHh4IDIyEklJSXfcfvXq1QgPD4eXlxeCg4MxYcIE5OfnmxWYiMgWpgxqi9eGdsa6yf0R5OchOg41gMmFbv369YiLi8PcuXNx8OBBhIeHIyYmBrm5ubVuv3v3bowdOxbPPPMMUlNTsXHjRiQlJWHSpEkNDk9EZEkX8ktxOLMAAKBSqTD5T+0Q6Msi5+hMLnQLFy7EpEmTMH78eHTt2hVLliyBl5cXli9fXuv2e/fuRVhYGF544QW0adMGd999N5599tk/HAUSEdnS+bxSxH6eiKeW7cOxS4Wi45AFmVToKisrkZycjOjo6JtvoFYjOjoae/furfU1UVFRyMzMxObNmyHLMnJycrBx40YMHTq0zs/R6/UoKiqq8SAishajJOOZVfuRXVSBIF8PjuKcjEmFLi8vD0ajEYGBgTWeDwwMRHZ2dq2vGThwIFavXo3Y2Fi4u7sjKCgI/v7+WLRoUZ2fM3/+fPj5+VU/QkNDTYlJRGQSjVqFdx7vgV6t/LFmUn8089GJjkQWZPWrLtPS0jBjxgzMmTMHycnJSEhIwPnz5zFlypQ6XzN79mwUFhZWPzIzM60dk4hc0NXSSsiyDADoG9YE300dwCLnhEyaRxcQEACNRoOcnJwaz+fk5CAoKKjW18yfPx8DBgzAzJkzAQA9evSAt7c37rnnHrz11lsIDg6+7TU6nQ46Hf+xEZH1nMktwRPxiRgd0Qov3d8RgHIBCjkfk0Z07u7u6NOnD7Zv3179nCRJ2L59O6Kiomp9TVlZGbTamvVUo1FucXHjNykiIls6l1eKJ+ITcaVYjy2p2Vyg2cmZfOgyLi4O8fHxWLVqFY4fP46pU6eitLQU48ePB6Acdhw7dmz19sOHD8e3336Lzz77DOnp6di9ezdeeOEFREREICQkxHI/CRFRPQU0ckdoY090DvLBmkn9uUCzkzP5bzc2NhZXrlzBnDlzkJ2djZ49eyIhIaH6ApWsrCxkZGRUb//000+juLgYn376KV5++WX4+/tjyJAheOeddyz3UxARmcDHww2rJkTAYJTR2NtddByyMpXsAMcPi4qK4Ofnh8LCQvj6+oqOQ87u3nuBnj2BDz8UnYQs6GR2Meb/fBwfxfaCn5eb6DhkBnNrAde6JNeRmQnMmwd06gSsWSM6DdnQyexijI5PxI6TV/CvzWmi45CN8cA0OTeDAfj5Z2DJEqVVqQBJAj77DBg9WnQ6spG9Z/OQX1qJu1r44bWhXUTHIRtjoSPnlJEBLFsGfP45kJMDaDSALCsPAEhOVgqeuo6DGgYD8PzzwFdfAW5uwNSpwD/+oRRKcjhPD2wDL50WMV2DeNjSBfHQJTkPgwH44QfgwQeBsDDgX/9SihwAGI01ty0vB06erPu9Vq0CtFogKQn46CNg4ULgiy+sFp0sL+1yEfacyav+88i+oSxyLoojOnJ8588rRWjpUuDKlZujt98Xt987cADoUsdhrNBQ4IMPlBFcp07A0aPKn3nXDYeQerkQT36xDxVVRqye2B99WjcWHYkE4oiOHNvJk0DbtsDbbytFDvjjAgcohyMPHKj7v/fvX/MwZVQUcPp0/d6bhDJKMqavOYSCsip0DvJFh8BGoiORYCx05Nh8fYFevZQRXF3n22pTVQUkJlovFwmjUavw6ejeGNK5Ob58JgK+Hjxc6epY6MixBQcrI7NvvlHOy5ni8GHlvF5t9u2r+efERKBDB+WwKNmlzKtl1csKdg3xxfKn+7HIEQAWOnIGKhXw+OPKYczPPweaNavf1ZF6PZBWx5yqjAwgLk55z7VrgU8+AWbMsGxuspgjFwsw7OPf8OaPqVxDl27DQkfOQ6sFJk9WLk55+22gUT3OzdR1nm7sWOXKzIgIYNo0pchNnmzRuGQZp3KK8eQX+1BUYUDq5SJUVEmiI5Gd4VWX5Hy8vIBZs4D8fGDBAuU5jeb2C0luXJAyYULN53fsuNn/7DOrRqWGa9XEC71aNUZ5pQErxkfA052Hl6kmrnVJzqmwEPD3V/qHDytz4VasUArereflevYEDh0Sk5EaRJbl6vvHVVQZYZRkeOv4u7sz41qXRLe6Md9txgygRw9llZTUVGDoUOX5GxeVHDsGVFaKyUhmS75wDX9ZvAe5xRUAAA83DYsc1YmFjpzPpUvAxo1Kf/78m8936aKsnLJnjzJPDlBGd0eP2j4jmS35wlWMW56Ew5kFWLj1lOg45AD4KxA5nwcfVNqPPwY8PW//71FRwG+/KYs8JyYCd91l23zUIGdzS1GiNyCqbVPMGd5VdBxyADxHR84lNRXo3l3pGwyc9+aktqXl4O72AbzwxMXwHB0RcLPI/fADi5wTSTp3FVtSs6v/fH/XQBY5qjcWOnIev/xysz98uLgcZFGJ6fl4ekUSpq0+iH3p+aLjkANioSPnIMvAkCFKf/9+3jfOSRiMEl777ijKKo2IatcU4aH+oiORA2KhI+fw9ddKGxgI9O0rNgtZjFajxrKn+2FEn5aIH9sXHm48XEmm48Uo5PiqqgB3d6Wfng60aSM2DzXY8awidAz0gUbNkTndxItRyHXdmCs3eDCLnBPYfSYPf1m0G3/beBhGye5/DycHwEJHjq2kBJg7V+lv2CA2CzVY6uVCTFi5H3qDhMLyKhgkLtBMDccJ4+TYnntOaadMAQICxGahBusU6IPoLoGoqDJi8VO9odPynBw1HM/RkePKyQGCgpR+WVntq6CQQzAYJWg1ygGmKqMESZZZ5Og2PEdHrufGXLkFC1jkHNiOk7l44IOdyLxaBgBw06hZ5MiiWOjIMZ0+rcyXA4AXXxSbhcz2y8lcTP4yGel5pfh851nRcchJsdCRYwoPV9pvvuFSXw6svNIIoywjplsg5jzcTXQcclK8GIUcz549QHm50n/sMbFZqEGG3hWMgEY69GrlDzcNf+8m6+C/LHIssgwMHKj0ExO51JcD+m9aDr5Jvlj954g2TVjkyKo4oiPH8s03SuvjA0RGis1CJtuWloPnVifDIMkI8ffAgHacEkLWx1+jyHEYjcDIkUr/0CGxWchkBqOE97acQJVRxrC7ghER1kR0JHIRHNGR43j/faUdMABo105sFjKZVqPGlxMisXz3ObwS06l63hyRtXHCODmGsjLA21vp5+YCzZqJzUP1lnTuKsJD/Tg3jhqME8bJud2YKzdhAoucA9l8NAtPxCdi2uqDqDRw3UoSg4WO7F9+PhAfr/Q/+URsFqq3IxcLMH3tIRglGT4ebrzlDgnDc3Rk/27MlXv7bcDLS2wWqrfuIX4Y2bcl9FUS3vtrOAsdCcNzdGTfzp+/eY+5qipAy9/N7F1ZpQFe7srfkyTJkAEWObIInqMj59S3r9KuW8ci5wB+SLmEe9/bgVM5xQAAtVrFIkfCsdCR/UpOVs7PATfnz5Hd+iHlEl5an4LcYj02HsgUHYeoGn9FJvt1YzS3ezeX+nIA/l7u0GrUeKxXC8x+qIvoOETVWOjIPv34o9K6uSkTxMnuDerYDD9MG4hOgT5Q83Al2REeuiT7I0nAI48o/dRUsVnojr5Jvohlu85V/7lLsC+LHNkdjujI/ixapLS9egEdOojNQnXaeCATr3x7BLIMdAnywYD2XKCZ7BNHdGRf9HrghReUfkKC2CxUJ4NRwqq95yHLwJj+rRHVrqnoSER14oiO7MvMmUo7ejTQvLnYLFSnGws0bziQiWf/1BYqXixEdowTxsl+FBQAjRsr/ZKSm4s4k93YmpqNge0D4K3j78hke5wwTo4vNlZp33yTRc4OrdmXgclfJWP8yv2oqDKKjkNUbyx0ZB8yM4GtW5X+66+LzUK3SckswGubjgJQ1rDUafnVQY6Dxx/IPkRFKe2XX3KpLzsU3tIPU+9th0qDhDeGdeE5OXIoPEdH4h05AoSHK31J4ioodiS/RI+mjXQAgBtfFSxyJArP0ZHjulHkduxgkbMjq/acx70LduBQxjUASoFjkSNHxEJHYm3ZcrM/aJC4HFTDV4kXMPfHVBRXGLDj5BXRcYgahCdDSBxZBh58UOkfPy42C9XQNdgH3u4ajB0QhhejuToNOTYWOhInPl5pO3dWHmQ3+rRugi0v/Qkt/D15uJIcHg9dkhhVVcCzzyr9HTuERiHFF7+lY8GWk9UXnbRs7MUiR06BIzoS4403lPaxx4DAQLFZCEt3nsXbm08AAAa0a8oFmsmpcERHtldcDLz7rtJftUpsFkKVUcLW1BwAwAv3deACzeR0OKIj2xszRmlnzwYaNRKbheCmUWPlhAhsPpKFkf1CRcchsjizRnSLFi1CWFgYPDw8EBkZiaSkpDtur9fr8frrr6N169bQ6XQICwvD8uXLzQpMDi47G/jhB6U/b57YLC5uw/5MFJRVAgAa6bQscuS0TB7RrV+/HnFxcViyZAkiIyPx4YcfIiYmBidPnkTzOm6rMnLkSOTk5GDZsmVo3749srKyIElSg8OTA7rnHqX94gvAzU1sFhf26f9OY8HWU/gy0RffTBkADzeN6EhEVmPyEmCRkZHo168fPv30UwCAJEkIDQ3F9OnTMWvWrNu2T0hIwKhRo5Ceno4mTZqYFZJLgDmJEyeALl2UPpf6EuZgxjU8tngPAGBmTCdMG9xecCKi+rHJEmCVlZVITk5GdHT0zTdQqxEdHY29e/fW+poff/wRffv2xbvvvosWLVqgY8eO+Nvf/oby8vI6P0ev16OoqKjGg5zAjSL33/+yyAnUu1VjzHqoM4scuQyTDl3m5eXBaDQi8HeXgwcGBuLEiRO1viY9PR27du2Ch4cHNm3ahLy8PDz33HPIz8/HihUran3N/PnzMY/nb5zLrXPl7rtPWAxXlpFfhlZNvQAAUwa1E5yGyHasPr1AkiSoVCqsXr0aERERGDp0KBYuXIhVq1bVOaqbPXs2CgsLqx+ZmZnWjknWJMvA4MFKPzVVbBYX9cG2U7j/g1+x+0ye6ChENmdSoQsICIBGo0FOTk6N53NychAUFFTra4KDg9GiRQv4+flVP9elSxfIsoyLFy/W+hqdTgdfX98aD3JgX32ltGFhQNeuQqO4oiW/nsVH209Db5BwIrtYdBwimzOp0Lm7u6NPnz7Yvn179XOSJGH79u2IunHjzN8ZOHAgLl++jJKSkurnTp06BbVajZYtW5oZmxyGwQCMG6f09+wRm8VFDerYDE283fHGsC545u42ouMQ2ZzJhy7j4uIQHx+PVatW4fjx45g6dSpKS0sxfvx4AMphx7Fjx1ZvP3r0aDRt2hTjx49HWloadu7ciZkzZ2LChAnw9PS03E9C9umf/1TaYcOA4GCxWVyILMvVa1Z2CfbF9rhBmHhPW8GpiMQweR5dbGwsrly5gjlz5iA7Oxs9e/ZEQkJC9QUqWVlZyMjIqN6+UaNG2LZtG6ZPn46+ffuiadOmGDlyJN566y3L/RRkn8rKgH/8Q+mvXSs2iwuRZRnvbjmJkgoD/vFIN6hUKjT2dhcdi0gYk+fRicB5dA7qiSeAdeuAl18GFiwQncYlyLKMdxJOYsmvZwEAayf159qV5DRsMo+OqN7y8pQiBwDz54vN4kIMkowT2cq80zeHd2WRIwIXdSZrGTJEaT/7jEt92ZCbRo0lT/XBzlNX8EC32q+EJnI1HNGR5Z09Cxw9qvRv3FyVrEaWZSz59SxyiyoAAB5uGhY5oluw0JHldeigtAkJXOrLymRZxj//cxz/9/MJPBGfiIoqo+hIRHaHhY4sa+9eZSUUAIiJEZvFBRzKLMDy3ecAAM/c3ZZ3ISCqBc/RkeXIMjBggNI/ckRsFhfRu1VjzH/sLsgyMDqyleg4RHaJhY4sZ+NGpQ0MBO66S2wWJybLMlIvF6F7C2VZvSciWOCI7oSHLskyjEYgNlbpHzwoNosTk2UZc35IxV8W7UbCsWzRcYgcAkd0ZBnvvqu0990HhISIzeLEPth2Cl8lXoBKBRRXVImOQ+QQOKKjhquoAF57Tel/953YLE7usd4t0cLfE++NCMdf+4aKjkPkEDiio4abOlVpp00DuESbxUmSDEmWodWoERbgjf/GDYKnO6+uJKovrnVJDXPtGtCkidLX6wF3Lh5sSZIk47VNR1FeZcTCkT2hUXNeIrkurnVJYjz4oNJ+9BGLnIVJkozZ3x3Fuv2Z+PfhyziYcU10JCKHxEJH5rtwAUhKUvrTp4vN4oQkWUZppQFqFfBBbE/0C2siOhKRQ+I5OjJfly5K++9/c6kvK9Bq1PgwtifGRoUhog2LHJG5OKIj8yQnA+XlSv/hh8VmcSJGScbbm48jI78MgFLsWOSIGoaFjszTt6/ScnK4xRglGX/beBhLd6bjqWX7oDdwgWYiS2ChI9P98IPS+vgAvXqJzeJEDl8swI+HL0OjVmHWQ52h03IKAZEl8BwdmUaSgL/8RemnpYnN4mR6t2qMj0f1gloFPHRXsOg4RE6DhY5M89FHSjtgANCypdgsTsBglJB0/ioGtAsAAAzrwQJHZGk8dEn1V1kJxMUp/c2bxWZxAgajhJc2HMaTX+zDxgOZouMQOS0WOqq/GTOUduJEwM9PbBYn8O6XVUVOAAAgAElEQVSWk/j34cvQqlXw9+JkeyJr4RJgVD9FRTeLW0UFoNOJzeMEcooqMHZZEmbGdEJ010DRcYjsnrm1gOfoqH7+/GelffddFrkGqDJKkGQZOq0Ggb4e+OmFu6HV8MAKkTXx/zD6Y5cvA7/+qvRffllsFgdWZZQwfc0hTP36YPUcORY5IuvjiI7+WI8eSvvdd4CaX8zmqDJKeH7NQWxJzYG7Ro3Uy0Xo3aqx6FhELoGFju7s6FEgP1/pP/qo2CwOTKNSoWkjHdy1aiwd04dFjsiGWOjozm6M5m7cpYDMolar8NYj3TE2qjU6B/GCKiJb4nEoqtuWLUqr1QL9+onN4oD0BiP+tvEwTmYXA1CKHYscke2x0FHtZPnmTVXPnBGbxQHpDUZM/fogvkm+iAkr93OBZiKBWOiodp9/rrS9egGtW4vN4oBOZBVj15k86LRqvPN4Dy7QTCQQJ4zT7QwGwM1N6V+7Bvj7i83joH49dQVatQoD2weIjkLkFMytBRzR0e1eeUVpx4xhkTNBRZURPx/Nqv7zoI7NWOSI7AALHdVUWgp88IHSj48Xm8WBVFQZMenLA5i6+iBW7D4nOg4R3YKFjmoaMUJp33qLS32Z4P9+PoHfTufBy12DrsE8vE5kTziPjm7KzQUSEpT+7NlisziYGfd1QOrlQsyM6YyINk1ExyGiW7DQ0U19+ijtunVc6qseyiuNkGQZ3jotGnu7Y8OzUVCpVKJjEdHv8NuMFCdPAhcvKv3YWLFZHEBZpQETVu7H+BX7Uao3AACLHJGdYqEjRefOSrtnj9gcDqDSIGHCyv3Ym56PtKwinMsrFR2JiO6AhY6AX3652Y+KEpfDQbhpVOjdqjEa6bRYNSEC3VvwbutE9owTxl2dLN88H3f2LNC2rdg8DkKWZVwqKEfLxl6ioxC5DE4YJ/OsXKm0nTuzyN1BqV45J3cw4xoA5XwcixyRY2Chc2VGIzBhgtLfvVtsFjtWojfg6RVJ+N+JXExfcwiVBkl0JCIyAQudK/v735V2xAigCed+1eVCfimOZxXD10OLz57qDXct/7chciQ8R+eqyssBL6+bfQ8PsXnsXPKFa3DTqNCjJdf+JBKF5+jINKNHK+3f/84iV4uiiiqs3ncBN34P7NO6MYsckYPiyiiu6OpV4Pvvlf6bbwqNYo+KKqowdlkSUjILUFBWhWmD24uOREQNwBGdK4qMVNovv+RSX7V4f8tJpGQWwN/LDYM6NhMdh4gaiCM6V3P2LHDmjNJ/6imxWezUKw92xuXCCrwY3QHdQjgZnMjRsdC5mvbXD8Pt3AlwbcZqhWVVkGQZjb3d4a3TIn5sX9GRiMhCeNzKldy6juU994jLYWcKyirx5LJEPLVsHwrKKkXHISILY6FzJQMHKu2pU2Jz2BG9wYinlu3DsUtFyC6sQF6JXnQkIrIwFjpXsXat0rZuDXToIDaLHdFpNXgkvAUCGrljzaT+aN/cR3QkIrIwThh3BZIEaDRK/8oVICBAbB47IMtyjfvHXSutRGNvd4GJiOiPcMI41e1f/1LaYcNY5ABcLa3EY5/twa7TedXPscgROS8WOmen1wNz5ij9jRvFZrED+SV6jI5PxKGMAszedIQLNBO5ABY6Zzd+vNLOnAl4eorNYgeullYit1iPZj46rHg6ggs0E7kAnqNzZoWFgP/19RkNhpvn6VzciewiuGnUaNeskegoRGQCnqOj2/3pT0obH+/SRe5KsR6f/u909QLNnYN8WeSIXIhZhW7RokUICwuDh4cHIiMjkZSUVK/X7d69G1qtFj179jTnY8kUGRnAkSNK/5lnxGYRKLe4Ak/EJ2LB1lNYuI3zB4lckcmFbv369YiLi8PcuXNx8OBBhIeHIyYmBrm5uXd8XUFBAcaOHYv77rvP7LBkgrAwpd2+3aWX+lr8y1mcyS1BsJ8HHu/dUnQcIhLA5HN0kZGR6NevHz799FMAgCRJCA0NxfTp0zFr1qw6Xzdq1Ch06NABGo0G33//PVJSUur9mTxHZ6IDB4B+/ZS+/Z+CtaqKKiPm/pCK5wa3Q+um3qLjEFED2OQcXWVlJZKTkxEdHX3zDdRqREdHY+/evXW+bsWKFUhPT8fcuXPr9Tl6vR5FRUU1HmSCG0Xu+HGxOQTJLapAdmEFAMDDTYN3RvRgkSNyYSYVury8PBiNRgQGBtZ4PjAwENnZ2bW+5vTp05g1axa+/vpraLX1u1nC/Pnz4efnV/0IDQ01JaZr27RJaZs3Bzp3FptFgJyiCoxamojR8YnIKaoQHYeI7IBVr7o0Go0YPXo05s2bh44dO9b7dbNnz0ZhYWH1IzMz04opnYgsA489pvSPHRObRQC9wYgn4hORnlcKvUHiZHAiAmDi/egCAgKg0WiQk5NT4/mcnBwEBQXdtn1xcTEOHDiAQ4cO4fnnnwegnNOTZRlarRZbt27FkCFDbnudTqeDTqczJRoBwIIFShsdDTRzvTtj67QaTPlTO3zyy2msmdgfoU28REciIjtg1sUoERER+OSTTwAohatVq1Z4/vnnb7sYRZIkpKWl1Xhu8eLF+N///odvvvkGbdq0gbf3H5874cUo9VBVBbhfX6+xtBTwcp0v+SqjBDfNzYMTFVVGeLi57rxBImdlbi0w+Q7jcXFxGDduHPr27YuIiAh8+OGHKC0txfjrS03Nnj0bly5dwpdffgm1Wo3u3bvXeH3z5s3h4eFx2/PUQFOmKO0LL7hUkbtUUI4xy/Zh5gOd8NBdwQDAIkdENZhc6GJjY3HlyhXMmTMH2dnZ6NmzJxISEqovUMnKykJGRobFg9IdlJQAy5cr/YULxWaxoUsF5Ri1dC8yr5bjva0nEd01sMbIjogI4FqXzqF/f2DfPmDRIuC550SnsZmM/DLELt0Ld60aayf1R4g/F60mcmbm1gIWOkd3+TLQooXSlySXWwXlQn4p3LVqBPuxyBE5Oy7q7Krat1fahASXKHKZV8vwz/+kwWBUpg60burNIkdEd2TyOTqyI0eOAOXlSj8mRmwWG8jIL8MT8Ym4VFAOrUaF2Q91ER2JiBwAR3SOLDxcaY8eFZvDRlbtPY9LBeVoG+CNCQPbiI5DRA6CIzpHtXmz0vr4AC4yVWP2Q53hplFjwsAwNPf1EB2HiBwEL0ZxRLIMqK8PxrOygFpWpXEW5/NKoVKBizITES9GcSnXb5GEAQOcusidyyvFqKWJGLU0ERn5ZaLjEJGDYqFzNEajsvoJAGzdKjaLFVVUGfHUF/uQXVSBRjotPN252gkRmYeFztHcKHKTJwP1WCfUUXm4afDGsC7oGuyLtZP7o5kPF/kmIvPwHJ0jKSu7WdyqqoB63t/PkZTqDfDW3fy5DEYJWi7rRUTgOTrXMGyY0i5c6JRF7kxuMYa8vwMb9t+8/yCLHBE1FL9FHMWVK8COHUr/xReFRrGGM7nFGLV0H3KK9Fi55zyqjLxpKhFZhvMNC5xVl+urgPz4o1Mu9eXj4QZfDy2a++iwemIk70JARBbDc3SO4PhxoGtXpW//f11myymqgLtGjcbe7qKjEJEd4jk6Z3ajyB06JDaHhZ3MLkbc+hToDUYAQKCvB4scEVkcD13au+3blVarBXr2FJvFgk5kF2F0/D5cLa1EgI8Orw3lAs1EZB0c0dm76GilPX9eaAxL+yHlMq6WVuKuFn6Ydm970XGIyIlxRGfPvvhCaXv1unlzVScx84FO8Pd0w6h+reDn5SY6DhE5MV6MYq8kCdBcX/aquBho1EhsHgtIvVwItUqFLsEu8ndIRBbFi1GczauvKu2YMU5R5I5dKsTo+H148ot9OJNbLDoOEbkQFjp7pNcDCxYo/eXLxWaxgIoqIyauOoDC8iq0burFe8kRkU2x0NmjRx9V2rffdoqlvjzcNFjw13BEtW2KLydEwNeD5+SIyHZ4js7eXLsGNGmi9CXJoVdBySvRI6DRzbsOyLIMlQP/PEQkFs/ROYsePZT2m28cusgdzizAkAU7EL8zvfo5FjkiEoGFzp6cPQtcvKj0H39cbJYGOHapEE99sQ9FFQZsTcvmAs1EJJTjnwByJu2vT5xOShKbo4FaNvZEq6Ze8HbXYsX4flygmYiEYqGzF7t23ez36ycuhwX4e7lX34Hg1puoEhGJwF+17cU99yjthQtic5gp+cJVjF+RhFK9AYBS7FjkiMgesNDZg6+/VtpOnYBWrcRmMcOB81cxdlkSfjl5BR9tPy06DhFRDSx0osmysvoJ4LDn5vaezUdppRED2jXFS9EdRcchIqqBx5ZEmzNHaf/6V8BB5wg+P6Q9Av08MLxHCDzdNaLjEBHVwAnjIlVVAe7XbzRaWQm4Oc6KIfvS86HVqNGndWPRUYjIRXDCuCOKjVXauXMdqsjtPZuPp1fsx7jlSTieVSQ6DhHRHbHQiVJYCGzapPTnzhWbxQQVVUbMWHcI5VVG9GndGG0CvEVHIiK6IxY6Ufr2Vdo1axxqqS8PNw2WjOmDh3sE4/MxfeDhxnNyRGTfeI5OhAsXgLAwpW//ux8AcD6vFK2benG9SiIShufoHMmNIrdnj9AY9bXrdB5iPtyJ97achAP8XkREVAMLna3dOlcuKkpcjnpKvnANz6zaD71BwsnsYhglFjoiciycR2drkZFKm55+5+3sROcgH4S39IevpxaLnuwNLRdoJiIHw0JnSxs3Km3r1kCbNmKz/AFJkqFWq+Ctu3kHAnctixwROR5+c9mKLAMjRyr9lBSxWf7ALydz8dhne3CttBIA4K3TssgRkcPit5etzJ+vtA8/DPj7i81yB7+cyMWzXyYjJbMAn+90jMOrRER3wkJnCwYD8PrrSv/bb8Vm+QPn8kpRaZTwUPcgvPwAF2gmIsfHc3S2MG6c0r766s21Le3UhLvbILSJF+7t1Ix3Bicip8AJ49ZWWgo0aqT0JckuV0HZlpYDnVaNP3VsJjoKEVGdOGHcXg0YoLQrVthlkduSmo3nVidj0pcHcOxSoeg4REQWx0JnTZcvA0eOKP2nnxYapTYVVUbM/SEVVUYZD3QLQucgH9GRiIgsjuforKllS6XdsUNojLp4uGmwakIEVu+7gDkPd+VkcCJySvxms5aUlJsLNg8aJDbL7xy5WADp+lJenYJ88I9HurPIEZHT4rebtfTqpbSnTonN8Ts/HcnCo4v34I0fjlUXOyIiZ8ZCZw3//rfSNm8OdOggNsstEtPz8cK6QzBKMioqjWCZIyJXwHN0libLwJ//rPSPHxeb5Xd6tfLHoI7N0NjLHe+O6AGN2v6uAiUisjQWOkv74AOlHTIEaNJEbJbrKg0S3LVq6LQafPZUb2jVahY5InIZPHRpSZIEvPyy0v/5Z7FZrvv+0CU89NFOZBdWAAB0Wg2LHBG5FBY6S5o8WWlfeMEulvradOgi4jak4OyVUqxJyhAdh4hICBY6SykvB5YtU/o3Dl8KVmWUIQN4IiIUL95nPxfFEBHZEs/RWcrgwUq7ZAmgto/fH0b2DUWbAG/0adUYah6uJCIXxUWdLSE3FwgMVPqCd+fGA5nw1mkx9K5goTmIiCzN3FrAEZ0ltGmjtFu3Co2xYX8mXv3uCNQqFVo18UL3Fn5C8xAR2QOzjrEtWrQIYWFh8PDwQGRkJJKSkurc9rvvvsP999+PZs2awdfXF1FRUdiyZYvZge1OWhpQVqb0779fWIyKKiM+/t9pyDLwVGQrdAuxw5EvEZEAJhe69evXIy4uDnPnzsXBgwcRHh6OmJgY5Obm1rr9zp07cf/992Pz5s1ITk7G4MGDMXz4cBw6dKjB4e1Ct25Km5YmNIaHmwZrJvbHi9Ed8Oafu0Flh7cEIiISweRzdJGRkejXrx8+/fRTAIAkSQgNDcX06dMxa9aser1Ht27dEBsbizlz5tRre7s9R7d1KxATo9xYtbhYSITfTl9BVNumXJSZiJyeTW68WllZieTkZERHR998A7Ua0dHR2Lt3b73eQ5IkFBcXo8kdVg3R6/UoKiqq8bBLMTFKm54u5ONX77uAMcuS8OL6FBi5QDMRUa1MKnR5eXkwGo0IvHGF4XWBgYHIzs6u13ssWLAAJSUlGDlyZJ3bzJ8/H35+ftWP0NBQU2LaxmefKe2AAUCzZjb/+F2n8/D6pmMAgEBfD3D2ABFR7Wx6vGvNmjWYN28eNmzYgObNm9e53ezZs1FYWFj9yMzMtGHKepAk4LnnlP727UIiRLVrikd6hmDi3W3wxrAuPCdHRFQHk6YXBAQEQKPRICcnp8bzOTk5CAoKuuNr161bh4kTJ2Ljxo01Dn3WRqfTQafTmRLNtmbMUNrJkwEPD5t+dInegEY6LTRqFRaO7Am1CixyRER3YNKIzt3dHX369MH2W0YxkiRh+/btiIqKqvN1a9euxfjx47F27VoMGzbM/LT2QK8Hrl+IU3340kZW7D6H+xf+ivN5pQAAjVrFIkdE9AdMPnQZFxeH+Ph4rFq1CsePH8fUqVNRWlqK8ePHA1AOO44dO7Z6+zVr1mDs2LF4//33ERkZiezsbGRnZ6OwsNByP4UtPfig0n70kU2X+lqx+xzm/TsNWYUV+PlY/c6HEhGRGSujxMbG4sqVK5gzZw6ys7PRs2dPJCQkVF+gkpWVhYyMmyvlL126FAaDAdOmTcO0adOqnx83bhxWrlzZ8J/AlvLzgR07lP4LL9j0o4P9PKFVq/DsoLaYMqitTT+biMiRca1LUwQEKMXup5+AoUNt/vEns4vRMbARD1cSkUuyyTw6l3bqlFLkAJsVufid6Vh3y33kOgX5sMgREZmIizrXV6dOSnv4sE0+7vNfz2L+zycAAN1b+HGBZiIiM3FEVx83zstptUCPHlb/uPJKI9YfUOYOvhjdgUWOiKgBOKKrjxs3VbXRxHVPdw3WTeqPhNRsjI0Ks8lnEhE5K47o/sjy5UrbqxfwB5PiG+rHw5dRUWUEADT39WCRIyKyABa6O5Fl4JlnlP7u3Vb9qE+2n8YLaw9hytfJXKCZiMiCWOju5NVXlXbsWMDT02ofs+NkLt7fdgoA0C+sCTRcoZmIyGI4j64uVVWAu7vSNxqtugqKLMuY80MqQvw9MfXedlb7HCIiR2ZuLeDFKHV55BGl/b//s0qRk2UZV0r0aO7jAZVKhX88wruCExFZAw9d1qawEPj5Z6V/4/ClBcmyjA+2nULMBztxPEu5qSyLHBGRdbDQ1aZ7d6X97jurvP3H28/g4/+dwbWyKhw4f9Uqn0FERAoeuvy9c+eAixeV/qOPWuUj+oY1hk6rxsyYThjDKQRERFbFQvd7ba/fGeDAAat9xMD2Afjlb/cixN96V3ISEZGChy5vtWfPzX6fPhZ7W1mW8U7CCSz59Wz1cyxyRES2wRHdrQYOVNobhy4tQJZl/F/CCXz+azoA4O72AVy7kojIhjiiu2HNGqXt1Alo0cJib1tRJWH3mTwAwD8e6cYiR0RkY5wwDihLfd2YK1dSAnh7W/TtC8oq8dvpPAwPD7Ho+xIRuRLeeLUh3nxTaUeMsEiRk2UZq/acR4neAADw93JnkSMiEoQjOqNRuc8cABgMgEbToLeTZRn/+E8aVuw+j4iwJlg7uT/XriQisgCO6Mw1cqTSvvlmg4scAOw4eQUrdp8HADzauwWLHBGRYK49oispAXx8lL4kARZYhuvG8l4h/p4YFdGqwe9HREQKLupsjt69lXbt2gYVOVmWcSG/DGEB3lCpVIh7oJOFAhIRUUO57qHLixeB06eV/qhRZr+NJMn4+w/HMOzj35B8getWEhHZG9ctdKGhSrt3b4PeZv7Px/F1YgbKqoy4kF9mgWBERGRJrlnobl3Hsn//Br3VsB4h8PN0w/t/DcdjvVs2MBgREVmaa56j69dPac+fN+vlkiRDpVLuIdcz1B87XxkMP083y+UjIiKLcb0R3bffKm2rVkDr1ia/XJJkzPruCN5JOIkbF6yyyBER2S/XG9GNGKG0R4+a/FJJkvHqt0ewMfki1Crgz+Eh6BpihSXJiIjIYlxrRDd/vtIOGwaYMR+v0ijhwtUyqFXAh6N6scgRETkA15kwLkk3Vz6pqrq57JeJSvUGHMy4hns6NDMvBxERmYVLgP2RsWOVdtYsk4qcUZKxcNspXC2tBAB467QsckREDsQ1RnRlZTfvSmDCUl9GScbLG1Lwfcpl9Gjph03PDeTalUREgnBEdydRUUq7apVJS339dvoKvk+5DK1ahefubcciR0TkgJx/RJeVBYRcvxecGT/qsl3n0MLfEw92DzL5tUREZDlc1LkuLVoo7c6d9drcYJRwIrsY3Vv4AQCeubuNtZIREZENOPehyyNHbo7i7rnnDzc3GCXMWJ+Cxz7bg19PXbFyOCIisgXnLnTh4Up75ky9Np/7Yyp+OpIFWZZRaZCsGIyIiGzFeQvdf/6jtM2bA+3a1esl4waEIcTPA0ue6oP7uwZaMRwREdmK816McuPqymvXAH//OjerMkrQqFRQX7+isqLKCA83TUMjExGRhXF6wa0+/FBphwy5Y5GrNEiYtvogXv/+KCRJqfcsckREzsX5Cp0sAy+9pPQTEurcrNIgYdqag9ialoNvD17CqdxiGwUkIiJbcr5CN3my0s6YAbjd+fY5sizDXatG/Ni+6BzEBZqJiJyRc52j0+sBDw+lX4+lvvQGI05kFSM8tO7Dm0REZB94jg4ABg1S2s8/r7XI6Q1G/P37Y8gqLAcA6LQaFjkiIifnPIXuyhVg3z6lf+Pw5S0qqoyY8lUyvkq8gPEr9sMo2f1AloiILMB5Cl1YmNL+97+1/uekc1ex49QVeLipMefhrlygmYjIRTjHWpfHjyu34gGA++6rdZM/dWyG90aEI8TPAwPaB9gwHBERieQcha5rV6U9frzG0xVVRhy8cK26sI3o09LWyYiISDDHP3S5bZvSNmoEdO5c/XRFlRGTvjyAMcuTsPlolqBwREQkmuMXugceUNpz52o8/cb3x/Db6TzotGoENNIJCEZERPbAsQvdkiVKGxUFBNQ87zbjvg7oGNgIqyZEIKJNEwHhiIjIHjjuhHFZBtTX63RFBaDTobzSCK1GBTeN8rxRknl1JRGRk3C9CeMzZijt5MmAToeySgPGr0zCi+tTYDAq95JjkSMiIse86rKqCvjkE6W/ZAnKK40Yv2I/9p27ikY6Lc7nl6J9cx+xGYmIyC445oguJkZpP/oIUKngrlUj2M8DPjotvnwmgkWOiIiqOd45OqMRaHL94pJbohslGefzS9GuWSNBKYmIyJpc5xxdx44AgJIf/oNpqw/iXF4pAOV8HIscERH9nmMVurNngbw8FLt7YtylxvjpaBamfJVcfXdwIiKi33OsQte7NwDg6H9+xeHMAvh6aPHeX3tAzasriYioDmYVukWLFiEsLAweHh6IjIxEUlLSHbffsWMHevfuDZ1Oh/bt22PlypXmfKxCo8GA+/th8ZO9sXpif/RoyfvJERFR3UwudOvXr0dcXBzmzp2LgwcPIjw8HDExMcjNza11+3PnzmHYsGEYPHgwUlJS8OKLL2LixInYsmWLyWG3tYsALl0CADzQLQh3tfQz+T2IiMi1mHzVZWRkJPr164dPP/0UACBJEkJDQzF9+nTMmjXrtu1fffVV/PTTTzh27Fj1c6NGjUJBQQESEhLq9Zk3rrQJfXEDFjzZHyP7hpoSmYiInIC5V12aNGG8srISycnJmD17dvVzarUa0dHR2Lt3b62v2bt3L6Kjo2s8FxMTgxdffLHOz9Hr9dDr9dV/LiwsBAA0UunR2kf5YYmIyLXc+O43dVacSYUuLy8PRqMRgYGBNZ4PDAzEiRMnan1NdnZ2rdsXFRWhvLwcnp6et71m/vz5mDdv3m3PH/9gDPp/YEpiIiJyNsXFxfDzq/+pK7tcAmz27NmIi4ur/nNBQQFat26NjIwMk344Z1dUVITQ0FBkZmaaNIx3Bdw3teN+qRv3Te3sab/Isozi4mKEhISY9DqTCl1AQAA0Gg1ycnJqPJ+Tk4OgoKBaXxMUFFTr9r6+vrWO5gBAp9NBp7v9HnJ+fn7Cd7Q98vX15X6pA/dN7bhf6sZ9Uzt72S/mDHZMuurS3d0dffr0wfbt26ufkyQJ27dvR1RUVK2viYqKqrE9AGzbtq3O7YmIiCzJ5OkFcXFxiI+Px6pVq3D8+HFMnToVpaWlGD9+PADlsOPYsWOrt58yZQrS09Pxyiuv4MSJE1i8eDE2bNiAl156yXI/BRERUR00b7755pumvKB79+7w9/fHv/71LyxYsAAAsHr1anTq1AkA8PXXX+PChQt4+umnAQCNGzfGwIEDsXjxYvzzn/9Eamoq3n//fYwYMcK0oBoN7r33Xmi1dnlaURjul7px39SO+6Vu3De1c/T94hB3LyAiIjKXY611SUREZCIWOiIicmosdERE5NRY6IiIyKnZTaETeusfO2bKfvnuu+9w//33o1mzZvD19UVUVJRZd4lwFKb+m7lh9+7d0Gq16Nmzp5UTimHqftHr9Xj99dfRunVr6HQ6hIWFYfny5TZKazum7pfVq1cjPDwcXl5eCA4OxoQJE5Cfn2+jtLazc+dODB8+HCEhIVCpVPj+++//8DUO9/0r24F169bJ7u7u8vLly+XU1FR50qRJsr+/v5yTk1Pr9unp6bKXl5ccFxcnp6WlyZ988oms0WjkhIQEGye3LlP3y4wZM+R33nlHTkpKkk+dOiXPnj1bdnNzkw8ePGjj5NZn6r654dq1a3Lbtm3lBx54QA4PD7dRWtsxZ7/8+c9/liMjI+Vt27bJ586dk/fs2SPv2rXLhqmtz9T9smvXLlmtVssfffSRnJ6eLv/2229yt27d5EcffdTGya1v8+bN8uuvvy5/9913MgB506ZNd9zeEb9/7aLQRUREyNOmTav+s9FolENCQuT58+fXuv0rr7wid+vWrcZzsbGxckxMjFVz2pqp+wHs0TgAAARkSURBVKU2Xbt2lefNm2eNeEKZu29iY2PlN954Q547d65TFjpT98vPP/8s+/n5yfn5+baKKISp++W9996T27ZtW+O5jz/+WG7RooVVc4pWn0LniN+/wg9d3rj1z6238jH31j91be+IzNkvvydJEoqLi9GkSRNrxRTC3H2zYsUKpKenY+7cubaIaXPm7Jcff/wRffv2xbvvvosWLVqgY8eO+Nvf/oby8nJbxbY6c/ZLVFQUMjMzsXnzZsiyjJycHGzcuBFDhw61VWy75Yjfv8IL3Z1u/ZOdnV3ra/7o1j/OwJz98nsLFixASUkJRo4caY2Iwpizb06fPo1Zs2bh66+/dtjVHf6IOfslPT0du3btwrFjx7Bp0yZ8+OGH+Oabb/Dcc8/ZIrJNmLNfBg4ciNWrVyM2Nhbu7u4ICgqCv78/Fi1aZIvIds0Rv3+FFzqyjjVr1mDevHnYsGEDmjdvLjqOUEajEaNHj8a8efPQsWNH0XHsiiRJUKlUWL16NSIiIjB06FAsXLgQq1atstsvLVtIS0vDjBkzMGfOHCQnJyMhIQHnz5/HlClTREcjMwj/1dZWt/5xNObslxvWrVuHiRMnYuPGjbcdYnAGpu6b4uJiHDhwAIcOHcLzzz8PQPmCl2UZWq0WW7duxZAhQ2yS3ZrM+TcTHByMFi1a1Lj1SZcuXSDLMi5evIgOHTpYNbMtmLNf5s+fjwEDBmDmzJkAgB49esDb2xv33HMP3nrrLQQHB1s9t71yxO9f4SM63vqndubsFwBYu3Ytxo8fj7Vr12LYsGG2iGpzpu4bX19fHD16FCkpKdWPKVOmoFOnTkhJSUFkZKQt41uNOf9mBg4ciMuXL6OkpKT6uVOnTkGtVqNly5ZWz2wL5uyXsrKy2w5xazQaAMrNP12ZQ37/ir0WRrFu3TpZp9PJK1eulNPS0uTJkyfL/v7+cnZ2tizLsjxr1ix5zJgx1dvfuLx15syZ8vHjx+VFixbZ/eWt5jB1v6xevVrWarXyokWL5KysrOpHQUGBqB/BakzdN7/nrFddmrpfiouL5ZYtW8ojRoyQU1NT5V9//VXu0KGDPHHiRFE/glWYul9WrFgha7VaefHixfLZs2flXbt2yX379pUjIiJE/QhWU1xcLB86dEg+dOiQDEBeuHChfOjQIfnChQuyLDvH969dFDpZluVPPvlEbtWqlezu7i5HRETIiYmJ1f9t3Lhx8qBBg2ps/8svv8g9e/aU3d3d5bZt28orVqywbWAbMWW/DBo0SAZw22PcuHG2D24Dpv6buZWzFjpZNn2/HD9+XI6OjpY9PT3lli1bynFxcXJZWZmNU1ufqfvl448/lrt27Sp7enrKwcHB8pNPPilfvHjRxqmt75dffrnj94YzfP/yNj1EROTUhJ+jIyIisiYWOiIicmosdERE5NRY6IiIyKmx0BERkVNjoSMiIqfGQkdERE6NhY6IiJwaCx0RETk1FjoiInJqLHREROTUWOiIiMip/T92DeuTYOMFwAAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x3294d8050>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(5,5))\n",
    "arrow(0,0,b[1],b[2],head_width=0.05, head_length=0.03,color=\"r\")\n",
    "plot([0,1.1],[0,1.1],\":\")\n",
    "text(b[1]+.03,b[2],\"b\",color=\"r\")\n",
    "text(1.03,1.06,\"a\")\n",
    "axis([0,1.1,0,1.1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAGwCAYAAADfbKDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XlcVOX+B/DPLDAssgnKoigiaiqJK0ja1YoiNW/rFZe0NNfUNK6WtmjeupdfuVSWVlKiXXGtLG+a5qXMRBBFqNwXVFBZRIVhHWDm+f0xOVcSjBln5swMn/frNS+YM+fMfOc4zofnOec8j0wIIUBEROSg5FIXQEREZEkMOiIicmgMOiIicmgMOiIicmgMOiIicmgMOiIicmgMOiIicmgMOiIicmgMOiIicmgMOiIicmgMOiIicmhKqQtoCp1Oh8uXL8PDwwMymUzqcoiISAJCCJSVlSEoKAhyedPbaXYRdJcvX0ZwcLDUZRARkQ3Iy8tD27Ztm7y+XQSdh4cHAP2b8/T0lLgaIiKyJiEEpqzLxL6jebj00bOGTGgqmT1M06NWq+Hl5YXS0lIGHRFRM1RVo8Xan45h2kM9jM4CnoxCREQ26dtfL+NKmQYA4OqswOio9iY9D4OOiIhszoaMXMxYn4XRielQV9fe0XMx6IiIyKacKSrHK1t/AwAM7OQHD9WdnU5iFyejEBFR8xHWugX+8Wg4zl2pwOuPdL3jy8rYoiMiIrPbuXMnBg4cCG9vb/j6+uKRRx7B2bNnb7vNpZIqw+9j+7fHguHdzHLtNIOOiIjMrqKiAvHx8Th06BBSUlIgl8vx+OOPQ6fTNbj+v9PO477Fe/DfY4Vmr4Vdl0REZHZPPvlkvfurV69Gq1atcOzYMYSHh9d7bNPBXLz+zVEAQGbudcR08zdrLWzRERGR2Z0+fRqjRo1CaGgoPD09ERISAgDIzc29Zd1+IS3R2kOFKYNC8VJsF7PXwhYdERGZ3fDhw9G+fXskJiYiKCgIOp0O4eHhqKmpMawjhIBMJkNoqxbYMete+Lo7W2Q8YwYdERGZ1dWrV3Hy5EkkJibi3nvvBQDs27ev3jqf/pyD4/lleOepHlDIZfBrobJYPQw6IiIyKx8fH/j6+mLVqlUIDAxEbm4u5s2bZ3j8059z8Nb24wCAB7u1xsPhgRath8foiIjIrORyOTZu3IjMzEyEh4fjxRdfxOLFiwHouyszL1wHALxwfxhiuwdYvB4O6kxERFZVq9Vh55ECPNIj0KhjcqZmAVt0RERkcf9Ov4C8a5UAACeFHMMjgqw2kTaDjoiILGrFj2fw+tdHMMoMAzSbgkFHREQWc6aoDEu/PwkAiOsbDE8XJ6vXwLMuiYjIYsJae+DduJ64eL0K0+8Lk6QGBh0REZndiQI17grQnzDyaM82ktbCrksiIjKr9/57CkPe/xlbsy5KXQoAtuiIiMiMVu87h/f+exoAUKTWSFyNHlt0RERkNg+HB6BdSze8MvQuTBnUUepyALBFR0REd0gIAa1OQKmQI8jbFd/NuhfuKtuJF7boiIjIZEIILPn+JKYlH0ZNnX5SVVsKOYAtOiIiMpEQAot3ncTKPWcBAHtPXTH7pKnmwBYdERGZRAjgarl+frmFw7vZZMgBbNEREZGJ5HIZEp64G8MjgjCwk5/U5TSKLToiImoyIQSWp5zG6cIyAPqws+WQAxh0RETUREII/HP7cSzbfQqjPz2AMgkGaDaF0UG3d+9eDB8+HEFB+ikWvv766z/dZs+ePejduzdUKhXCwsKwZs0aU2olIiIJnSkqx+dpFwAAsx7oBA8JBmg2hdFBV1FRgYiICKxYsaJJ6587dw7Dhg3Dfffdh+zsbMyePRsTJ07Erl27jC6WiIik08nfA5+M7YN/PX43nu7fXupymuyOZhiXyWTYunUrHnvssUbXefnll7F9+3YcOXLEsGzkyJEoKSnBzp07m/Q6nGGciEgaQgik51xDdEdfqUux3RnG09LSEBMTU29ZbGws0tLSGt1Go9FArVbXuxERkXUJIbBw21GMSkzHmtRzUpdjMosHXUFBAfz9619b4e/vD7Vajaqqqga3SUhIgJeXl+EWHBxs6TKJiOgPPv4pB5+nXYBMBrg52+/VaDZ51uX8+fNRWlpquOXl5UldEhFRszOib1t0DfTE20/2wIh+9tvgsHhEBwQEoLCwsN6ywsJCeHp6wtXVtcFtVCoVVCqVpUsjIqI/0OkEarQ6uDgp4NtChW0zBsBJYZNtoiazePXR0dFISUmpt2z37t2Ijo629EsTEZERdDqBV7/+DeOTDqKqRgsAdh9ygAlBV15ejuzsbGRnZwPQXz6QnZ2N3NxcAPpux3HjxhnWnzp1KnJycvDSSy/hxIkTWLlyJTZv3owXX3zRTG+BiIjulE4n8MrW37AhIw8Hzl1FxvlrUpdkNkYH3aFDh9CrVy/06tULABAfH49evXphwYIFAID8/HxD6AFAhw4dsH37duzevRsRERFYunQpPv30U8TGxprpLRAR0Z2SyQBPVyfIZcCyET0xqHMrqUsymzu6js5aeB0dEZHlCSFw9LIa4W28pC6lQTZ7HR0REdkmrU7gjW1H8evFEgD6QUBsNeTuBIOOiKgZ0uoE5n7xC9bsP49nkw7azQDNpmDQERE1Q+eKy7HzSAEUchnefDTcbgZoNoX9XupOREQmC2vtgbUTIlFcpsGQuwOlLseiGHRERM1EnVaH3ccKDcHWL6SlxBVZB7suiYiagTqtDvGbf8G05MP4IOW01OVYFYOOiKgZWLnnLLb9chlKuQydAzykLseq2HVJRNQMjB8QgtQzxXhuYAc81D1A6nKsikFHROSgarU6VNVq4eniBA8XJ2yc3B8ymUzqsqyOXZdERA6oVqvDCxuy8PSnB1Bapb9GrjmGHMCgIyJyODqdwMz1WfjuSAFO5Jfh2GW11CVJikFHRORg5HIZegR7wVkhxydj+yC6o6/UJUmKx+iIiBzQ84PD8MjdQWjn6yZ1KZJji46IyAHU1Okwc0MW9p8tNixjyOkx6IiI7JymTovnkzPxn18uY3ryYZRr6qQuyaYw6IiI7Nyl61U4dOE6VEo5lo/qhRYqHpW6GfcGEZGdC23VAuuei0JpVS0GhPlJXY7NYYuOiMgOVddqkXzgAoQQAIDwNl4MuUawRUdEZGeqa7WY/O9M7D11BbnXKjF/SFepS7JpbNEREdmZlXvOYu+pK3B1UmBw59ZSl2Pz2KIjIrIzzw/uiJMFaowf0AH9Q5v3xeBNwaAjIrIDVTVaVNTUwa+FCi5OCnwytq/UJdkNdl0SEdm4qhotnlt7ECNXpeNKmUbqcuwOg46IyIbpdALPrT2I/WevIr+kChevV0pdkt1h0BER2TC5XIZhPQLhoVLi8+ci0audj9Ql2R0eoyMiskFCCMP8cWOi2uPh7gHwbaGSuCr7xBYdEZGNqdDUYexnGdh9rNCwjCFnOrboiIhsSIWmDuOTDiLj/DUcz1fjno73wZ1jV94RtuiIiGzItYoaXLhWAQ8XJVY/248hZwbcg0RENiS4pRs2TOqPck0derT1lroch8AWHRGRxMqqa/HhD6eh1ekHaA5t1YIhZ0Zs0RERSUhdXYtnVmcgK7cEV8o0WPRouNQlORy26IiIJPTZz+eQlVsCL1cn/K1vsNTlOCS26IiIJDTj/jAUlVVjTFR7hLfxkroch8SgIyKystKqWpRV16KtjxucFHIkPNFD6pIcGrsuiYisqLSyFmM/O4C4T9KRd43jVloDg46IyEq0OoFnkjLw68VSVNXqp90hy2PQERFZiUIuw3MDO6CVhwrrJ0XhrgBPqUtqFniMjojIwuq0OigV+nbF8IggPNC1Ndyc+fVrLWzRERFZ0LWKGjy2MhVfZl40LGPIWReDjojIQq5V1GB0YjqOXFLjnV0nUKHhMTkpMOiIiCykpk6HqlotWnmokDyxPwdolgj3OhGRhQR4uWDDpP6oqtWiY6sWUpfTbLFFR0RkRsXlGrz57THU1OkAAEHergw5ibFFR0RkJlfKNBidmI7TReWorNEi4Ym7pS6JwBYdEZHZrD+Qi9NF5QjwdMHkv4RKXQ79ji06IiIzmXl/GKrrtIjrG4wQP3epy6HfMeiIiO5Akboa6upahLX2gFwuw8sP3yV1SfQH7LokIjJRoboaI1elY+SqdJwuLJO6HGoEg46IyARancAzqzOQU1wBlVIBFyeF1CVRIxh0REQmUMhlmDfkLoT6uWPj5P4IbukmdUnUCB6jIyIyQlWNFq7O+tbb4C6tMSDMD04KthlsGf91iIia6HJJFYa8vxer950zLGPI2T7+CxERNcHlkiqMXJWO81crkbT/HCo5aardYNARETWBi5MCbs4KtGvpho2ToznVjh3hvxQRURO0dHfG+kn9UV2rRZC3q9TlkBHYoiMiakTetUrEb8o2dFO2dHdmyNkhtuiIiBqQd60SI1el41JJFZQKGd55KkLqkshEbNERETXg21/zcamkCh383BH/YBepy6E7YFLQrVixAiEhIXBxcUFUVBQyMjJuu35ycjIiIiLg5uaGwMBATJgwAVevXjWpYCIia5g6KBSvDL0LGyf3R4CXi9Tl0B0wOug2bdqE+Ph4LFy4EIcPH0ZERARiY2NRVFTU4PqpqakYN24cnnvuORw9ehRbtmxBRkYGJk2adMfFExGZ04WrFfglrwQAIJPJMPkvHeHvyZCzd0YH3bJlyzBp0iSMHz8e3bp1w8cffww3NzesXr26wfXT0tIQEhKCF154AR06dMDAgQMxZcqUP20FEhFZ0/niCsR9ko6nPzuAI5dKpS6HzMiooKupqUFmZiZiYmL+9wRyOWJiYpCWltbgNtHR0cjLy8OOHTsghEBhYSG2bNmCoUOHNvo6Go0GarW63o2IyFK0OoHn1h5EgboaAZ4ubMU5GKOCrri4GFqtFv7+/vWW+/v7o6CgoMFtBgwYgOTkZMTFxcHZ2RkBAQHw9vbGihUrGn2dhIQEeHl5GW7BwcHGlElEZBSFXIa3n+yBXu28sX5Sf7TyUEldEpmRxc+6PHbsGGbNmoUFCxYgMzMTO3fuxPnz5zF16tRGt5k/fz5KS0sNt7y8PEuXSUTN0LWKGgghAAB9Q1riq2n3MOQckFHX0fn5+UGhUKCwsLDe8sLCQgQEBDS4TUJCAu655x7MnTsXANCjRw+4u7vj3nvvxVtvvYXAwMBbtlGpVFCp+GEjIss5U1SOUYnpGB3ZDi8+2BmA/gQUcjxGteicnZ3Rp08fpKSkGJbpdDqkpKQgOjq6wW0qKyuhVNbPU4VCP8XFjb+kiIis6VxxBUYlpuNKmQa7jhZwgGYHZ3TXZXx8PBITE7F27VocP34c06ZNQ0VFBcaPHw9A3+04btw4w/rDhw/Hl19+iY8++gg5OTlITU3FCy+8gMjISAQFBZnvnRARNZFfC2cE+7jirgAPrJ/UnwM0Ozij/3Xj4uJw5coVLFiwAAUFBejZsyd27txpOEElPz8fubm5hvWfffZZlJWV4cMPP8Tf//53eHt74/7778fbb79tvndBRGQEDxcnrJ0QiTqtgI+7s9TlkIXJhB30H6rVanh5eaG0tBSenp5Sl0NEduhkQRkSvjuO9+N6wcvNSepyyASmZgHHuiQih3eyoAyjE9Ox5+QV/HPHManLIStj0BGRw0s7W4yrFTW4u40XXhnaVepyyMp4BJaIHN6zAzrATaVEbLcAdls2Q2zREZFDOnZZjf1nig33R/QNZsg1Uww6InI4Ry+XYvSn6Ziw9iAyL1yXuhySGIOOiByKVicwc30WSiprcVeAJzr5t5C6JJIYg46IHIpCLsOHo3vj/rta4/PnIuHpwu7K5o5BR0QOIe9apWFYwW5Bnlj9bD+GHAFg0BGRA/j1YgmGLf8Zb2w7yjF06RYMOiKya6cKyzDm0wNQV9fh6GU1qmt1UpdENobX0RGRXWvX0g292vmgqqYOSeMj4eqskLoksjEMOiKyS0IIyGQyuDgpsGpsH2h1Au4qfqXRrdh1SUR2J/PCdTy2cj+KyqoBAC5OCoYcNYpBR0R2JfPCNTyzOgO/5JVg2fenpC6H7ACDjojsytmiCpRr6hAd6osFw7tJXQ7ZAbb1iciujOgXDB93ZwwM8+OJJ9QkbNERkc3LOHcNu44WGO4/2M2fIUdNxqAjIpuWnnMVzyZlYHryYRzIuSp1OWSHGHREZLPqtDq88tVvqKzRIrqjLyKCvaUuiewQg46IbJZSIcdnz/bDU33aInFcX7g4sbuSjMegIyKbczxfDa1OP2ZlBz93LPlbBEOOTMagIyKbknqmGI+tSMWcLb8Ywo7oTjDoiMhmHL1ciglrDkJTp0NpVS3qdBygme4cr6MjIpvRxd8DMV39UV2rxcqne0OlZHcl3TkGHRFJrk6rg1Ihh1Ihx3sje0InBEOOzIZdl0QkqT0ni/DQu3uRd60SAOCkkDPkyKwYdEQkmR9PFmHy55nIKa7AJ3vPSl0OOSgGHRFJpqpGC60QiO3ujwWPdJe6HHJQPEZHRJIZencg/Fqo0KudN5wU/LubLIOfLCKyqv8eK8QXmRcN9yM7tGTIkUWxRUdEVrP7WCGeT85EnU4gyNsF93T0k7okagb4ZxQRWUWdVofFu06gVisw7O5ARIa0lLokaibYoiMiq1Aq5Ph8QhRWp57DS7FdoGR3JVkJP2lEZFEZ565BU6cFAAR4ueCVoV0ZcmRV/LQRkcXs+C0foxLTMT35MGrqOG4lSYNBR0QW8evFEszckAWtTsDDxQkKuUzqkqiZ4jE6IrKI8CAvjOjbFppaHRb/LYJBR5Jh0BGRWVXW1MHNWQm5XIZ/PnY3BMCQI0mx65KIzOab7EsYvHgPThWWAQDkchlDjiTHoCMis/gm+xJe3JSNojINthzKk7ocIgN2XRKRWXi7OUOpkOOJXm0wf0hXqcshMmDQEZFZDOrcCt9MH4Au/h6Qs7uSbAi7LonIZF9kXsRn+84Z7ncN9GTIkc1hi46ITLLlUB5e+vJXCAF0DfDAPWEcoJlsE1t0RGS0Oq0Oa9POQwhgbP/2iO7oK3VJRI1ii46IjHZjgObNh/Iw5S+hkMnYXUm2iy06Imqy748WoEJTBwBo6e6MqYM6MuTI5jHoiKhJ1h/IxeR/Z2L8moOortVKXQ5RkzHoiOhPZeeV4JWtvwHQj2GpUvKrg+wHj9ER0Z+KaOuFaYM7oqZOh9eGdWV3JdkVBh0RNepquQa+LVSQyWR4KbYLADDkyO6w/4GIGrR2/3kMXrIHWbnXAegDjiFH9ohBR0S3+Hf6BSzcdhRl1XXYc/KK1OUQ3RF2XRLRLboFesDdWYFx94RgdkwnqcshuiMMOiK6RZ/2LbHrxb+gjbcruyvJ7rHrkogAAJ/+nIMlu05CCAEAaOvjxpAjh8AWHRFh1d6z+NeOEwCAezr6coBmcihs0RE1c7VaHb4/WggAeOGBThygmRwOW3REzZyTQo41EyKx49d8jOgXLHU5RGZnUotuxYoVCAkJgYuLC6KiopCRkXHb9TUaDV599VW0b98eKpUKISEhWL16tUkFE5F5bD6Yh5LKGgBAC5WSIUcOy+gW3aZNmxAfH4+PP/4YUVFReO+99xAbG4uTJ0+idevWDW4zYsQIFBYW4rPPPkNYWBjy8/Oh0+nuuHgiMs2HP5zGku9P4fN0T3wx9R64OCmkLonIYmTixilWTRQVFYV+/frhww8/BADodDoEBwdj5syZmDdv3i3r79y5EyNHjkROTg5atmxpUpFqtRpeXl4oLS2Fp6enSc9BRHqHc6/jiZX7AQBzY7tg+n1hEldE1DSmZoFRXZc1NTXIzMxETEzM/55ALkdMTAzS0tIa3Gbbtm3o27cv3nnnHbRp0wadO3fGnDlzUFVV1ejraDQaqNXqejciMo/e7Xwwb8hdDDlqNozquiwuLoZWq4W/v3+95f7+/jhx4kSD2+Tk5GDfvn1wcXHB1q1bUVxcjOeffx5Xr15FUlJSg9skJCRg0aJFxpRGRH8i92ol2vm6AQCmDuoocTVE1mPxywt0Oh1kMhmSk5MRGRmJoUOHYtmyZVi7dm2jrbr58+ejtLTUcMvLy7N0mUQO7d3dp/Dguz8h9Uyx1KUQWZ1RQefn5weFQoHCwsJ6ywsLCxEQENDgNoGBgWjTpg28vLwMy7p27QohBC5evNjgNiqVCp6envVuRGSaj386i/dTTkNTp8OJgjKpyyGyOqOCztnZGX369EFKSophmU6nQ0pKCqKjoxvcZsCAAbh8+TLKy8sNy06dOgW5XI62bduaWDYRNdWgzq3Q0t0Zrw3riucGdpC6HCKrM7rrMj4+HomJiVi7di2OHz+OadOmoaKiAuPHjweg73YcN26cYf3Ro0fD19cX48ePx7Fjx7B3717MnTsXEyZMgKurq/neCREZCCEMY1Z2DfRESvwgTLw3VOKqiKRh9HV0cXFxuHLlChYsWICCggL07NkTO3fuNJygkp+fj9zcXMP6LVq0wO7duzFz5kz07dsXvr6+GDFiBN566y3zvQsiMhBC4J1dJ1FeXYd/PNodMpkMPu7OUpdFJBmjr6OTAq+jI2oaIQTe3nkSH/90FgCwYVJ/jl1JDsMq19ERkW2r0wmcKNBfd/rG8G4MOSJwUGcih+KkkOPjp/tg76kreKh7w2dCEzU3bNER2TkhBD7+6SyK1NUAABcnBUOO6CYMOiI7JoTAm98ex/99dwKjEtNRXauVuiQim8OgI7JjWXklWJ16DgDw3MBQzkJA1AAeoyOyY73b+SDhibshBDA6qp3U5RDZJAYdkZ0RQuDoZTXC2+iH1RsVyYAjuh12XRLZESEEFnxzFI+tSMXOIwVSl0NkF9iiI7Ij7+4+hX+nX4BMBpRV10pdDpFdYIuOyI480bst2ni7YvFTEfhb32CpyyGyC2zREdk4nU5AJwSUCjlC/Nzx3/hBcHXm2ZVETcUWHZEN0+kEXtn6G/6+5RdodfphaRlyRMZh0BHZKJ1OYP5Xv2HjwTz855fLOJx7XeqSiOwSg47IRumEQEVNHeQy4N24nugX0lLqkojsEo/REdkopUKO9+J6Ylx0CCI7MOSITMUWHZEN0eoE/rXjOHKvVgLQhx1DjujOMOiIbIRWJzBnyy9YtTcHT392AJo6DtBMZA4MOiIb8cvFEmz75TIUchnmDbkLKiXPriQyBx6jI7IRvdv5YPnIXpDLgCF3B0pdDpHDYNARSahOq0PG+Wu4p6MfAGBYDwYckbmx65JIInVaHV7c/AvGfHoAWw7lSV0OkcNi0BFJ5J1dJ/GfXy5DKZfB281Z6nKIHBaDjkgizw3sgC7+HvhoTB882M1f6nKIHBaP0RFZUa1WB50QUCkV8Pd0wfYXBkKp4N+bRJbE/2FEVlKr1WHm+ixMW3fYcI0cQ47I8tiiI7KCWq0OM9Yfxq6jhXBWyHH0shq92/lIXRZRs8CgI7IChUwG3xYqOCvlWDW2D0OOyIpkQgghdRF/Rq1Ww8vLC6WlpfD09JS6HCKT6HQCp4rKcFcAP8NEpjA1C3iAgMhCNHVazNnyC04WlAEA5HIZQ45IAgw6IgvQ1Gkxbd1hfJF5ERPWHOQAzUQSYtARWcCJ/DLsO1MMlVKOt5/swQGaiSTEk1GILCAi2BuJ4/pCKZdhQJif1OUQNWts0RGZSXWtFt/9lm+4P6hzK4YckQ1g0BGZQXWtFpM+P4RpyYeRlHpO6nKI6CYMOiIz+L/vTuDn08Vwc1agWyDPrCSyJTxGR2QGsx7ohKOXSzE39i5EdmgpdTlEdBMGHZGJqmq00AkBd5USPu7O2DwlGjKZTOqyiOgP2HVJZILKmjpMWHMQ45MOokJTBwAMOSIbxaAjMlJNnQ4T1hxEWs5VHMtX41xxhdQlEdFtMOiIjOSkkKF3Ox+0UCmxdkIkwtt4SV0SEd0Gj9ERGUkmk2FubBeMjmqHtj5uUpdDRH+CLTqiJqjQ6I/JHc69DkAfdgw5IvvAoCP6E+WaOjyblIEfThRh5vos1NTppC6JiIzAoCP6ExeuVuB4fhk8XZT46OnecFbyvw2RPeExOqI/0T3IC2snRMJJIUOPtt5Sl0NERuKfpkQNUFfXIvnABQghAAB92vsw5IjsFFt0RH+grq7FuM8ykJ1XgpLKWky/L0zqkojoDrBFR/QHS3edRHZeCbzdnDCocyupyyGiO8QWHdEfvPTwXbhcWo3ZMZ3QPYgXgxPZOwYdEYDSylrohICPuzPcVUokjusrdUlEZCbsuqRmr6SyBmM+S8fTnx1ASWWN1OUQkZkx6KhZ09Rp8fRnB3DkkhoFpdUoLtdIXRIRmRmDjpo1lVKBRyPawK+FM9ZP6o+w1h5Sl0REZsZjdNQsCSEM88dN+ksonurTFj7uzhJXRUSWwBYdNTvXKmrwxEf7se90sWEZQ47IcTHoqFm5Wq7B6MR0ZOWWYP7WXzlAM1EzwKCjZuVaRQ2KyjRo5aFC0rORHKCZqBngMTpqVjr5e2D9pCg4KeTo2KqF1OUQkRXwz1lyeFfKNPjwh9OGAZrvCvBkyBE1IyYF3YoVKxASEgIXFxdERUUhIyOjSdulpqZCqVSiZ8+eprwskdGKyqoxKjEdS74/hWW7T0ldDhFJwOig27RpE+Lj47Fw4UIcPnwYERERiI2NRVFR0W23Kykpwbhx4/DAAw+YXCyRsVb+eBZnisoR6OWCJ3u3lbocIpKATNzoz2miqKgo9OvXDx9++CEAQKfTITg4GDNnzsS8efMa3W7kyJHo1KkTFAoFvv76a2RnZzf5NdVqNby8vFBaWgpPT09jyqVmrrpWi4XfHMXz93VEe193qcshojtgahYY1aKrqalBZmYmYmJi/vcEcjliYmKQlpbW6HZJSUnIycnBwoULm/Q6Go0GarW63o2oqYrU1SgorQYAuDgp8PZTPRhyRM2YUUFXXFwMrVYLf3//esv9/f1RUFDQ4DanT5/GvHnzsG7dOiiVTTvJMyEhAV5eXoZbcHCwMWVSM1aorsbIVekYnZiOQnW11OUQkQ2w6FmXWq0Wo0ePxqJFi9C5c+cmbzd//nyUlpYabnl5eRaskhyFpk6LUYmKL3zQAAAgAElEQVTpyCmugKZOx4vBiQiAkdfR+fn5QaFQoLCwsN7ywsJCBAQE3LJ+WVkZDh06hKysLMyYMQOA/pieEAJKpRLff/897r///lu2U6lUUKlUxpRGBJVSgal/6YgPfjyN9RP7I7ilm9QlEZENMCronJ2d0adPH6SkpOCxxx4DoA+ulJQUQ5DdzNPTE7/99lu9ZStXrsQPP/yAL774Ah06dLiD0on0arU6OCn0nRMj+gXjrz2D4OKkkLgqIrIVRo+MEh8fj2eeeQZ9+/ZFZGQk3nvvPVRUVGD8+PEA9N2Oly5dwueffw65XI7w8PB627du3RouLi63LCcyxaWSKoz97ADmPtQFQ+4OBACGHBHVY3TQxcXF4cqVK1iwYAEKCgrQs2dP7Ny503CCSn5+PnJzc81eKNEfXSqpwshVaci7VoXF359ETDd/Q8uOiOgGo6+jkwKvo6OG5F6tRNyqNDgr5dgwqT+CvF2lLomILMjULOCgzmS32vm6YePk/nBWyhHoxZAjooaxn4fsSt61Srz57THUafWXDrT3dWfIEdFtsUVHdiP3aiVGJabjUkkVlAoZ5g/pKnVJRGQH2KIju7E27TwulVQh1M8dEwbw0hQiahq26MhuzB9yF5wUckwYEILWni5Sl0NEdoJBRzbtfHEFZDL9sTilQo55Q+6SuiQisjPsuiSbda64AiNXpWPkqnTkXq2UuhwislMMOrJJ1bVaPP3pARSoq9FCpYSrM0c7ISLTMOjIJrk4KfDasK7oFuiJDZP7o5UHB/kmItPwGB3ZlApNHdxV+o/lkLsD8WA3fyg5rBcR3QF+g5DNOFNUhvuX7sHmg/+bf5AhR0R3it8iZBPOFJVh5KoDKFRrsGb/edRqOWkqEZkHuy7JJni4OMHTRYnWHiokT4ziLAREZDYMOrIJ/p4u2DC5P5wVcvi4O0tdDhE5EP7ZTJI5WVCG+E3Z0NRpAejDjiFHRObGFh1J4kSBGqMTD+BaRQ38PFR4ZSgHaCYiy2CLjiTxTfZlXKuowd1tvDB9cJjU5RCRA2OLjiQx96Eu8HZ1wsh+7eDl5iR1OUTkwBh0ZDVHL5dCLpOha6An5HIZpgzqKHVJRNQMsOuSrOLIpVKMTjyAMZ8ewJmiMqnLIaJmhEFHFlddq8XEtYdQWlWL9r5unEuOiKyKQUcW5+KkwJK/RSA61BefT4iEpwuPyRGR9fAYHVlMcbkGfi30sw4M7OSHAWG+kMlkEldFRM0NW3RkEb/kleD+JXuQuDfHsIwhR0RSYNCR2R25VIqnPz0AdXUdvj9WwAGaiUhS7Loks2vr44p2vm5wd1YiaXw/DtBMRJJi0JHZebs5G2YguDGJKhGRVPinNplF5oVrGJ+UgQpNHQB92DHkiMgWMOjojh06fw3jPsvAjyev4P2U01KXQ0RUD4OO7lja2auoqNHino6+eDGms9TlEBHVw74lumMz7g+Dv5cLhvcIgquzQupyiIjqYdCRSQ7kXIVSIUef9j6QyWQY0TdY6pKIiBrErksyWtrZq3g26SCeWZ2B4/lqqcshIrotBh0ZpbpWi1kbs1BVq0Wf9j7o4OcudUlERLfFoCOjuDgp8PHYPnikRyA+GdsHLk48JkdEto3H6KhJzhdXoL2vG2QyGXq380Hv0T5Sl0RE1CRs0dGf2ne6GLHv7cXiXSchhJC6HCIiozDo6LYyL1zHc2sPQlOnw8mCMmh1DDoisi/suqTbuivAAxFtveHpqsSKMb2h5ADNRGRnGHTUIJ1OQC6XwV31vxkInJUMOSKyP/zmolv8eLIIT3y0H9cragAA7iolQ46I7Ba/vaieH08UYcrnmcjOK8EnN80OTkRkrxh0VM+54grUaHUYEh6Avz/EAZqJyP7xGB3VM2FgBwS3dMPgLq04MzgROQR+kxF2HyvE3lNXDPcf7ObPkCMih8Fvs2Zu19ECPJ+ciUmfH8KRS6VSl0NEZHYMumasulaLhd8cRa1W4KHuAbgrwEPqkoiIzI7H6JoxFycF1k6IRPKBC1jwSDdeDE5EDonfbM3QrxdLoPt9KK8uAR74x6PhDDkiclj8dmtmtv+aj8dX7sdr3xwxhB0RkSNj0DUj6TlX8cLGLGh1AtU1WjDmiKg54DG6ZqRXO28M6twKPm7OeOepHlDIZVKXRERkcQy6ZqCmTgdnpRwqpQIfPd0bSrmcIUdEzQa7Lh3c11mXMOT9vSgorQYAqJQKhhwRNSsMOge2Nesi4jdn4+yVCqzPyJW6HCIiSTDoHFitVkAAGBUZjNkPdJK6HCIiSfAYnQMb0TcYHfzc0aedD+TsriSiZootOgez5VAedvyWb7jfL6QlQ46ImjW26BzI5oN5ePmrXyGXydCupRvC23hJXRIRkeRMatGtWLECISEhcHFxQVRUFDIyMhpd96uvvsKDDz6IVq1awdPTE9HR0di1a5fJBVPDqmu1WP7DaQgBPB3VDt2DPKUuiYjIJhgddJs2bUJ8fDwWLlyIw4cPIyIiArGxsSgqKmpw/b179+LBBx/Ejh07kJmZifvuuw/Dhw9HVlbWHRdP/+PipMD6if0xO6YT3vhrd8hk7K4kIgIAmRDCqJGgoqKi0K9fP3z44YcAAJ1Oh+DgYMycORPz5s1r0nN0794dcXFxWLBgQZPWV6vV8PLyQmlpKTw92VK52c+nryA61JeDMhORwzM1C4z6dqypqUFmZiZiYmL+9wRyOWJiYpCWltak59DpdCgrK0PLli0bXUej0UCtVte70a2SD1zA2M8yMHtTNrQcoJmIqEFGBV1xcTG0Wi38/f3rLff390dBQUGTnmPJkiUoLy/HiBEjGl0nISEBXl5ehltwcLAxZTYL+04X49WtRwAA/p4u4ImVREQNs2p/1/r167Fo0SJs3rwZrVu3bnS9+fPno7S01HDLy8uzYpX2IbqjLx7tGYSJAzvgtWFdeUyOiKgRRl1e4OfnB4VCgcLCwnrLCwsLERAQcNttN27ciIkTJ2LLli31uj4bolKpoFKpjCmt2SjX1KGFSgmFXIZlI3pCLgNDjojoNoxq0Tk7O6NPnz5ISUkxLNPpdEhJSUF0dHSj223YsAHjx4/Hhg0bMGzYMNOrbeaSUs/hwWU/4XxxBQBAIZcx5IiI/oTRXZfx8fFITEzE2rVrcfz4cUybNg0VFRUYP348AH2347hx4wzrr1+/HuPGjcPSpUsRFRWFgoICFBQUoLS01HzvohlISj2HRf85hvzSanx3pGnHQ4mIyISRUeLi4nDlyhUsWLAABQUF6NmzJ3bu3Gk4QSU/Px+5uf8bKX/VqlWoq6vD9OnTMX36dMPyZ555BmvWrLnzd9BMBHq5QimXYcqgUEwdFCp1OUREdsPo6+ikwOvo9E4WlKGzfwt2VxJRs2SV6+jIuhL35mDjTfPIdQnwYMgRERmJgzrbqE9+OouE704AAMLbeHGAZiIiE7FFZ4OqarTYdEh/7eDsmE4MOSKiO8AWnQ1ydVZg46T+2Hm0AOOiQ6Quh4jIrrFFZ0O2/XIZ1bVaAEBrTxeGHBGRGTDobMQHKafxwoYsTF2XyQGaiYjMiEFnA/acLMLS3acAAP1CWkLBEZqJiMyGx+hswKDOrTC2f3sEebti2uCOUpdDRORQGHQSEULgSrkGrT1cIJPJ8I9HOSs4EZElsOtSAkIIvLv7FGLf3Yvj+fpJZRlyRESWwaCTwPKUM1j+wxlcr6zFofPXpC6HiMihsetSAn1DfKBSyjE3tgvG8hICIiKLYtBJYECYH36cMxhB3q5Sl0JE5PDYdWkFQgi8vfMEPv7prGEZQ46IyDrYorMwIQT+b+cJfPJTDgBgYJgfx64kIrIitugsrLpWh9QzxQCAfzzanSFHRGRlbNFZmKuzAuuei8LPp4sxPCJI6nKIiJodtugsQAiBtfvPo1xTBwDwdnNmyBERSYRB90e1tY0/NngwMHv2bTcXQuAf3x7Dwm1HMSHpIAdoJiKSGIMOACorgS+/BOLiAE9PYNYsk59qz8krSEo9DwB4vHcbDtBMRCSx5nuMrrIS2LED2LwZ+M9/gOpqQKkE6uqArVuB99836WkHd2mFF+4PQ5C3K0ZGtjNz0UREZKzm1aKrqAC2bAH+9jfA11f/c+tWfcgB+pADgEuXGu/CrKsDZswAvLwAPz/g9dchdDqcL64AoB+zMv6hLgw5IiIb4fgtuooKYPt2fctt+3Z9qCkUgFY/k7ch3G6m0wHnzgGdO9/62Nq1wHPPARkZwKFDEJMn45urcrzqF43Pn4tEn/YtLft+iIjIKI7doktN1bfc4uKAr7/+X8vtRsjdzpkzDS8PDgbefRfo0gUYMwZpQ0eh++YkVNZqceFqpflql9CPP/6Ijz76SOoyiIjMwrGDTi7Xt97k8qaF283bnT7d8GP9+wM3TakT8PD9CLl+GUufCMcTvdveYcFNFxISgvfee8/sz1tUVITJkydj3bp1+Prrr83+/ERE1ubYQRcdrQ+sJ5/U35c38e0qFI0HHQCdTkAI/WUDoX7uUMplZg258+fPQyaTGW6+vr546KGHkJWVZbbXaMyMGTOwePFifPnll3jjjTdw7RqnESIi++bYQQcAQUH643O7dgHt2tVrjTWqthY4darBh8SBA5j31a94e+dJfdilp0PWqZM+HM3sv//9L/Lz87Fr1y5UVFRgyJAhKC0tNfvr3Gzz5s147LHHEBAQgOzsbLRsyWOORGTfHD/obnjoIeD4cWDBAv1lBMrfz8Nx1v8QAEoAnASwF8CWw4fxwQcf4LXXXsOkSZMwfPhfcSjzMCpOnUHntxdi99afcOmjJOCDD4BZszB48GDMmDEDM2bMgJeXF/z8/PD6668bWn6m8PX1RUBAAPr27YslS5agsLAQ6enphsfLysowatQouLu7o02bNlixYoXJrzVhwgT06NEDGo0GAFBTU4NevXph3LhxhnVefvlldO7cGW5ubggNDcXrr7+O2ttdYE9EZAuEHSgtLRUARGlpqXme8NQpIQBx+QGVODejhYjwVApnyAT0eWe4yRVK4eLdSrgFdRKuHfuKve4+IqnbIPF5xMOiytVNCB8fIV55RQidTgwaNEi0aNFCzJo1S5w4cUKsW7dOuLm5iVWrVgkhhJgyZYpwd3e/7e2Gc+fOCQAiKyvLsOzw4cMCgNi2bZsQQoj27dsLDw8PkZCQIE6ePCmWL18uFAqF+P777w3bPPzww7d9vW7duhnWLSsrE6GhoWL27NlCCCHmzJkjQkJC6u3zN998U6Smpopz586Jbdu2CX9/f/H222+b59+EiOhPmJoFMiHuoMlhJWq1Gl5eXigtLYWnp+edP2FmJjCwL8qntEALbzlOXZfh0dOPotgpEAp3byjcvSF394Fc5Q7ZH7o6qy8dR+m+ZBz6dh26d+9uWD548GAUFRXh6NGjhm3mzZuHbdu24dixYygqKoJarb5tWWFhYQD0x+g6dOiArKws9OzZEyUlJZgwYQJ2796NM2fOwN/fHyEhIejatSu+++47w/YjR46EWq3Gjh07AACXLl1CVVVVo6/n5OSE9u3bG+6npaVh0KBBmDdvHhISEvDjjz9i4MCBjW6/ZMkSbNy4EYcOHbrt+yIiMgdTs8Dxr6P7I60W6NsXAOAe8zGuZ7yEAnUJyjQ6uPe49883L7+G6vPZCAgIuOWx/v371wvG6OhoLF26FFqtFq1bt0br1q2NKvWee+6BXC5HRUUFQkNDsWnTJvj7+9d7/ptFR0fXOxOzTZs2Rr1edHQ05syZgzfffBMvv/zyLSG3adMmLF++HGfPnkV5eTnq6urM84cHEZEFNZ9jdDcsWKD/2a0bZI+Mgc/fD+Dndi/g0g/rUJq2+U8311WUQKFUwsfHx6iXnTp1Klq0aHHb2x9t2rQJv/zyC65fv46zZ89i6NChRr3mkCFDbvt6N7dIAUCn0yE1NRUKhQJn/nAdYVpaGsaMGYOhQ4fi22+/RVZWFl599VXU1NQYVRMRkbU1rxbd5cvAv/6l/33fPv1Pr7Z4dVECamUqfPDOP7C8+wEkeL6CfPg2+BTaiuvw9WsNeQOXKhw4cKDe/fT0dHTq1AkKhQL/+Mc/MGfOHKPKDQ4ORseOHRt9/OYTU27c79q1q+H+p59++qddlzdbvHgxTpw4gZ9++gmxsbFISkrC+PHjAQD79+9H+/bt8eqrrxrWv3DhglHvh4hICs0r6Dp00P/86CPgDy2yN954A3H4Fl1xGndXvowx8v/DZfjd8hTaihIE3NR9eLPc3FzEx8djypQpOPz7WZtLly4FAJO6Lv9Mamoq3nnnHTz22GPYvXs3tmzZgu3btxseN6brMisrCwsWLMAXX3yBAQMGYNmyZZg1axYGDRqE0NBQdOrUCbm5udi4cSP69euH7du3Y+vWrWZ9P0REltB8ui43bQJudLNNmdLgKl1f/BrX4I0ObpXoffKDBtfRVlxHm6Bbj88BwLhx41BVVYXIyEhMnz4ds2bNwuTJk81SfkP+/ve/49ChQ+jVqxfeeustLFu2DLGxsUY/T3V1NZ5++mk8++yzGD58OABg8uTJuO+++zB27FhotVr89a9/xYsvvogZM2agZ8+e2L9/P15//XVzvyUiIrNrHmddVlYC7u7638+fB2460/BmdVodRrzyPjrJLmD55p+hCuoC74Fj6q1zJXkO4h6MxurVq+stHzx4MHr27GmRYbmIiMj0LGgeLbphw/Q/p09vNOQA4MF5iTgs64yN2sF4/MmnUJq6ASU/J9e76FtXWVLvzEciIrJtjn+MLisL2LNH//ttJlOdOXMmUjduQ+u/vYF34vpi9KDHEOEnx7x58wAIeP3esqstu9bgpQVERGSbHDvodDqgd2/97wcONDgeZa1WhxnTn8eqTz4BAOyZez/COuhbfS+//DLkcjleeuklCCHgFfUEtLU1Dbbo9twIUyIisimOHXTz5ul/RkUBkZG3PFxTp8OAuZ8h57wTABkuXbqIoKCgeuvMnTsXMpkMc+fORV1JAQCwRUdEZEcc9xhdYSGweLH+9x9+uOXhmjodov+eiCuqILQIvx/7jpy9JeRumDNnDpYuXYrK4z8BAI/RERHZEcdt0d1oda1fD7i53fLw2LFjkVfTEa4dWuODuHAM6N7htk8XHx8PpVKJQ4ez0a5dO0tUTEREFuCYQbd+/f9+HzXqloefeuopfPnll4BCiZ+yz+Av4Y2fiXmzF154wVwVEhGRlThe0FVXA2N+v/atqKjeQ5o6LaKfX4Jfv98DALhefAXe3t5WLpCIiKzJ8Y7RDRqk/zl/PtCqlWFxda0WvWd+jGstw9H6qTdw7XoJQ46IqBlwrKA7fBjIyND//s9/1nvovhGTUO7ZDrraavz7xUfh4+0lQYFERGRtjtN1KQTQp4/+9yNHgN/nhRNC4N5770V6aircwy/gmw1r8UB4WwkLJSIia3KcoJs1S/8zJgb4fZ61qpo63DN8DLJTUwEAVzK+haurq1QVEhGRBByj67KwEPjg99kGduwAoA+5u59fgWsRY+DWZQCqqqoYckREzZBjBN2Na+b+8x/AyQlCCIQ/+xbq/MIg6mrw/bYv4eLiIm2NREQkCfvvuvz8c/1PV1fgkUcghEBYWBhyr1ag1ZOt8eX8pxAdZt4JT4mIyH7Y93x0Gg1wo6VWWopKlRs6h4Xi0sU8AEC1pgYqZycJKiYiInNrnvPR3ZiZICEBFSpXdJn8HjS9RwEyOerq6hhyRERkx12XGRnAsWMAgMoX56DzxHfh1KYbXFu2x8nL16FoYEoeIiJqfuyzRSeEfuodAOLMGbT09oSmpBA6TQW+nv0AOgc0vUlLRESOzT5bdFOmAADEo49CHhamX/bflUj/bCHCWntIWBgREdka+wu6/HwgMRFlzq4IxV1Q+hyEm7YcJSUlkP0+GgoREdEN9hd0QUFQO7ui898Wwb1tN7j6d8C5lZMYckRE1CD7OkaXlAQBoHVAGFSBnSFqKvHt63FQKOzrbRARkfWYlBArVqxASEgIXFxcEBUVhYwbMwY0Ys+ePejduzdUKhXCwsKwZs0aU14WutmzIQegyf0NThlr8W38Q+jRllPtEBFR44wOuk2bNiE+Ph4LFy7E4cOHERERgdjYWBT9YZLTG86dO4dhw4bhvvvuQ3Z2NmbPno2JEydi165dRhcb1DESABAeHo4zP32Fu9tyqh0iIro9o0dGiYqKQr9+/fDhhx8CAHQ6HYKDgzFz5kzMmzfvlvVffvllbN++HUeOHDEsGzlyJEpKSrBz584mveaNq+GDZ2+G55ldOPKfT40pmYiIHICpI6MYdTJKTU0NMjMzMX/+fMMyuVyOmJgYpKWlNbhNWloaYmJi6i2LjY3F7NmzG30djUYDjUZjuF9aWgoAEBVX8dmShVCr1caUTUREDuDGd7+xI1caFXTFxcXQarXw9/evt9zf3x8nTpxocJuCgoIG11er1Y1OnZOQkIBFixbdsvxi4jT0TzSmYiIicjRlZWXw8mr6oSubvLxg/vz5iI+PN9wvKSlB+/btkZuba9Sbc3RqtRrBwcHIy8szqhnfHHDfNIz7pXHcNw2zpf0ihEBZWRmCgoKM2s6ooPPz84NCoUBhYWG95YWFhQi4MSfcHwQEBDS4vqenZ6MToapUKqhUqluWe3l5Sb6jbZGnpyf3SyO4bxrG/dI47puG2cp+MaWxY9RZl87OzujTpw9SUlIMy3Q6HVJSUhAdHd3gNtHR0fXWB4Ddu3c3uj4REZE5GX15QXx8PBITE7F27VocP34c06ZNQ0VFBcaPHw9A3+04btw4w/pTp05FTk4OXnrpJZw4cQIrV67E5s2b8eKLL5rvXRARETVC8cYbb7xhzAbh4eHw9vbGP//5TyxZsgQAkJycjC5dugAA1q1bhwsXLuDZZ58FAPj4+GDAgAFYuXIl3nzzTRw9ehRLly7FU089ZVyhCgUGDx4MpdImDytKhvulcdw3DeN+aRz3TcPsfb/YxQzjREREpuIgkURE5NAYdERE5NAYdERE5NAYdERE5NBsJuikmvrH1hmzX7766is8+OCDaNWqFTw9PREdHW3SLBH2wtjPzA2pqalQKpXo2bOnhSuUhrH7RaPR4NVXX0X79u2hUqkQEhKC1atXW6la6zF2vyQnJyMiIgJubm4IDAzEhAkTcPXqVStVaz179+7F8OHDERQUBJlMhq+//vpPt7G7719hAzZu3CicnZ3F6tWrxdGjR8WkSZOEt7e3KCwsbHD9nJwc4ebmJuLj48WxY8fEBx98IBQKhdi5c6eVK7csY/fLrFmzxNtvvy0yMjLEqVOnxPz584WTk5M4fPiwlSu3PGP3zQ3Xr18XoaGh4qGHHhIRERFWqtZ6TNkvf/3rX0VUVJTYvXu3OHfunNi/f7/Yt2+fFau2PGP3y759+4RcLhfvv/++yMnJET///LPo3r27ePzxx61cueXt2LFDvPrqq+Krr74SAMTWrVtvu749fv/aRNBFRkaK6dOnG+5rtVoRFBQkEhISGlz/pZdeEt27d6+3LC4uTsTGxlq0Tmszdr80pFu3bmLRokWWKE9Spu6buLg48dprr4mFCxc6ZNAZu1++++474eXlJa5evWqtEiVh7H5ZvHixCA0Nrbds+fLlok2bNhatU2pNCTp7/P6VvOvyxtQ/N0/lY+rUP42tb49M2S9/pNPpUFZWhpYtW1qqTEmYum+SkpKQk5ODhQsXWqNMqzNlv2zbtg19+/bFO++8gzZt2qBz586YM2cOqqqqrFW2xZmyX6Kjo5GXl4cdO3ZACIHCwkJs2bIFQ4cOtVbZNssev38lD7rbTf1TUFDQ4DZ/NvWPIzBlv/zRkiVLUF5ejhEjRliiRMmYsm9Onz6NefPmYd26dXY7usOfMWW/5OTkYN++fThy5Ai2bt2K9957D1988QWef/55a5RsFabslwEDBiA5ORlxcXFwdnZGQEAAvL29sWLFCmuUbNPs8ftX8qAjy1i/fj0WLVqEzZs3o3Xr1lKXIymtVovRo0dj0aJF6Ny5s9Tl2BSdTgeZTIbk5GRERkZi6NChWLZsGdauXWuzX1rWcOzYMcyaNQsLFixAZmYmdu7cifPnz2Pq1KlSl0YmkPxPW2tN/WNvTNkvN2zcuBETJ07Eli1bbulicATG7puysjIcOnQIWVlZmDFjBgD9F7wQAkqlEt9//z3uv/9+q9RuSaZ8ZgIDA9GmTZt6U5907doVQghcvHgRnTp1smjN1mDKfklISMA999yDuXPnAgB69OgBd3d33HvvvXjrrbcQGBho8bptlT1+/0reouPUPw0zZb8AwIYNGzB+/Hhs2LABw4YNs0apVmfsvvH09MRvv/2G7Oxsw23q1Kno0qULsrOzERUVZc3yLcaUz8yAAQNw+fJllJeXG5adOnUKcrkcbdu2tXjN1mDKfqmsrLyli1uhUADQT/7ZnNnl96+058Lobdy4UahUKrFmzRpx7NgxMXnyZOHt7S0KCgqEEELMmzdPjB071rD+jdNb586dK44fPy5WrFhh86e3msLY/ZKcnCyUSqVYsWKFyM/PN9xKSkqkegsWY+y++SNHPevS2P1SVlYm2rZtK5566ilx9OhR8dNPP4lOnTqJiRMnSvUWLMLY/ZKUlCSUSqVYuXKlOHv2rNi3b5/o27eviIyMlOotWExZWZnIysoSWVlZAoBYtmyZyMrKEhcuXBBCOMb3r00EnRBCfPDBB6Jdu3bC2dlZREZGivT0dMNjzzzzjBg0aFC99X/88UfRs2dP4ezsLEJDQ0VSUpJ1C7YSY/bLoEGDBIBbbs8884z1C7cCYz8zN3PUoBPC+P1y/PhxERMTI1xdXUXbtm1FfHy8qKystHLVlmfsflm+fLno1q2bcHV1FYGBgWLMmDHi4sWLVq7a8n788cfbfm84wvcvp+khIiKHJvkxOiIiIkti0BERkUNj0BERkVoCjMcAAAA9SURBVENj0BERkUNj0BERkUNj0BERkUNj0BERkUNj0BERkUNj0BERkUNj0BERkUNj0BERkUNj0BERkUP7f0SyB61T7SDFAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x329944310>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x329310c50>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = rand(2)      # random red vector\n",
    "\n",
    "\n",
    "figure(figsize=(5,5))\n",
    "arrow(0,0,b[1],b[2],head_width=0.05, head_length=0.03,color=\"r\")\n",
    "text(b[1]+.03,b[2],\"b\",color=\"r\")\n",
    "plot([0,1.1],[0,1.1],\":\")\n",
    "axis([0,1.1,0,1.1]);\n",
    "\n",
    "a = ones(2)    # target direction\n",
    "x̂ = (a'b)/(a'a)\n",
    "p = a * x̂      # projection\n",
    "plot([b[1],p[1]],[b[2],p[2]],\":\")\n",
    "arrow(0,0,p[1],p[2],head_width=0.05, head_length=0.03)\n",
    "text(p[1]+.03,p[2],\"p=Pb=x̂a\")\n",
    "text(1.03,1.06,\"a\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us break this into steps\n",
    "\n",
    "1. Find $\\hat{x}$ \n",
    "2. Find $p$\n",
    "3. Find matrix $P$ such that $Pb=p$\n",
    "\n",
    "To do this form the \"error\" vector  $e = b - p = b - \\hat{x}a$ where $\\hat{x}$ is the unknown.  We choose $\\hat{x}$ specifically to make $e \\perp a$.\n",
    "\n",
    "We want $a \\cdot (b-\\hat{x}a) = 0$ (where $\\cdot$ denotes the dot product) so that \n",
    "\n",
    "1. $\\hat{x} = (a \\cdot b)/(a \\cdot a) = (a^Tb)/(a^Ta) $  we then have\n",
    "2. $p = \\hat{x}a = a\\hat{x} = a (a^Tb)/(a^Ta)$\n",
    "3. $Pb = a (a^Tb)/(a^Ta)$ gives $P = (aa^T)/(a^Ta)$\n",
    "\n",
    "Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.5  0.5\n",
       " 0.5  0.5"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = (a*a')/(a'a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [],
      "text/plain": [
       "Interact.Slider{Int64}(Signal{Int64}(2, nactions=1),\"\",2,1:15,\"horizontal\",true,\"d\",true)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Rational{Int64},2}:\n",
       " 1//2  1//2\n",
       " 1//2  1//2"
      ]
     },
     "execution_count": 6,
     "metadata": {
      "comm_id": "f3733903-cde7-4d0b-ba65-93e75acac018",
      "reactive": true
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@manipulate for n=slider(1:15,value=2)\n",
    "    a = ones(Rational,n)\n",
    "    P = (a*a')/(a'a)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the special case of a being the ones vector, $\\hat{x}$ is the mean of $b$.  If only one number is used to\n",
    "summarize a large data vector $b$, it is commonly the mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider more general $a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAGwCAYAAAApJSV7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtcVHX+x/HXzHATFVBU8ELipUzL1DSJtF0rzc1qq21Xu2yZlmVZ6ZKV/izLarO7VlpulmllZXbfNMvobqZ5K8tbanlFFJWLgODMnN8fXwFZkRxlOHDm/Xw85sGH45yZzwwyH77f8z2f47Isy0JEREQq5LY7ARERkZpMhVJERKQSKpQiIiKVUKEUERGphAqliIhIJVQoRUREKqFCKSIiUgkVShERkUqoUIqIiFRChVJERKQSKpQiIiKVCLM7gaPh9/vZvn079evXx+Vy2Z2OiIjYwLIs8vLyaNasGW539Y3zakWh3L59O0lJSXanISIiNcCWLVto0aJFtT1frSiU9evXB8ybExMTY3M2IiJih9zcXJKSkkprQnWpFYWyZLo1JiZGhVJEJMRV9yE4LeYRERGphAqliIhIJVQoRUREKqFCKSIiUgkVShERkUqoUIqIiFRChVJERKQSKpQiIiKVUKEUERGphAqliIhIJVQoRUREKqFCKSIiUgkVShERkUqoUIqIiFRChVJERKQSKpQiIiKVUKEUERGphAqliIhIJVQoRUREKqFCKSIiUgkVShERkUqoUIqIiFRChVJERKQSKpQiIiKVUKEUERGphAqliIhIJVQoRUREKqFCKSIiUomAC+XXX3/NxRdfTLNmzXC5XLz//vt/uM+XX37J6aefTmRkJG3btmX69OnHkquIiEi1C7hQ5ufn06lTJyZPnnxU9//tt9+48MILOeecc1ixYgUjRozghhtu4JNPPgk4WRERkeoWFugOF1xwARdccMFR33/KlCm0atWKJ598EoD27dvz7bffMmHCBPr27Rvo04uIiFSroB+jXLhwIb179y63rW/fvixcuPCI+xQVFZGbm1vuJiIiYoegF8odO3aQkJBQbltCQgK5ubkUFhZWuM/48eOJjY0tvSUlJQU7TRERkQrVyFWvo0ePJicnp/S2ZcsWu1MSEZEQFfAxykAlJiaSmZlZbltmZiYxMTHUqVOnwn0iIyOJjIwMdmoiIiJ/KOgjytTUVNLT08ttmz9/PqmpqcF+ahERkeMWcKHct28fK1asYMWKFYA5/WPFihVs3rwZMNOm1157ben9hw4dysaNG7nrrrtYs2YNzz33HG+99Rb/+te/qugliIiIBE/AhXLJkiV06dKFLl26AJCWlkaXLl0YO3YsABkZGaVFE6BVq1bMmTOH+fPn06lTJ5588klefPFFnRoiIiK1gsuyLMvuJP5Ibm4usbGx5OTkEBMTY3c6IiJiA7tqQY1c9SoiIlJTqFCKiIhUQoVSRESkEiqUIiIilVChFBERqYQKpYiISCVUKEVERCqhQikiIlIJFUoREZFKqFCKiIhUQoVSRESkEiqUIiLVrVcvGDHC7izkKKlQioiIVCLM7gRERBzD74eMDNi4ETZsMF//+U846SS7M5PjoEIpIhKI/fvht9/KF8P162HdOti0CYqLy9//hx/g448PfxyvF269FV59FcLD4eab4YEHwOWqntchR02FUkTkaFx7LXzyCezcWbbN7QaPBw4cqHgflwvq1q3432bMgOuvh8WLYckSuPFGOOEEGDKk6nOX46JCKSLyR/x++PLL8kWyZLvff+T93O4jF8qkJJgwwRTTdu1g5UrzvQpljaPFPCIif8TthoUL4eSTzQgykP3q1av43848s/w0a2oq/Por+HzHl6tUORVKEZGj0bw5fPcddO9uCuDROlKhlFpDhVJE5Gg1aADp6dCv39EvujlSoVy0qPz3338PJ54Y2IhVqoUKpYhIIOrUgffeg+uu++P7+v1HLpSbN0NaGqxdC2+8Ac8+C8OHV2mqUjW0mEdEJFBhYfDSS1C3Drz2vNmWbR1+P5/vyIXy2muhsNBM5Xo8pkjeeGPwcpZjpkIpInIs8vJMkRxeH4otGJ9X8f0qKpRfflkWP/98UNKTqqOpVxGRQO3fD7GxUDKIjIwwo8KKjltqMU+tp0IpIhIIrxcSEkyc2AzqNIR6jWHOHIiMPHxFrAplraepVxGRo2VZcPrpkJtrvv9pI0RElP37V19B375mWrbkfEgVylpPI0oRkaN1ySWmgw7Avn3liySYhTnff29GnCXTsCqUtZ4KpYjI0bjtNvjvf028e7dpTWdZ4DsA3qKy+7VrZ/q3tmtnmp3HxtqTr1QZTb2KiPyRRx+FSZNMvHUrNGxo4vwseKKtie/LLhtFNm9uWt4VF0OTJtWfr1QpFUoRkcq88gqMGmXidetMESzhOmRSzrLKr3qNi6ue/CToNPUqInIkH38MAweaeMkS02LuUOVOB6mg4YA4gkaUIiIVWbTI9HQFmD8funY9/D51GsCYHYCr/OhSHEWFUkTkf61ZYy6DBaYPa+/eFd/P5YLwOtWXl9hCfwKJiBxqyxZo397EEyfCFVfYm4/YToVSRKREVhaccIKJR47846t5FBfA423hsTZQnB/8/MQWmnoVEQHTTadxYxP/4x/w+ON/vI/LBfm7TGz5g5eb2EojShGRoiKIiTFx9+4wa9bR7Vfu9BAVSqfSiFJEQpvPB0lJJk5IgO++q/gqIBXxRMCtS839I9SqzqlUKEUkdFmWWd266+D06ebN5nJZR8vlgkZtg5Ob1BiaehWR0NW/v2kkAOYY5f82ORdBhVJEQtUdd8Dbb5s4K+vYr/LxQi+YcjYU7Kmy1KRm0dSriISeCRPgqadMvHUrxMcf+2Nl/GgW8vgOVE1uUuNoRCkioeWNNyAtzcRr15Zvcn4sSla+atWrY2lEKSKh47PP4KqrTPzDD3DSScf/mP981yzqiW54/I8lNZIKpYiEhmXLoE8fE3/yCXTrVjWP2/rPVfM4UmNp6lVEnG/9+rKrf7z2Gpx/vr35SK2iQikizrZ9e9l1JB9/HK6+umof/+3B8MaVkJdZtY8rNYYKpYg41969ZYt1hg83jc6r2rpPYe1cKN5X9Y8tNYIKpYg4U0EBNDy4wOayy8wls4JBF2x2PC3mERHnOXAA6tY1cZcu8M47wXuui54CvxfqNg7ec4itVChFxFn8fmjVysQNGpjTQI62yfmx6Pj34D221AiaMxAR57As+NOfYNs2831GRmBNzkUqoEIpIs5xzTWwYIGJ8/IgMjL4z5n+AMy5A3K3B/+5xBYqlCLiDGPGwMyZJj6eJueBWj4TfngRCnZXz/NJtVOhFJHa77nn4OGHTbxly/E1OQ+Uer06nhbziEjt9s47MGyYiVevhhYtqvf5e44w51DWS6ze55Vqo0IpIrXXV1/B3w+uOv3+ezj55OrPIeWm6n9OqVaaehWR2unHH6FXLxPPnQspKbamI86lEaWI1D6//QadO5t4+nS44AL7clnyMhTuhU5XQEwz+/KQoDmmEeXkyZNJTk4mKiqKlJQUFi9eXOn9Z86cSadOnYiOjqZp06YMHjyY3bu1QkxEjkFmJrRubeLx42HgQHvzWfA0pI+D7C325iFBE3ChnDVrFmlpadx3330sW7aMTp060bdvX3bu3Fnh/RcsWMC1117L9ddfzy+//MLs2bNZvHgxQ4YMOe7kRSTE5ORA4sFFM8OGwahR9uYDWvUaAgIulE899RRDhgxh0KBBdOjQgSlTphAdHc20adMqvP/ChQtJTk7m9ttvp1WrVvTs2ZObbrrpD0ehIiLlFBZCXJyJL7wQJk2yN58SHf8O3a6Hek3szkSCJKBCWVxczNKlS+ndu3fZA7jd9O7dm4ULF1a4T2pqKlu2bGHu3LlYlkVmZiazZ8+mX79+R3yeoqIicnNzy91EJIR5vRAdbeIOHeC//7U3n0Od83+mMXp8G7szkSAJqFBmZWXh8/lISEgotz0hIYEdO3ZUuE+PHj2YOXMmAwYMICIigsTEROLi4pg8efIRn2f8+PHExsaW3pKSkgJJU0ScxLKgXTsTR0fDTz8Ft8m5yP8I+ukhq1atYvjw4YwdO5alS5cyb948fv/9d4YOHXrEfUaPHk1OTk7pbcsWHSQXCUmWBX36wMaN5vs9e2pek/N1n8KPsyCv4sGC1H4BnR7SqFEjPB4PmZmZ5bZnZmaSmFhxV4rx48dz1llnceeddwJw2mmnUbduXc4++2weeughmjZtetg+kZGRRFZHM2MRqdmGDIH0dBPn5lZPk/NApT8AmSvhmvegvrrzOFFAI8qIiAi6du1Kesl/XMDv95Oenk5qamqF+xQUFBAWVr4eew7+RWhZVqD5ikioGDcOXnrJxLt2Qf369uZzJCXTwFr16lgBNxxIS0tj4MCBdOvWje7duzNx4kTy8/MZNGgQYKZNt23bxiuvvALAxRdfzJAhQ3j++efp27cvGRkZjBgxgu7du9OsmU7OFZEKTJ0K999v4s2boVEjW9OpVMseENMcoquxEbtUq4AL5YABA9i1axdjx45lx44ddO7cmXnz5pUu8MnIyGDz5s2l97/uuuvIy8tj0qRJ3HHHHcTFxXHuuefy6KOPVt2rEBHn+PBDuPFGE69aBTV9Md8Fj9idgQSZy6oF85+5ubnExsaSk5NDTEyM3emISLB89x306FEWH+GQjoQmu2qBer2KSM2walVZkfzvf2tPkcz4CfbnQJP2ULcGTxHLMdPVQ0TEflu2wCmnmPill+Cii+zNJxBzR8KMi2BzxU1XpPZToRQRe2VlwQknmPjBB2HwYHvzCZR6vTqepl5FxD55edC4sYlvvBHuucfefI5FwzZQvA8itX7CqVQoRcQeRUVQsiDj/PPhP/+xN59jdemR23GKM2jqVUSqn89XdiWQE0+EefPszUekEhpRikj1sizo2BH274ewMLPatTY3Oc/LhAP5ULcxRNbQ7kFyXDSiFJHqddFFsHq1ifPyTLGszd67CZ7pAmvm2p2JBIkKpYhUn1tugbkHC0puLkRF2ZtPVdCqV8dToRSR6vHII/D88yauyU3OAxUVY/q8esLtzkSCpJbPeYhIrTBjBowebeJNm2p2k/NA/WO63RlIkGlEKSLB9fHHcN11Jv7557LmAiK1hEaUIhI8ixdDv34m/uabsjZ1TuLzguUDlwc8+kh1Io0oRSQ41q6FlBQTv/8+9Oxpbz7BMnsgPNQEls2wOxMJEhVKEal627fDySeb+D//gUsusTefYCo9B7TGX7FQjpEKpYhUrb17oXlzE48dW3YRZqcqPT1EhdKpdOFmEak6+flQr56Jr7sOXn7Z1nSqhbcYsMAdBm6P3dk4mi7cLCK1W3FxWZH8859Do0gChEXYnYEEmaZeReT4+f1ll8s64QT44gt78xGpQiqUInJ8LAu6djUt6QA2bKjdTc4DNWckPNYGFk+1OxMJEhVKETk+l18OK1aYuLCw9jc5D1TxPijIggOFdmciQaJCKSLHLi0N3nvPxDk5zmhyHrCDo2c1RXesEPvTT0SqzFNPwYQJJt65E0J1RXrv++FPIyG6od2ZSJCoUIpI4N54A+64w8S//Va2kCcU1U8AEuzOQoJIU68iEpjPPoOrrjLxTz9BcrKt6YgEmwqliBy95cuhTx8Tf/UVdOxobz41wTdPwpSzYck0uzORIFGhFJGjs2EDnH66id9+G/70J3vzqSlytsGOnyAv0+5MJEhUKEXkj2VmQtu2Jp482ZwSIkZJr1c1RXcsLeYRkcrl5EBioolHj4ZbbrE3n5omZSh0uATikuzORIJEhVJEjqywEOLiTHz11fDww/bmUxM1amtu4liaehWRinm9EB1t4rPOgldftTcfEZuoUIrI4fz+smtKJibCt9+GVv/WQCyfCa9fYb6KI2nqVUTKsywzgty503y/ZYuKZGV2rYZ1H2v61cE0ohSR8q6+GhYtMnEoNjkPVMmqV0urXp1KvwEiUmbUKNOeDkK4yXmATrkMGreHxu3szkSCRIVSRIxJk+DRR02cmRm6Tc4D1ayLuYljaepVROCdd+C220y8cSM0aWJvPiI1iAqlSKj76iv4+99NvGIFtGplbz61zbpPYc5I+PkduzORIFGhFAllK1dCr14m/vxz6NTJ1nRqpe3L4Iep8PsCuzORIFGhFAlVmzbBaaeZeNYsOOcce/OprUpXvfrtzUOCRot5RELRrl1l15F8+mno39/WdGq15LPh3HuhqUbjTqVCKRJq8vLKFuvceSfcfru9+dR2LVPNTRxLU68ioaSoqOy0j3/8Ax57zN58RGoBjShFQoXPV9ZAoFs3c1xSjt/WJfDb19CkA7T7i93ZSBBoRCkSCiyr7Jhkw4amRZ36t1aNTd9B+jhY9b7dmUiQqFCKhIJzzoGtW02cmQlu/epXGa16dTxNvYo43aBBpqkAQEGBmpxXtcRTodv10Px0uzORINFvjIiT3XcfTJ9u4uxsqFPH1nQcqXUvcxPH0vyLiFO98AI88ICJd+yA2Fh78xGppTSiFHGiDz+Em24y8YYNkJBgbz5OlrXetLGLaQ7JPezORoJAI0oRp/nuO7jkEhMvWwatW9ubj9P99iW8OwQWTbE7EwkSFUoRJ1m9GnocHNXMnw9ddJ3EoNOqV8dToRRxiq1boUMHE8+cCb1725tPqIg9AU76CzTrbHcmEiQ6RiniBHv2QFKSiZ94Aq66yt58QsmJvc1NHEsjSpHaLj8f4uNNPGIE3HGHvfmIOIxGlCK1WXEx1Ktn4ksvhQkT7M0nFO3bCVnrICoWEjvanY0EgUaUIrWV31/W5Py00+Ddd+3NJ1Rt+BymXwjz77M7EwkSFUqR2siyoF0787VuXVi+XE3ObXPwfdeqV8dSoRSpjf7yF1i/3sR796rJuZ3qNDBTrg1b2Z2JBMkx/XZNnjyZ5ORkoqKiSElJYfHixZXev6ioiDFjxtCyZUsiIyNJTk5m2rRpx5SwSMgbOhQ+/dTEBQUQHm5vPqHupPNh6LdwkY4PO1XAi3lmzZpFWloaU6ZMISUlhYkTJ9K3b1/Wrl1LkyZNKtynf//+ZGZm8tJLL9G2bVsyMjLw+zVNIRKwhx+G//zHxGpyLlItXJZlWYHskJKSwhlnnMGkSZMA8Pv9JCUlcdtttzFq1KjD7j9v3jyuuOIKNm7cSMOGDY8pydzcXGJjY8nJySEmJuaYHkOk1ps+3VwyCyAjAxITbU1HDirOh32ZEBYFMc3szsbR7KoFAU29FhcXs3TpUnof0vHD7XbTu3dvFi5cWOE+H374Id26deOxxx6jefPmnHTSSYwcOZLCwsIjPk9RURG5ubnlbiIh7eOPy4rkr7+qSNYk69PhmS7w9mC7M5EgCWjqNSsrC5/PR8L/XIkgISGBNWvWVLjPxo0b+fbbb4mKiuK9994jKyuLW265hd27d/Pyyy9XuM/48eMZN25cIKmJONfixdCvn4l/+AHatrU3HylPvV4dL+hL5fx+Py6Xi5kzZ9K9e3f69evHU089xYwZM444qhw9ejQ5OTmlty1btgQ7TZGaad06SEkx8bx50K2bvfnI4cIioU5DiKxvdyYSJAGNKBs1aoTH4yEzM7Pc9szMTBKPMBXUtGlTmjdvTuwhF41t3749lmWxdetWTjzxxMP2iYyMJDIyMpDURJwnI8OcKwnwyivQt6+9+UjFTuwDd/9mdxYSRAGNKCMiIujatSvp6eml2/x+P+np6aSmpla4T48ePdi+fTv79u0r3bZu3TrcbjctWrQ4xrRFHC47G5odXBjyyCNwzTX25iMSwgKeek1LS2Pq1KnMmDGD1atXc/PNN5Ofn8+ggwsNRo8ezbXXXlt6/6uuuor4+HgGDRrEqlWr+Prrr7nzzjsZPHgwdbS0XeRwhYXQoIGJhw2Du++2Nx+pnN8P3mLwFtmdiQRJwOdRDhgwgF27djF27Fh27NhB586dmTdvXukCn4yMDDZv3lx6/3r16jF//nxuu+02unXrRnx8PP379+ehhx6qulch4hReL0RHm7hfPzh4GpbUYBs/h9cuN915hn5rdzYSBAGfR2kHnUcpIcHvN31b9+83xyZXr1b/1tpgw+fw6mWQcCrcvMDubBytVpxHKSJBYlnmCiD790NYGKxapSJZa5Q0Ra/xYw45RroepUhNcOml8MsvJi4oUJPz2qTVn2BMZtn5lOI4KpQidhs+HD780MRqcl77uD3mJo6lP4FE7PTEE/DMMybeu1dNzkVqIBVKEbu8/jrceaeJt2+HuDh785Fjs3UJPNYGpp5rdyYSJJp6FbHDZ5/B1VebeM0aaNrU3nzk2Pl9UJAFUVqR71QaUYpUt+XLoU8fEy9aVNamTmonNUV3PI0oRarTxo1w+ukmnjMHune3Nx85fokd4bZl4NEiLKdSoRSpLjt3Qps2Jn7ppbJLZ0ntFh4F8W3szkKCSFOvItUhLw9KruP64IMwWBf5FaktVChFgq2oCErabQ0ZAvfcY28+UrWy1sN//gSvXGp3JhIkmnoVCSafD6KiTNy7N7zwgr35SNXzFkLGj1Avwe5MJEg0ohQJFsuC+HgTt2oFn35qbz4SHFr16ngaUYoEyxlnQE6OidevV5Nzp4prCQP/C54IuzORIFGhFAmG/v1h6VITFxerybmTRdYzjdHFsfTbK1LV7roLZs82cX6+mpyL1HIqlCJV6Zln4PHHTbxnD0RH25uPBF9uBrxxJcweZHcmEiSaehWpKrNnm0tmAWzdCg0a2JuPVA9vIaydCxH17c5EgkQjSpGq8OWX5rgkwKpV0Ly5relINdKqV8fTiFLkeK1cCeecY+IFC6B9e3vzkeoV3Qgu+w+49XHqVPrJihyPTZvgtNNM/MEHcNZZ9uYj1S+yHnS6wu4sJIg09SpyrHbvhuRkE7/wAvz1r7amIyLBoUIpcizy86FRIxOPHWt6uEpo2p8Dc+6AOSPtzkSCRIVSJFDFxVCvnokHDoRx4+zNR+x1YD/88KK5iSOpUIoEwu+HyEgTn302TJ9uazpSA5SsesUy/X3FcbSYR+RoWRY0a2biZs3gq6/szUdqhoi6cO696uXrYCqUIkerZ0/IzDTxli36YBQjIhr+pOOTTqapV5Gjcc018N13JlaTc5GQohGlyB+55x547TUTq8m5/C9vESycbDrz9BgBHn2sOo1+oiKVmTIF/v1vE+/erSbncjhfMaQfXPmceqsKpQNp/kjkSD74AG6+2cSbN0PDhvbmIzWT65CPUfV7dST96SNSkQUL4NJLTfzzz5CUZG8+UnO5w6HbYFMw3R67s5EgUKEU+V+rV5sVrgDffAOnnGJvPlKzhUXARRPszkKCSFOvIofatg06dDDxO++UFUwRCVkaUYqU2LsXWrQw8eTJ8Le/2ZuP1A5+P6ycDVhwymUQFml3RlLFVChFAAoLyxbrjB4Nt9xibz5Su7x3o/nato8KpQNp6lXkwIGy0z6uvBIeftjefKR2ObRDk1a9OpJGlBLa/H6IiDBxSgq8/rq9+Ujt43LBSReYWOdQOpJ+qhK6LAtatTJxfDwsXGhvPlJ7XfWm3RlIEGnqVULXeeeZRgIAO3eqybmIVEgjSglNN9wAX3xh4qIiNTmX47NpIfi90OIMCI+yOxupYvp0kNDz4IPw0ksm3rev7BilyLF69TKYcRHk77Q7EwkCFUoJLdOmwdixJs7Kgrp17c1HnKGk36tWvTqSpl4ldMyZA9dfb+JNm8wCHpGqkHgqHCgEtz5SnUg/VQkNixbBRReZ+Kef4IQT7M1HnOX6T+3OQIJIU6/ifOvWwZlnmviLL6BjR3vzEZFaRSNKcbYdO6BdOxPPmgW9etmajjjU3k1m1WtsC7WwcyCNKMW5cnOhaVMTT5wI/fvbm48419Rz4NnTYc9GuzORIFChFGfavx9iY008ciQMH25vPuJsWvXqaJp6FefxeqFOHRNffjk8/ri9+YjzRceblogujT2cSIVSnMWyIDzcxJ07w9tv25uPhIZhi+zOQIJIf/6Is5x8svlarx4sW2ZvLiLiCBpRinNceKE5FQQgO1tNzqX6eIvN8UlPhPoGO5B+ouIMw4bB3LkmLioCj8fefCS0TOoK/06A7cvtzkSCQIVSar9HH4XnnjNxXp6anIsNDs5eaNWrI6lQSu326qswapSJd+0yxyZFqlvpalfL1jQkOHSMUmqvTz+Fa6818caN0KiRvflI6Lrle3NM3B1udyYSBCqUUjstWwZ9+5p4+XJo1crefCS06WLNjqapV6l9Nm6Erl1N/Nln5nxJEZEgOaZCOXnyZJKTk4mKiiIlJYXFixcf1X4LFiwgLCyMzvpgk2O1axe0aWPimTPhvPPszUcEYNoF8Fhr2PSd3ZlIEARcKGfNmkVaWhr33Xcfy5Yto1OnTvTt25edO3dWul92djbXXnst5+mDTY7Vvn3QpImJH38crrrK3nxESuzPhoLd4Cu2OxMJgoAL5VNPPcWQIUMYNGgQHTp0YMqUKURHRzNt2rRK9xs6dChXXXUVqampx5yshLDiYqhf38S3324anYvUFGqK7mgBFcri4mKWLl1K7969yx7A7aZ3794sXLjwiPu9/PLLbNy4kfvuu++onqeoqIjc3NxyNwlhPh9EHrzG31//Ck8/bW8+Iv/ryjfh9uWQdKbdmUgQBFQos7Ky8Pl8JCQklNuekJDAjh07Ktzn119/ZdSoUbz22muEhR3dItvx48cTGxtbektKSgokTXESyyo7N7JDB/jgA3vzEalIXBI0bA0R0XZnIkEQ1FWvPp+Pq666inHjxnHSSScd9X6jR48mJyen9LZly5YgZik1WufO5tqS4eHw8892ZyMiISig8ygbNWqEx+MhMzOz3PbMzEwSExMPu39eXh5Llixh+fLl3HrrrQD4/X4syyIsLIxPP/2Uc88997D9IiMjiSyZapPQddll8NNPJi4sVJNzqbneGwqZP8NfHoHknnZnI1UsoBFlREQEXbt2JT09vXSb3+8nPT29wkU6MTExrFy5khUrVpTehg4dSrt27VixYgUpKSnH/wrEmdLS4P33Tawm51LT7V4PO1bCfq2ncKKAO/Nn/YkfAAAgAElEQVSkpaUxcOBAunXrRvfu3Zk4cSL5+fkMGjQIMNOm27Zt45VXXsHtdnPqqaeW279JkyZERUUdtl2k1IQJ5gaQm6sm51LzadWrowVcKAcMGMCuXbsYO3YsO3bsoHPnzsybN690gU9GRgabN2+u8kQlRMyaZUaTAJmZZaeEiNRkFzwGRXnQpL3dmUgQuCzLqvHt7nNzc4mNjSUnJ4eYmBi705Fg+eILKDlmvX59WQceERHsqwXq9So1w48/lhXJJUtUJEWkxlChFPtt2lTW2HzevLKG5yK1xecPwetXwO8L7M5EgkCFUuy1ezckJ5t4xoyyS2eJ1CZbFsO6jyF3u92ZSBCoUIp9CgrKLrb88MNlF2EWqW1KVr1S45d8yDHQhZvFHgcOQN26Jr75Zhg92t58RI5Hj+HQ6UpI6m53JhIEKpRS/fz+snMj+/aF556zNx+R49XmHLszkCDS1KtUL8uChg1N3Lq1WbwjIlKDqVBK9UpJgZwcE69fb28uIlVlyTT4KA02f293JhIEKpRSfa68En74wcRer5qci3P8+hkseQl2rbE7EwkCFUqpHqNHw5tvmnj/fjU5F2cp+aNPvV4dSYt5JPgmT4ZHHjFxTg7oEmriNKf1h+ZdzU0cR4VSguvdd+HgtUjJyAD16hUn6nCJ3RlIEGnqVYLn22/h8stNvHYtVHBxbxGRmk4jSgmOVavg7LNN/P33cNJJ9uYjEkxr5sKu1dC6l6ZfHUiFUqre1q1wyikmnjPHnBIi4mQ/vwM/vw1hdVQoHUhTr1K1srMhKcnEL70E/frZm49IdSjp9apVr46kEaVUnf37oUEDE48bB4MH25uPSHVp3QuiYiDxVLszkSBQoZSq4fVCnTomHjwYxo61Nx+R6tTlanMTR9LUqxw/y4LwcBOfc46ZchURcQiNKOX4NWtmvrZoAZ9/bm8uInbY8gPs2QCJHSHhFLuzkSqmEaUcn7PPhh07TLx5s725iNhl+avw3k2wdq7dmUgQqFDKsRs40DQVADU5l9BWuurVsjcPCQoVSjk2990Hr7xiYjU5l1CXeCq06wcNW9udiQSBjlFK4KZOhQceMHF2tpqci5xxg7mJI2lEKYH573/hxhtNvG0bxMbam4+ISJBpRClH7/vv4a9/NfHq1WWrXUVC3e4NkLsNYlto+tWBNKKUo7N2LaSmmnjBAjj5ZHvzEalJfngRZlwMy16xOxMJAhVK+WMZGWWF8f334ayz7M1HpKZRr1dHU6GUyuXmlk2xPv88XKIL1IocJqYZJJ4G9ZvanYkEgY5RypEVFZUt1rnnHhg61N58RGqq1GHmJo6kEaVUzOeDqCgTX3MNPPigvfmIiNhEI0o5nGVB2MH/Gj16lDUWEJGKFeyB/dkQGQN1G9mdjVQxjSjlcK1ama+NGpW1qBORI/v+OXimC3z1mN2ZSBCoUEp5ffrApk0m3rnT3lxEagutenU0FUopc9NN8NlnJj5wQE3ORY5WeDREx0NEtN2ZSBDoGKUYDz8ML7xg4sLCsmOUIvLHeo4wN3EkjSgFpk+HMWNMvHdv2WpXERHRiDLkffwxDBpk4i1bIC7O3nxEaiO/D/xewAVhEXZnI1VMI8pQtmQJ9Otn4p9/hhYt7M1HpLZa8DQ81AQ++pfdmUgQqFCGqg0b4IwzTPz113DKKfbmI1KblS58s2xNQ4JDhTIU7dwJbdua+O234eyz7c1HpLbT6SGOpmOUoWbfPkhIMPGzz8Lll9ubj4gTnDkMUoaCy2N3JhIEKpSh5MABqF/fxHfdBbfeam8+Ik7hCUMfp86lqddQ4fdDxMHVeP37w6OP2puPiEgtoUIZCiwLPAenhLp2hVmz7M1HxGmWvQKPtYYPNEvjRCqUoaBDB/O1fn1zSoiIVC1vERTshqI8uzORIFChdLqLLoI1a0yck2NvLiJOVXJ6iFa9OpKOPjvZbbfBnDkmVpNzkeA59e/Q+hyIqGd3JhIEKpRO9fjjMGmSiQsK1ORcJJjqxJmbOJKmXp3o9dfN6R8Au3dDnTr25iMiUoupUDpNejpcfbWJN22Chg3tzUckFKyZC1POhrl32p2JBIHm45xkxQro3dvEP/4IJ5xgbz4ioaJwL+z4Ceol2J2JBIFGlE7x++/QpYuJP/8cTjvN1nREQkpJr1c1RXckjSidYPduaNXKxG++CeecY28+IqGmzbkw8COo08DuTCQIVChru4ICaNTIxE8+CQMG2JuPSCiqn2Bu4kiaeq3NvF6oW9fEI0ZAWpq9+YiIOJAKZW1lWRAebuJLL4UJE+zNRySUbf4eXr8CPrvf7kwkCDT1WltFR5uvp54K771nby4ioS5vB6z7GIpy7c5EgkAjytqoUyfYv9+MKH/6ye5sRKRk1at6vTrSMRXKyZMnk5ycTFRUFCkpKSxevPiI93333Xfp06cPjRs3JiYmhtTUVD755JNjTjjkXX55WXEsKlL/VpGaoFkXuOw/8KeRdmciQRBwoZw1axZpaWncd999LFu2jE6dOtG3b1927txZ4f2//vpr+vTpw9y5c1m6dCnnnHMOF198McuXLz/u5EPOyJHw7rsmVpNzkZojLgk6XQFte9udiQSBy7KsgM6QTUlJ4YwzzmDSwYbbfr+fpKQkbrvtNkaNGnVUj3HKKacwYMAAxo4de1T3z83NJTY2lpycHGJiYgJJ1zmeftqsbAXIzy87RikiEiLsqgUBjSiLi4tZunQpvXuX/dXkdrvp3bs3CxcuPKrH8Pv95OXl0bCSHqRFRUXk5uaWu4W0t94qK5K7dqlIitQ0O1fDnDvgm6fszkSCIKBCmZWVhc/nIyGh/Im1CQkJ7Nix46ge44knnmDfvn3079//iPcZP348sbGxpbekpKRA0nSWL78sayKwcWNZcwERqTlytsIPL8KqD+zORIKgWle9vv7664wbN4633nqLJk2aHPF+o0ePJicnp/S2ZcuWasyyBlm5sqwd3dKlZW3qRKRmKVkvoFWvjhTQeZSNGjXC4/GQmZlZbntmZiaJiYmV7vvmm29yww03MHv27HJTtxWJjIwkMjIykNScZ8uWssbmn34Kp59ubz4icmQN28B5Y3X1EIcKaEQZERFB165dSU9PL93m9/tJT08nNTX1iPu98cYbDBo0iDfeeIMLL7zw2LMNFXv3ll0i69VXoU8fe/MRkco1bAVn3wFd/ml3JhIEAXfmSUtLY+DAgXTr1o3u3bszceJE8vPzGTRoEGCmTbdt28Yrr7wCmOnWgQMH8vTTT5OSklJ6LLNOnTrExsZW4UtxiP37yy62/Mgj8E/94omI2CngQjlgwAB27drF2LFj2bFjB507d2bevHmlC3wyMjLYvHlz6f1feOEFvF4vw4YNY9iwYaXbBw4cyPTp04//FTiJzwd16ph42DC4+2578xGRo5O9BVa+ZS6z1W2w3dlIFQv4PEo7hMR5lJYF7oMz4RdcAHPn2puPiBy93xfA9H7Q6CS49Qe7s3GsWnEepQRRyXRr27YqkiK1jXq9OpquHlITpKRAdraJ162zNxcRCVz9BDjjBqjb2O5MJAhUKO129dVQ0lTe71f/VpHaqGFruPBJu7OQINHUq53+7//g9ddNXFysIikiUgNpRGmX556D8eNNnJdnri0pIrVTwR749VMIrwMdLrE7G6liKpR2eO89c/oHQGYm1Ktnbz4icnyyN8F7N0FMcxVKB9LUa3VbsAD+9jcTr18PlfS8FZFaQqteHU0jyuq0ejX07GnixYuhTRt78xGRqhEVCyddANFHvnyg1F4qlNVl+3bo0MHEc+fCGWfYm4+IVJ0GyXDVm3ZnIUGiqdfqkJMDzZubeNo003lHRERqBY0og62oCOLiTPzAA3CwebyIOEhxAWxbCm4PtDzL7mykiqlQBpPfD1FRJr7hBrj3XnvzEZHgyN0GMy6CqDgYtcnubKSKaeo1WCwLPB4Tn3ceTJ1qbz4iEkQHm4XU/GtMyDHQiDJYmjUzX5OS4LPP7M1FRIIrLBIST4MInRPtRCqUwfDnP8PBC1SzSdMwIo4XlwRDv7E7CwkSTb1WtUGD4OuvTezzqX+riEgtpxFlVbr/fpg+3cRFRWUXYhYRZ/N5TRs7gHg1EnEaFcqq8uKLMG6ciXNzISLC3nxEpPrsy4RnTwd3OIzNsjsbqWIa8lSFjz6CIUNMnJEB9evbm4+IVK+SXq9o1asTaUR5vBYtgosvNvHatZCYaG8+IlL93B6Ijge3PlKdSD/V47FuHZx5pokXLoSTTrI3HxGxR70mcNdGu7OQINHU67HKzIR27Uz84YdlBVNERBxFI8pjkZdXNsX6n/+UTb2KSGiyLPAVm+tRhkXptDCH0YgyUAcOQEyMie+9F2680d58RMR++7PhoSbw70Twe+3ORqqYCmUgLKvstI9rrzVXAxERcR3yUWr57ctDgkKFMhAlDQR69IAZM+zNRURqkEOmWtUY3XF0jPJotWplvjZuDN9+a28uIlKzRNaHe3YCLvCE252NVDEVyqNx/vnw++8mzsy0NRURqYFcLnMFEXEkTb3+kaFDYf58E6vJuYhIyFGhrMzDD5vTPwD271eTcxGpmLcYHmsNj7aC/Tl2ZyNVTFOvRzJjBowZY+LsbIjUtIqIHIHLDQW7TaxVr46jIVJFPvkErrvOxFu3QmysremISA1X7vQQrXp1Go0o/9fSpfCXv5h41Spo3tzefESk5nO54LZlpmBG6Q9rp1GhPNTGjdCtm4m/+Qbat7c3HxGpHVwuXbDZwTT1WmLXLmhz8D/6O+9Az5725iMiIjWCCiVAfj40aWLiSZPgb3+zNx8RqX2mngvP94R9O+3ORKqYpl69XqhXz8R33w3Dhtmbj4jUTjtWmiuI+IrtzkSqWGiPKC0Lwg+2mxowAB55xN58RKT2Kln5qtNDHCe0R5Qej/l6xhnw5pv25iIitds/3zFf6zaxNw+pcqFbKNu3NyPKmBhYvNjubESktkvWAkCnCs2p14svhjVrTJydbW8uIiJSo4Veobz9dvjoIxN7vWpyLiJV450b4PUBkLPN7kykioVWoXziCXj2WRMXFpYdoxQROV7rP4N186A43+5MpIqFTqF8/XW4804T79kDUVH25iMizqJVr44VGot50tPh6qtNvHkzNGhgbz4i4jwXPgm+A1A/0e5MpIo5v1CuWAG9e5v4p58gKcnefETEmU65zO4MJEicPfW6aRN06WLiL76Ajh3tzUdERGod5xbK3bshOdnEs2ZBr152ZiMiTpf+AHz0L8jeYncmUsWcWSgLC6FRIxNPmAD9+9ubj4g434+zYMk0yN9ldyZSxZxXKH0+iI42cVoajBhhbz4iEhpKV71a9uZhky+++ILnn3/e7jSCwlmF0rIg7OD6pMsugyeftDcfEQkdPW6H8+6DmGZ2Z3JEycnJTJw4scofd+fOndx444289tprvP/++1X++HZzVqGsW9d87dgR3n3X3lxEJLR0H0Je5yEQ09SWp//9999xuVylt/j4eM4//3yWL18e9Oe+9dZbefzxx3nnnXe4//772bNnT9Cfszo5p1B27myOTUZEmNNARESq0csvv0x8o0Y89thjWDZOv3722WdkZGTwySefkJ+fzwUXXEBOTk5Qn/Ott97i0ksvJTExkRUrVtCwYcOgPl91c0ahvPxy+PFHE+/fb28uIhKSGm36iJEpLib9exRXXX01hYWFf7hPr169uPXWW7n11luJjY2lUaNG3HvvvcdVaOPj40lMTKRbt2488cQTZGZm8v3335f+e15eHldeeSV169alefPmTJ48+Zifa/DgwZx22mkUFRUBUFxcTJcuXbj22mtL73P33Xdz0kknER0dTevWrbn33ns5cODAMT/nsZg3bx49e/YkLi6O+Ph4LrroIjZs2HDU+9f+QjlyZNk0q5qci4hNUlnGw+dG0vnCfzL77Xc5q0dPtm7d+of7zZgxg7CwMBYvXszTTz/NU089xYsvvgjA0KFDqVevXqW3ykQdbNVZXFxcuu3xxx+nU6dOLF++nFGjRjF8+HDmz59f+u8XXHBBpc93yimnlN73mWeeIT8/n1GjRgEwZswYsrOzmTRpUul96tevz/Tp01m1ahVPP/00U6dOZcKECUfxjlad/Px80tLSWLJkCenp6bjdbi677DL8/qNrN+iy7JwjOEq5ubnExsaSk5NDTExM2T88/XTZqtb8/LLVriIi1WzX/a1pzG6uKL6HrzOi2PP+Q9QPhw/ef4+zzjqrwn169erFzp07+eWXX3Ad/CN/1KhRfPjhh6xatYqdO3eSm5tb6fO2bdsWMMcoW7VqxfLly+ncuTPZ2dkMHjyY+fPns379ehISEkhOTqZ9+/Z8/PHHpftfccUV5ObmMnfuXAC2bdtW6Wg4PDycli1bln6/cOFC/vznPzNq1CjGjx/PF198Qc+eR7425xNPPMGbb77JkiVLKn1dFTliLQhQVlYWjRs3ZuXKlZx66ql/eP/a28Ju9uyyIpmVpSIpIrZa5enAqqXfktm+AREJTWn8z6fY8+Ej/PnPvZgy5Xmuv/76Cvc788wzS4skQGpqKk8++SQ+n48mTZrQpEmTgPI466yzcLvd5Ofn07p1a2bNmkVCQkK5xz9UampquZWwzZs3D+j5UlNTGTlyJA8++CB33333YUVy1qxZPPPMM2zYsIF9+/bh9XqPq8gdi19//ZWxY8eyaNEisrKySkeSmzdvPqpCWTunXr/6qqyJwG+/QXy8vfmISMhbGN6TEekufrPMqldP3Tga9X+QqFPO44YbbuC2224L+NjcsUy9zpo1ix9//JG9e/eyYcMG+vXrF9BzBjL1CuD3+1mwYAEej4f169eXf08WLuTqq6+mX79+fPTRRyxfvpwxY8aUmwquDhdffDF79uxh6tSpLFq0iEWLFgEcdR61b0T5889l7eiWLy9rUyciYiO3233YJbZcnnDi/3IrEQmtee7551n58y+88/Zs4g/5477kQ7vE999/z4knnojH4+GBBx5g5MiRAeWRlJREmzZtjvjvhy7sKfm+ffv2pd+/+OKLfzj1eqjHH3+cNWvW8NVXX9G3b19efvllBg0aBMB3331Hy5YtGTNmTOn9N23aFNDrOV67d+9m7dq1TJ06lbPPPhuAb7/9NqDHOKYR5eTJk0lOTiYqKoqUlBQWL15c6f2//PJLTj/9dCIjI2nbti3Tp08/lqeFrVvLGpvPn29OCRERqQHa+DdwdXuLRHYf9m/1u/Sjcf+H+O6HZZzetRsrV64s/bfNmzdz+4gR/LhyFa+8NpNnn32W4cOHU+T1EV43jgaJSTRp0ZK2bdvStm1bWrRsRZ1GzYls2Ix6jVuUPk6R1wfApt35rN+5r3S71+fnxy3ZLN+8l2KvGf099thjrFu3jrseeIzZs2fT+x/Xsf+A2b958+ZkuRvwe3F9NhbVo1GzE0qfO8vdgB+zw/lgxTa27i1g+fLljB07lrv/PYGtEUlccds93Hb7cDZu3AiAJ64pmzZt5qb7n+HtL5bwzDPP8N577+G3YPIX65n0+a98vDKjNNffsvJ54pO1PDZvDa8u/L10+87c/Yz77y+Mn7s64J9LgwYNiI+P54UXXmD9+vV8/vnnpKWlBfQYAY8oZ82aRVpaGlOmTCElJYWJEyfSt29f1q5dW+Fc+m+//caFF17I0KFDmTlzJunp6dxwww00bdqUvn37BvbkJUP+114ru3SWSA1iWRZ+C/yWhd+yiAzzlP7bviIvfsvC8kO9qDA8bnNcave+Ig74LCwsGtaNKN1nR85+Coq9+C1IiImkflQ4lmWxdU8+u/cV4fP7aVwvnPjocPx+P9uzC9mevR+v30/jumE0i4nA7/eTmbufDVmF+Px+4qI8tI2PxO/3szu/mJ93mO11w12clmC25+73smhrAZYF4W6LM5uZxyn2+vh8UxF+vx+/BeclefD7/fh8Pj7b7KXY68fnt0htfIAwy4vP5+O7nR5yi8370bFuHvVcxXi9Xn7cV4/dB8Lw+y3auLOIYx8+n49fDzRgl1Ufn9+iqTeDhsU78fl8ZLibkBmeiN+yiC3YRsO8DXi9XrIjE9gV1x7Lgsh924ndvhiv10tBVCNyWv4Jy3Lhzt9F5M8f4PP5OBARQ3HXq8CysPbv48AXk/B6vfhcHupcdA+4XFjArldG4PV6sSyLxIETcUdEg8tF5sy78OXvBSDhiocJi2+By+Vi59sPsODi7Uy70MN1xZvZ4T/8cFDUCR1pOuw1XC4XvYY+xLcv3g9Ap3P+ystfrWPS82cQFRHGv4YP58Ybb+StJVu4+x1TUPt1TOS5q7sCsOi3PQycZgYnnZPieH9YDwC27zWjwJtfW8YJJ+7n+/87D4DCAz4umbwAgJ15Rdx7cPXnuHHjKHZHEtPrep5ZW5cr8otpFlcHgHvf/5lfDxbbd24+i64tIwBT3L75NQuApy5vz73X/ZPrrruOzLhTeeadlUBH2nZK4ZprruHrr7/GOqEbdbv+lRcfu5fpj3m57K8Xc++99zLm3vt4/JO1pa/tgo5munrzngImfbG+9LVdk5oMQHbhAV5e8Dv+ooIAfyPNSP/NN9/k9ttv59RTT6Vdu3Y888wz9ArgQhkBF8qnnnqKIUOGlA6tp0yZwpw5c5g2bVrpEuFDTZkyhVatWvHkwXZy7du359tvv2XChAmBF0qARx8tuwhzNSpZHFxy0P2Az4/XZz4MPW4XUeHmw23/AR/5RebDLcLjJjY6vHT7rrwi8+ET5qJprPkPWeT18XtWAX7Ljwto0ygav9/PAa+Pldty8Pr8+P1+TmtWD8uy8Pl8LNmUTbHXh9fn45QmUXgwH1bLtu5jX5EXn99PuwYe6nj8eL1eft5ZzJ5CL16fn9b1/cSFm+3rsy0yCiy8Ph/NI4tpHF6Mz+djc76brYXh+Px+4t2FJLpz8fl87CyOYNOBevgsi3r+fJr7MvF6veT4I9nkborfsojwFtB03zp8Ph8F/jC2xpyC37Jwe4tovON7fD4fxT6LzBa9sDAfoDGrP8DnNR+s+zpcgt8djt8C19JZ+Atz8Xq9uLoNwKobj2VBwXevcSBrM16vl7pnXUV44kngcpHz7evs32yaTcSk9if6xDPB5SZvyQfk//IFAPW69KP+6RfhcrnY9/Pn5H4/G4Dok3vSoNdgcLko3LiUPZ+Y5e2RJ3Sk8SWjwOWmOHMDO2fdA0BYwxY0G/QM4MKXv4dtU8xCDVdEHU741+zS/zebHr2oNE4aMQt3pOketfW5QfjyTPPspoMnE9HYrCLc8dpIiratAaDJP8ZRp7X5cNz1waMUrPkGgIZ/uY36nczvzp7PXyLvh/cAiO1xJXE9ze9G7g/vs/dzc4pB3Y59aNRvOAD5axeQ9f54AKKSu5Aw4EHz/3D7Gna8aqb4wuOTaHaD6dnpzdvNtucGVvjaHr6u4tc2ZfT1+PKyDnttbz87rsLX9vkH0w5/bR5Y/vU8chcffG1nXUHc2SkArF/1I3vT3zrktV1o3rsdGWR9/vnB19aZhC6DzWsr8rLj4AguLL4Fzfu2BsCXt5uMjIzS11a/cXLp6zlwwAuY3/nwuETcUQePBbrLPjLd0bGE1TMn1rs8YVz0egHhHoi+ufwxvEOVfH74fb7S42OesDDi+w4jvu8wrjsrmfv/avZ3UbbA59CzGMLcLiLD3Lhd5muJlsnJdB73CW6Xi7josulRj9tF87g6uFyQNOYN7h15Du6Df6RdMnkBhcVe3C5X6R9uAO0S6xMdGYbbBdERZX/sdWgWg9dn4XZD04ax/PLLLwBMX/Abvds3weVy8fdrXqXvKYmlj3PN8DG4XffQq11j/na6GQH3uPRa3l6yFZcLTm0eW/r4zePqcN1ZybhdLpo3qFO6vWHdCG7p1Yaign2MPYYOfL1792bVqlXltgVywkdAhbK4uJilS5cyevTo0m1ut5vevXuzcOHCCvdZuHAhvf9n9Ne3b19GVNKsvKioqPQEVqC0q0Sby+7B+m4vu0tWiHkiSLjiIXPupMtF5qt3UvKfu8mVj+COjMblcpE5eyz+feavwMaXjyU83vywdn3wKAcyzUmnDS9Mo07LTuByseeTyRSuN8cNGpw3hHod+wCQ/c0r5C39CICY7pcTe9YAAPKWzyX7q+kARHfoRfz5twCQ/+v37JnzFGA+cJv87V7z+jLWs3PW/5mX0KAZzQaan7x33x4yXhxqXltYJEm3vlr6HmyZWHYFlOZDp5X+4m576Wb8eWaqJ+GfjxPRyHwoZc66h+KMdQA0unQUdZJPByBr7kQK13138LXdSL2OvQEX2d/MOsJr++qQ13YO8efffPC1rS97bUmn0uTyseCCoqz1fDprxiGvrV/pa1s8Z87B1xZB0q1lKwDXLFpc+nNrdubteEpe27btZa/tnKalr62g2E9xfj4A9WMSiEgwx2NcUWULGzzRcYQ3ND9n9yHb3eF1CKtvrixT8sFuvgkrvZ8rLJJDuTwlHzplH0qW5cfyeQHwH/wKYPl9//NXr6v0tfn25x/8HizKfkn9+/fhK8g1o5xDPhG9BTl4c3ZiWRbWgbLfByt/Lwd2bwUswv3FxMXF4Xa7iXAdwJf1O2BR3+OjQevWhIWFYTWoy4Gsg//Pw4ppdcYZhIWF4a3fjL271+N2QYw7ny4XXIDH48EXWZ/N2RtwuyDCOsC511yDx+MBTzg/5/+OywVul4tLR4wgPCwMj8fDD77dWP5s3C4X/7zzDuqFg8fjYcm+cAqsnXhcLlLu+heN6rjweDysyI5k94E8PC4Xp991G0kxppCv2gvb8vfjcrlof/tA2sRdj9vtZkO2j99yfbhcLpK7/I0O/3cFLpeLHfkWa7KKcbvdJJzRm053XYzb7Sa70MdPu4rwuFzERHWn67+W4Xa72e/182NGIS63m6iwxnS/fQ0ulws/LpZv3YfH7WYD8twAAAi1SURBVMLtctFt6EY8Hg8ul4tVO/KxAI/Hzcm3/0RURBgul4vfdxfgtUwhajHyW16Z9iKj7xlL8yIvUPZ/osSB7AxyPp1MpK+Ql180Awifz0f7JnV4Pe1M3C4X4R536ekg558Yw/mjeh78eHOVbu/YJIIf7upR+rgl2xtFwlcjzjxsO8DHt3QrjfftyyuNX73m0Gv0FpOba4r3wxe1LZd7yWMN69EcaH7Y9r91jOdvHeMP296zZV16tmx72PZ2DcMYc37yYdubREFar6TDtkcAQ89qRm5uLmMJrMhVCSsA27ZtswDru+++K7f9zjvvtLp3717hPieeeKL18MMPl9s2Z84cC7AKCgoq3Oe+++6zMJ8uuummm2666VbutmXLlkBK13GrkateR48eXe5ga3Z2Ni1btmTz5s3ExsZWsmdoyc3NJSkpiS1btlT7eUk1nd6biul9OTK9NxWrSe+LZVnk5eXRrFn1XqEloELZqFEjPB4PmZmZ5bZnZmaSmJhY4T6JiYkV3j8mJoY6depUuE9kZCSRkZGHbY+NjbX9B1UTxcTE6H05Ar03FdP7cmR6bypWU94XOwZLAZ0eEhERQdeuXUlPTy/d5vf7SU9PP6zbQ4nU1NRy9weYP3/+Ee8vIiJSkwR8HmVaWhpTp05lxowZrF69mptvvpn8/PzSVbCjR48u1zl+6NChbNy4kbvuuos1a9bw3HPP8dZbb/Gvf/2r6l6FiIhIkHjuv//++wPZ4dRTTyUuLo5///vfPPHEEwDMnDmTdu3aAfDaa6+xadMmrrvuOsCc7NmjRw+ee+45HnzwQX755ReefPJJ/v73vweWqMdDr169CAurkYdVbaP35cj03lRM78uR6b2pWKi/L7Xi6iEiIiJ2qZ1N0UVERKqJCqWIiEglVChFREQqoUIpIiJSiRpTKG27dFcNF8j78u6779KnTx8aN25MTEwMqampfPLJJ9WYbfUK9P9MiQULFhAWFkZnh16mLdD3paioiDFjxtCyZUsiIyNJTk5m2rRp1ZRt9Qn0fZk5cyadOnUiOjqapk2bMnjwYHbvPvwSWrXd119/zcUXX0yzZs1wuVy8//77f7hPqHz+lqrWhnlH8Oabb1oRERHWtGnTrF9++cUaMmSIFRcXZ2VmZlZ4/40bN1rR0dFWWlqatWrVKuvZZ5+1PB6PNW/evGrOPLgCfV+GDx9uPfroo9bixYutdevWWaNHj7bCw8OtZcuWVXPmwRfoe1Ni7969VuvWra3zz///9u4vpKk2jgP4Nx0LCUq6SNdGF4KKCtWFTFREkHgDgyAQB4pIICEVCKEiKMpAEExEjXW7gpaSkeCFRl6IOKKLYkKlULQUJScUiBO9mt/3QhT/veftHDjPmdvvA7vw+AjffZnPT+fB5x9eu3ZNUVp1jPRy+/ZtFhUVcWpqij9//uT79+8ZDAYVpjaf3l6CwSBTUlI4ODjIcDjM2dlZFhQU8M6dO4qTm29iYoLt7e188+YNAXBsbExzfbLsvwfFxaB0u9188ODB/sexWIyXL19mT0/PietbW1tZUFBw6JrH4+HNmzdNzama3l5Okp+fT6/Xa0Y8SxntxuPxsKOjg11dXQk5KPX2Mjk5yQsXLvDPnz+qIlpCby+PHz9mVlbWoWtDQ0N0Op2m5rTa3wzKZNl/D7L8rde9o7sOHsVl9Oiu/1p/Ghnp5aidnR1Eo1FcvHjRrJiWMNqN3+9HOBxGV1eXipjKGellfHwchYWF6O3thdPpRE5ODpqbm7G9va0qtumM9FJcXIzl5WVMTEyAJNbW1jA6OorKykpVseNWMuy/R1k+KH///o1YLIaMjIxD1zMyMhCJRE78mkgkcuL6jY2NhPkGN9LLUX19fdjc3ER1dfX/Lz5FjHTz/ft3tLW14cWLFwn730WM9BIOhxEMBvHlyxeMjY1hYGAAr1+/xv3791VEVsJIL6WlpQgEAvB4PLDb7cjMzER6ejp8Pp+KyHEtGfbfoywflMIcL1++hNfrxatXr3Dp0iWr41gqFouhpqYGXq8XOTk5VseJKzs7Ozhz5gwCgQDcbjcqKyvR39+P58+fJ+ym9zfm5+fR1NSEzs5OfPr0CW/fvsXi4iIaGxutjiYsYPmP1qqO7jptjPSyZ2RkBA0NDRgdHT32Fkki0NtNNBrFx48fEQqF8PDhQwC7A4IkbDYb3r17h4qKCiXZzWTkNeNwOOB0Og8dXZSXlweSWFlZQXZ2tqmZVTDSS09PD0pKStDS0gIAuHr1Ks6dO4eysjJ0d3fD4XCYnjteJcP+e5Tlv1HK0V0nM9ILAAwPD+Pu3bsYHh7GrVu3VERVTm8358+fx+fPnzE3N7f/aGxsRG5uLubm5lBUVKQyvmmMvGZKS0vx69cvbG5u7l/79u0bUlJS4HK5TM+sgpFetra2jr1Fn5qaCmD38OBklgz77zHW3ku0a2RkhGfPnuWzZ884Pz/Pe/fuMT09nZFIhCTZ1tbGurq6/fV7tye3tLRwYWGBPp8vIW9P1ttLIBCgzWajz+fj6urq/mN9fd2qp2Aavd0clah3vertJRqN0uVysaqqil+/fuXMzAyzs7PZ0NBg1VMwhd5e/H4/bTYbnz59yh8/fjAYDLKwsJBut9uqp2CaaDTKUCjEUChEAOzv72coFOLS0hLJ5N1/D4qLQUmST5484ZUrV2i32+l2u/nhw4f9z9XX17O8vPzQ+unpaV6/fp12u51ZWVn0+/1qAyuip5fy8nICOPaor69XH1wBva+ZgxJ1UJL6e1lYWOCNGzeYlpZGl8vFR48ecWtrS3Fq8+ntZWhoiPn5+UxLS6PD4WBtbS1XVlYUpzbf9PS05r6RzPvvHjlmSwghhNBg+d8ohRBCiHgmg1IIIYTQIINSCCGE0CCDUgghhNAgg1IIIYTQIINSCCGE0CCDUgghhNAgg1IIIYTQIINSCCGE0CCDUgghhNAgg1IIIYTQIINSCCGE0PAvzGKJy0yoLlUAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x329310e90>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x3296e27d0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = rand(2)                       # random red vector\n",
    "a = rand(2); a *= 1.1/maximum(a)  # target direction\n",
    "P = (a*a')/a'a\n",
    "p = P*b\n",
    "\n",
    "figure(figsize=(5,5))\n",
    "arrow(0,0,b[1],b[2],head_width=0.05, head_length=0.03,color=\"r\")\n",
    "text(b[1]+.03,b[2],\"b\",color=\"r\")\n",
    "plot([0,a[1]],[0,a[2]],\":\")\n",
    "axis([0,1.1,0,1.1]);\n",
    "\n",
    "plot([b[1],p[1]],[b[2],p[2]],\":\")\n",
    "arrow(0,0,p[1],p[2],head_width=0.05, head_length=0.03)\n",
    "text(p[1]+.03,p[2],\"p=Pb=x̂a\")\n",
    "text(a[1]+.03,a[2],\"a\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.999947    0.00729359\n",
       " 0.00729359  5.31993e-5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.999947    0.00729359\n",
       " 0.00729359  5.31993e-5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.999947    0.00729359\n",
       " 0.00729359  5.31993e-5"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Powers of P remain equal.  Explain why geometrically?\n",
    "# Answer: once you project, projecting again keeps you where you were\n",
    "display(P)\n",
    "display(P^2)\n",
    "display(P^3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Relationship to least squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Array{Float64,2}:\n",
       " 0.770082"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " x̂= (a'b)/(a'a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Float64,1}:\n",
       " 0.770082"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection on a subspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 Array{Float64,2}:\n",
       " 0.857652  0.517855   0.280477\n",
       " 0.418829  0.556069   0.964192\n",
       " 0.376845  0.64954    0.692936\n",
       " 0.839149  0.0958249  0.297049\n",
       " 0.533046  0.988303   0.900709"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(5, 3) # consider the subspace spanned by the columns of A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.65487 \n",
       " 0.536363\n",
       " 0.847276\n",
       " 0.527487\n",
       " 0.667052"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = rand(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our problem\n",
    "1. Find the vector $p$ that is in the column space of $A$ that is closest to $b$\n",
    "2.  Project $b$ onto the column space of $A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Find the linear combination of the columns of ($m \\times n$) $A$ closest to $b$\n",
    "\n",
    "In other words, find an $\\hat{x}$ in $\\Re^n$ such that $A\\hat{x}$ is closest to $b$.\n",
    "\n",
    "How do we find $\\hat{x}$?  Idea is the same as the line. Make $e=b-A\\hat{x} \\perp $ to every column of $A$:\n",
    "\n",
    "$A^T(b-A\\hat{x})=0$ is equivalent to the first column of $A$ is orthogonal to $e$, and the second column is orthogonal to $e$, ..., and the last column of $A$ is orthogonal to $A$.\n",
    "\n",
    "$A^TA\\hat{x} = A^Tb$. (known as the **normal equations**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. $\\hat{x} = (A^TA)^{-1}A^Tb$\n",
    "2. $p = A\\hat{x} = A(A^TA)^{-1}A^Tb$\n",
    "3. $P =A(A^TA)^{-1}A^T $ (is the projection matrix)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 Array{Float64,2}:\n",
       " 0.857652  0.517855   0.280477\n",
       " 0.418829  0.556069   0.964192\n",
       " 0.376845  0.64954    0.692936\n",
       " 0.839149  0.0958249  0.297049\n",
       " 0.533046  0.988303   0.900709"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       "  0.680624   -0.275029   0.0718287   0.273674    0.24835 \n",
       " -0.275029    0.726659   0.227343    0.235613    0.125645\n",
       "  0.0718287   0.227343   0.233565   -0.0612794   0.344112\n",
       "  0.273674    0.235613  -0.0612794   0.765488   -0.212956\n",
       "  0.24835     0.125645   0.344112   -0.212956    0.593663"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = A * inv(A'A) * A'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       "  0.680624   -0.275029   0.0718287   0.273674    0.24835 \n",
       " -0.275029    0.726659   0.227343    0.235613    0.125645\n",
       "  0.0718287   0.227343   0.233565   -0.0612794   0.344112\n",
       "  0.273674    0.235613  -0.0612794   0.765488   -0.212956\n",
       "  0.24835     0.125645   0.344112   -0.212956    0.593663"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P^10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.458361\n",
       " 0.601126\n",
       " 0.721415\n",
       " 0.940134\n",
       " 0.816022"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = rand(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.658412\n",
       " 0.798797\n",
       " 0.561274\n",
       " 0.768752\n",
       " 0.721845"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = P*b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       "  0.200051\n",
       "  0.197671\n",
       " -0.16014 \n",
       " -0.171382\n",
       " -0.094177"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e = p - b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " -1.49186e-15\n",
       " -1.22125e-15\n",
       " -1.67921e-15"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A'e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  0.707741\n",
       " -0.265997\n",
       "  0.674437"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x̂ =  inv(A'A)*A'b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  0.707741\n",
       " -0.265997\n",
       "  0.674437"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A\\b # in matlab and in julia, to solve the least squares system\n",
    "# Ax=b for the best vector x̂, type A\\b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Math: $(A^TA)$ is invertible when $A$ has linearly independent columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that $A^TA$ is not invertible.  Then there is a nonzero x $x$ such that $A^TAx=0$.  Then\n",
    "$x^TA^TAx=0=\\|Ax\\|^2$.  Then $Ax=0$ meaning $A$ does not have linearly independent columns.\n",
    "Taking the contrapositive, if $A$ has linearly independent columns $A^TA$ is invertible.\n",
    "\n",
    "Note logically one should prove the converse too.  This is implied in the \"when.\" \n",
    "If $A$ does not have  linearly independent columns, there is a nonzero $x$ with $Ax=0$.\n",
    "Multiplying by $A^T$ we have $A^TAx$ is then $0$\n",
    "so $A^TA$ is not invertible.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Briefly mentioned:\n",
    "\n",
    "* Chebychev Approximation = polynomial fitting = linear equations\n",
    "\n",
    "* Machine learning = nonlinear fitting = nonlinear equations\n",
    "\n",
    "* In high school stats classes , students are told to divide by $n-1$, not $n$, for sample variance.\n",
    "\n",
    "  - Some argument about degrees of freedom usually appeases the masses. In fact, the projection matrix `P = I - ones(n,n)/n` can be viewed as \"removing the mean\" or projection orthogonal to the \"ones\" vector. Removing the true mean creates a vector whose element squares have expectation $\\sigma^2$ and cross terms have expectation 0.\n",
    "\n",
    "  - You might check that the sample variance numerator is $\\|Pb\\|^2$. This is the same as $b^TPb$, which is readliy checked to have average\n",
    "$\\sigma^2$ times the sum of the diagonal elements of $P$, which is  $n \\times \\left( 1-\\frac{1}{n}\\right)=n-1$.\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  },
  "widgets": {
   "state": {
    "fe29c940-49a2-4f03-9b57-b95c7d087987": {
     "views": [
      {
       "cell_index": 8
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}