{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# EART 70013  <a class=\"tocSkip\">\n",
    "    \n",
    "# Geophysical Inversion <a class=\"tocSkip\">    \n",
    "    \n",
    "## Lecture 4 - Homework solutions <a class=\"tocSkip\">  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Homework\" data-toc-modified-id=\"Homework-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Homework</a></span><ul class=\"toc-item\"><li><span><a href=\"#Homework---Row-operations-on-an-over-determined-problem\" data-toc-modified-id=\"Homework---Row-operations-on-an-over-determined-problem-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Homework - Row operations on an over-determined problem</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.1.1\"><span class=\"toc-item-num\">1.1.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li><li><span><a href=\"#Homework---An-over-determined-system-with-(by-construction)-an-exact-solution\" data-toc-modified-id=\"Homework---An-over-determined-system-with-(by-construction)-an-exact-solution-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Homework - An over-determined system with (by construction) an exact solution</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.2.1\"><span class=\"toc-item-num\">1.2.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li><li><span><a href=\"#Homework---Least-squares-solution-as-a-compromise-between-all-constraints\" data-toc-modified-id=\"Homework---Least-squares-solution-as-a-compromise-between-all-constraints-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Homework - Least squares solution as a compromise between all constraints</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.3.1\"><span class=\"toc-item-num\">1.3.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li><li><span><a href=\"#Homework---An-even-simpler-over-determined-case\" data-toc-modified-id=\"Homework---An-even-simpler-over-determined-case-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>Homework - An even simpler over-determined case</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.4.1\"><span class=\"toc-item-num\">1.4.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li><li><span><a href=\"#Homework---Outer-product\" data-toc-modified-id=\"Homework---Outer-product-1.5\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>Homework - Outer-product</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.5.1\"><span class=\"toc-item-num\">1.5.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li><li><span><a href=\"#Homework---Matrix-rank-and-RREF-(a-non-square-case)\" data-toc-modified-id=\"Homework---Matrix-rank-and-RREF-(a-non-square-case)-1.6\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>Homework - Matrix rank and RREF (a non-square case)</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.6.1\"><span class=\"toc-item-num\">1.6.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li><li><span><a href=\"#Homework---Minimal-norm-solution-to-under-determined-problem\" data-toc-modified-id=\"Homework---Minimal-norm-solution-to-under-determined-problem-1.7\"><span class=\"toc-item-num\">1.7&nbsp;&nbsp;</span>Homework - Minimal-norm solution to under-determined problem</a></span><ul class=\"toc-item\"><li><span><a href=\"#Solution\" data-toc-modified-id=\"Solution-1.7.1\"><span class=\"toc-item-num\">1.7.1&nbsp;&nbsp;</span>Solution</a></span></li></ul></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.linalg as sl\n",
    "from pprint import pprint\n",
    "import scipy.interpolate as si"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - Row operations on an over-determined problem\n",
    "\n",
    "Consider the following example from the lecture\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "  2x + 3y &= 7 \\\\[5pt]\n",
    "   x - 4y &= 3 \\\\[5pt]\n",
    "  -3x - 10y & = -11\n",
    "\\end{align*}\n",
    "   \\quad \\iff \\quad\n",
    "  \\begin{pmatrix}\n",
    "    2 & 3 \\\\\n",
    "    1 & -4  \\\\\n",
    "    -3 & -10 \n",
    "  \\end{pmatrix}\n",
    "  \\begin{pmatrix}\n",
    "    x \\\\\n",
    "    y \n",
    "  \\end{pmatrix}=\n",
    "  \\begin{pmatrix}\n",
    "    7 \\\\\n",
    "    3 \\\\\n",
    "    -11\n",
    "  \\end{pmatrix}   \n",
    "$$\n",
    "\n",
    "Use row operations on the augmented system in an attempt to solve this problem.\n",
    "\n",
    "<br>\n",
    "\n",
    "In doing this you will be able to also establish the rank of $A$ and the rank of the augmented matrix.\n",
    "\n",
    "What relationship do you need between these two values in order for the problem to have a solution?\n",
    "\n",
    "Can you think of an example where this wouldn't be the case and you have a system without an exact solution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "For the augmented system\n",
    "\n",
    "$$\n",
    "[A \\, | \\, \\boldsymbol{b}] = \n",
    "\\left[\n",
    "  \\begin{array}{rr|r}\n",
    "    2 & 3 & 7\\\\\n",
    "    1 & -4 & 3 \\\\\n",
    "    -3 & -10 & -11 \n",
    "  \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "Swap rows as an easy way to get a \"1\" top left\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "  \\begin{array}{rr|r}\n",
    "    1 & -4 & 3 \\\\\n",
    "    2 & 3 & 7\\\\\n",
    "    -3 & -10 & -11 \n",
    "  \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "use the \"1\" to set the values below to \"0\"\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "  \\begin{array}{rr|r}\n",
    "    1 & -4 & 3 \\\\\n",
    "    0 & 11 & 1\\\\\n",
    "    0 & -22 & -2 \n",
    "  \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "Now scale a row to get a \"1\" as the leading entry of the next row\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "  \\begin{array}{rr|r}\n",
    "    1 & -4 & 3 \\\\\n",
    "    0 & 1 & 1/11\\\\\n",
    "    0 & -1 & -1/11 \n",
    "  \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "use if to set entries above and below to zero\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "  \\begin{array}{rr|r}\n",
    "    1 & 0 & 3+4/11 \\\\\n",
    "    0 & 1 & 1/11\\\\\n",
    "    0 & 0 & 0 \n",
    "  \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "and so this problem has the unique solution $ x = 3+4/11 = 37/11$ and $y = 1/11$.\n",
    "\n",
    "Key here was that $\\text{rank}([A|\\boldsymbol{b}]) = \\text{rank}(A)$.\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "What would have been the outcome if we ended up with the final form\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "  \\begin{array}{rr|r}\n",
    "    1 & 0 & 3+4/11 \\\\\n",
    "    0 & 1 & 1/11\\\\\n",
    "    0 & 0 & 1 \n",
    "  \\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "i.e. a situation where $\\text{rank}([A|\\boldsymbol{b}]) > \\text{rank}(A)$?\n",
    "\n",
    "Our system would have had no solution as the three equations would have been inconsistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - An over-determined system with (by construction) an exact solution\n",
    "\n",
    "Recall the simple over-determined problem from the lecture\n",
    "\n",
    "`A = np.array([[2, 3], [1, -4], [1, 10]])`\n",
    "\n",
    "You were asked to think about how you could change the RHS vector only in order to come up with a version of the over-determined problem that has an exact solution. You were given the hint to think about the range of the LHS matrix $A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "All we need to do is ensure that the RHS vector is in the range of $A$. If it is then we can take an appropriate weighted sum of the columns in order to reach that point, these weights in turn are the solution to the problem $\\boldsymbol{x}$.\n",
    "\n",
    "As an example let's just take $\\boldsymbol{b}$ to be the sum of the columns.\n",
    "\n",
    "If our argument is right the three lines visualising the three equations/constriants should all coincide at a single location, e.g here that should be at the location $\\boldsymbol{x} = (1,1)^T$.\n",
    "\n",
    "Let's plot this situation to check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x24e10dc0d88>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAFSCAYAAABG/JyrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hVVdbH8e9KILTQSwi9xYINDWNDIaEJCGJBxQKWUWbGPupgey1jr1hGxzKKmqjEMhYYUUQMoI4NRkERIYCICIKEIqEkhKz3j31D2k1Iwr055yTr8zx5knvP4WYR5Zd999lnL1FVjDHGRFaM1wUYY0xtZOFqjDFRYOFqjDFRYOFqjDFRYOFqjDFRYOFqjDFRYOFqTBgicq+IzCvvsTF7Y+FqfEVEdC8fL3hU2p3ACR59bxNA9bwuwJhSEot9PRL4V6nndtRsOY6q5gA5XnxvE0w2cjW+oqq/Fn4Am0s/p6pbROQdEXm48M+IyEOhUe1hocciIr+JyGmhx7EicruIrBaRXBH5RkRGVKWuMNMEGSLyhoj8TUTWiki2iPxLRBoUOydGRG4SkR9FZIeILBSRM4odFxG5Q0RWhepaIyLPVvuHZ3zFwtUE0WwgtdjjFGBDsecOBloDc0KPJwJXANcAhwIzgHdE5MB9rGMI0C30fccBY4FLih1/ADgb+BPQG3gIeFFEBoeOnw1cClwMJAGjgfn7WJPxCQtXE0SzgUNFpJWINAcOAR6mKFxTgEWquiH0+FrgblV9VVWXqOp1uBC7Zh/r2ABcrqo/qOp04G1gEICItAAuB85X1Q9U9UdVfRF4kaIA7gr8AnyoqqtU9StVfXIfazI+YeFqgmgBbspgANAf+A54C+gvIjG4cJ0NICLtgFbAp6Ve4xPcaHJffKeqBcUerwHahb4+BKgPZIpITuEHcAHQM3RORqi2H0NTCqeJSP19rMn4hF3QMoGjqgUi8jFupLoLyFTVxSKSCxyOC9w/h06Xwj8W7qX2sZRdYV6vcMBS+HkY8Gup8/IAVHWFiPTCTS8MAh4DbhKRY1V15z7WZjxmI1cTVLNx4ZoS+hrcHOvlFJtvVdV1QDZwXKk/fxzwfRTr+xbIBzqr6rJSH6sKT1LVHao6VVWvBI7F/XI4Mop1mRpiI1cTVLNxF4jygbnFnvsHJedbAR4EbhSRH3FTChcCybi36FGhqhtF5FHg0dBb/U+BZrgA3aGqk0Xk4lD9XwHbcBfF8oDl0arL1BwLVxNUhfOuK1R1S+i5TCCWopFsoQeAxsAjQFtgMXCyqi6Oco0TcfOwNwI9QvV+DdwbOr4Zd7HtkVDdi4DRqvpLlOsyNUCsE4ExxkSezbkaY0wUWLgaY0wUWLgaY0wU+CpcRaSziGSKyGIRWSQiV4Y5J0VEtoTuD/9GRG7xolZjjKmI31YL5APXqOr/RKQpMF9EZqpq6fWIH6vqSA/qM8aYSvFVuKrqWmBt6OutIrIY6EgEFnu3adNGu3XrVunzt23bRpMmTfb129aIoNRqdUZWUOqE4NRanTrnz5+/QVXbljmgqr78wO02tApoVur5FNwdNwuA94CDKvN6ycnJWhWZmZlVOt9LQanV6oysoNSpGpxaq1MnME/DZI4v17mKSDzu9sW7VPXNUseaAQWqmhPak/NRVU0q53UmABMAEhISkjMyMipdQ05ODvHx8dX9K9SooNRqdUZWUOqE4NRanTpTU1Pnq2rfMgfCJa6XH7idhGYAV1fy/JVAm72dZyNX71mdkRWUOlWDU2skR65+Wy0gwHPAYlWdVM457UPnISJH4lY8ZNdclcYYs3e+uqAF9MNtXvGtiHwTeu5GoAuAqj4FjAH+IiL5uH5KY0O/PYwxxjd8Fa6q+glF+2+Wd87jwOM1U5ExxlSPr6YFjDGmtrBwNcaYKLBwNcaYKLBwNcaYKLBwDeOXX2DatESvyzDGBJiFaxiPPQaTJu3PnXd6XYkxJqh8tRTLL+6+G77++lduvrk9eXnw97+DVLhAzBhjSrJwDSM2Fq677ge6dGnPHXdAbi7ce68FrDGm8ixcyxEbC888A3FxcP/9kJcHkyZZwBpjKsfCtQIxMfDEEy5gH3nEBew//uGeN8aYili47oUIPPwwNGhQNIJ9+mkLWGNMxSxcK0HEzbk2aAB33OECdvJkN3VgjDHhWLhWkgjcfjvUrw+33OICNj0d6tlP0BgThkVDFd18sxvBXncd7NoFr7zi5mSNMaY4mzmshokT3Tzsv/8NY8a4pVrGGFOchWs1XXUV/POfMG0anHwy7NjhdUXGGD+xcN0Hf/kLPPsszJgBo0bB9u1eV2SM8QsL1330xz/CCy9AZiYMHw5bt3pdkTHGDyxcI2D8eHj5Zfj0Uxg2DLZs8boiY4zXLFwjZOxYePVV+PJLGDIENm3yuiJjjJcsXCPotNPgzTdhwQIYNAg2bPC6ImOMVyxcI2zUKHjnHfj+exg4ENav97oiY4wXfBWuItJZRDJFZLGILBKRK8OcIyLymIgsE5GFInKEF7VWZNgwePddWLYMUlJg7VqvKzLG1DRfhSuQD1yjqgcCRwOXikjvUucMB5JCHxOAJ2u2xMoZNAjeew9WrYIBA2D1aq8rMsbUJF+Fq6quVdX/hb7eCiwGOpY6bTSQps7nQAsR8WXDqwED4IMPYN069/VPP3ldkTGmQqoReylfhWtxItINOBz4otShjsDPxR6vpmwA+8axx8KHH8LGjdC/Pyxf7nVFxpgStm+HKVNgxAi6Pf98xF5WNIJJHSkiEg/MAe5S1TdLHXsXuEdVPwk9ngVMVNX5YV5nAm7qgISEhOSMjIxK15CTk0N8fHz1/xKlZGXFc+21hxEXV8BDD31Dly6Ru1820rVGi9UZWUGpE3xYa0EBLRYuJGHGDNrOnUu97dvZ2a4dy0eN4rdzz63SS6Wmps5X1b5lDqiqrz6A+sAM4Opyjj8NnFXs8RIgcW+vm5ycrFWRmZlZpfMrY+FC1bZtVdu3V120KHKvG41ao8HqjKyg1Knqo1oXL1a98UbVLl1UQTU+XvX881UzM1V3765WncA8DZM5vpoWEBEBngMWq+qkck6bCowPrRo4GtiiqoG4Hn/IITBnjvs6JQUWLvS0HGPqht9+c/2ZjjwSDjzQ7Xzfu7e7rXLdOnj+efcPMsLtRfy2n2s/YBzwrYh8E3ruRqALgKo+BUwHRgDLgO3ABR7UWW0HHugCduBASE2FmTPhCN8tJjMm4HbuhP/8x+1oP3065OdDnz7w0ENw1lmQGP1r4L4KV3XzqBX2Vw0Nwy+tmYqiY7/9igJ20CC3q9aRR3pdlTEBp+o2+EhPh9deg82boUMHtz/o+PHurWMN8lW41iU9e8LcuS5gBw+G9993KwuMMVW0bBm89JIL1RUroHFjOPVUF6gDB3rW7M5Xc651TdeubgSbmAhDhxbNxxpj9mLjRnjqKTciSUpyDe569IAXX3TzqOnpbgclD7uIWrh6rFMnmD0bunRx+8F++KHXFRnjU3l58PbbboekxES3W/3vv7sLVKtWuQsY48eDT5Z82bSADyQmuoAdPBhGjoS33nJBa0ydp+r28UxPh4wMyM6Gdu3gkktg3Dg4/HDXmtmHLFx9ol07181gyBDXk+v11+Gkk7yuyhiPrFxZNI+6dCk0bOj+YYwb5/6R1K/vdYV7ZeHqI61bw6xZblet005zd+SNGeN1VcbUkC1b4I03IC3NXe0FtynHdde5fxDNm3tbXxXZnKvPtGzpNns58kjX3WDKFK8rMiaK8vPdOtSxY6F9e7joIvj1V7jzTvjxRzdfduGFgQtWsJGrLzVv7pZmjRwJ557r5vHPO8/rqoyJEFX4+ms3Qp0yxe0o37q16/Y5fjz84Q++nUetCgtXn2ra1O0HO3o0XHAB7NrlfqkbE1irV7tbTtPTYdEiiItzI4jzznNzYXFxXlcYURauPta4MUyd6qabLr4YcnPh0kDfm2bqnK1b4c03Oeyxx9xoVdWtTX3qKTjjDDcPVktZuPpco0ZuadYZZ8Bll7kpgr/+1euqjKnA7t3uymx6uuvYuX07DTt0gFtucfNcvXp5XWGNsHANgAYN3NKss8+Gq692AXvddV5XZUwpCxe6QH3lFVizBlq0cGE6fjxf5OWRkprqdYU1ysI1IOLi3Brq8ePh+utdwN58s9dVmTrv119dmKaluZ7y9erBiBFuPerIkW59Krir/nWMhWuA1KvnBgZxce4dVm4u3HGH11WZOmf7dncbanq6WzdYUOCu8D/2mFtS1bat1xX6goVrwMTGwuTJLmDvussF7IgRXldlar2CArezUFqaW+ifk+M2xLj+ejdKPeAAryv0HQvXAIqJcRdb4+LgwQdhxYpepKTUiqWBxm9++MEF6ksvwc8/uzWCZ5zhArV//4jv3l+bWLgGVEyM61wRFwcPP9yJSy6BJ56w/9dNBPz2m1vcn54O8+a5t0tDh8J997mF140be11hIFi4BpiI61qxbt1PPPVUV/Ly4JlnPN3C0gTVzp0wbZoL1Pfec7elHn44TJrk2qK0b+91hYFj4RpwInDRRT/Sq1dXbr/d3ck1ebK7+GVMhVThk0+K2qJs2QIdO7r1fuPGwcEHe11hoNk/wVpABP7+dzdF8H//55ZppacHYlc244Vly9z/IOnpbnOUJk2K2qKkptpbnwixcK1FbrrJBezEiW4EO2VKrbtd21TXxo3w6qsuUD/7zP1GHjTI/VY+5RTf7N5fm1i41jJ/+5u7o+vKK92eBG+84R6bOigvz23nl5bm2kzv2uXe6t93n7vdr1Mnryus1XwXriIyGRgJrFfVMpM+IpICvAP8GHrqTVW9veYq9L8rrnAj1r/8xV3cfestt0eBqQNU4YsvitqibNwICQluY4rx4+Gww2zNXg3xXbgCLwCPA2kVnPOxqo6smXKC6c9/dgF70UXuLsSpU93UmqmlfvyxqC1KVlZRW5Tx411bFLvCWeN89xNX1bki0s3rOmqDCy90AXveea7h4bvvujXgpnaIzcmBZ591gVq8LcoNN7g5oWbNvC2wjvNduFbSMSKyAFgDXKuqi7wuyK/OPdetGjjnHDjhBLeEMYAdM0yhXbtgxgxIT+fYt99286r77+/aopx7LnTt6nWFJkRU1esaygiNXP9TzpxrM6BAVXNEZATwqKomlfM6E4AJAAkJCckZGRmVriEnJ4f4gFxBrUytH3/chttv703Pnjncf/9CmjXLr6HqigTlZ+q7OlWJX7qU9jNn0m7WLOI2byaveXN+Of54Np54Ilv339/386i++5mWozp1pqamzlfVvmUOqKrvPoBuwHeVPHcl0GZv5yUnJ2tVZGZmVul8L1W21mnTVOPiVPv0Uf3tt+jWFE5Qfqa+qXPVKtV77lHt3VsV3H+8005Tfftt1dxc/9RZCUGptTp1AvM0TOYEblpARNoD61RVReRIXAfbbI/LCoTCC1snnwwpKW6z+IQEr6syJYTaopCWBpmZ7ur/ccfB00/D6afX6rYotY3vwlVEpgApQBsRWQ3cCtQHUNWngDHAX0QkH9gBjA399jCVcMIJ7sLWqFFFAduhg9dV1XGFbVHS0lyw7tgBPXvCrbe6edSePb2u0FSD78JVVc/ay/HHcUu1TDUNHOhad48Y4S4uf/QRdO7sdVV10MKFLlBfeQXWrnVtUcaPd/f1H3us7+dRTcV8F66mZhx/vNtEftgwty1nZiZ06+Z1VXXA2rVFbVEWLnTrT088sagtit1OV2vY7p912DHHuHejmze7gF22zOuKaqlt2+Dll91vsk6d4Npr3SL/xx93Yfv2225dqgVrrWIj1zqub183ah08uGiKYP/9va6qFigocE350tLg3/92bVG6dnUL/M8919qi1AEWroY+fVwODBrkAnbWLDjoIK+rCqjvv3d3TL38clFblDPPdG/7jz/eWkXUIfZf2gBus6Q5c9y//ZQU1yXZVNL69a7zad++7rfSAw/AIYe4jVPWrXO3qA4YYMFax9h/bbPHAQe4gG3Y0K0omD/f64p8bMcOt3v/qFFuLduVV7qpgIcfhl9+cevdzjzTtiOrw2xawJSQlOQCduBAN00wYwYcdZTXVflEQQF8+qmbR33tNfj9dxes11xjbVFMGRaupowePYoCdsgQt9lLv35eV+WhrKyitigrV1pbFFMpFq4mrK5d3S52Awe6u7r+8x83F1tnZGcXtUX5/HM3XzpoENxxh2uLYpvjmr2wcDXl6tjRjWAHDXJ3c02d6pZs1Vq5uTB9OgdNmuR28y9si3L//a4tSseOXldoAsTC1VSofXu3TGvwYHcD0ZtvuqCtNVTdyLSwLcqmTTRv2dLaoph9ZuFq9qptW3dzwdChbket1193vbkCrbAtSlqauzWtUSP3lxs3js/i4hgwaJDXFZqAs6VYYXy99ms+WPcBX6z+go07Nnpdji+0bu1uLjj8cBgzxgVs4Gze7Nac9u/vrtrdcovbsWbyZPj1V3fP//DhqF2gMhFgI9cwXv/+de754R7u+eEeAFo3ak1S6yT2a70fSa2KPie1TiI+zv+7q0dKixYwc6abFhg71k1Jnn2211Xtxa5dbguw9HQ3aZyb6xb03nWX631jbVFMlFi4hnFbym0ckHsALXu2ZGn2UrI2ZrE0eymzVswibUHJprTt49uHDd2eLXvSqH7tW0DerJnLqlGj3C3yu3a5Boi+ourugEhPhylT4LffoE0bmDDBrUft29fmUU3UWbiGERcbR5fGXUjZP6XMsW1521i2cRlZG7PIys5i6calZGVnMXXJVH7b/tue8wShc/POZYJ3v9b70a1FN+rH1q/Bv1Fkxce7G5BOPhnOP9/1yLv4Yq+rAlatcvf0p6XBDz+41rcnneQCdfhw16nRmBpi4VpFTeKacFj7wzis/WFljm3ZuWXPKDcrO2vP11O+m8LmnZv3nFcvph7dW3QnqXWSG+kWhm/rJDo360xsjP/n/Bo3du+yTz3VDQjz8uDSSz0oZOtWt+tUWppb1mBtUYxPWLhGUPOGzenboS99O5RsBKmqbNi+oWi0G5pqyNqYxZyVc9i2a9uecxvENqBnq55lRrtJrZNIjE9EfPR2tmFDeOstdwv9ZZe56cyrr66Bb5yfDx9+6N72v/VWUVuU225zcxU9etRAEcZUzMK1BogIbZu0pW2Tthzb+dgSx1SVtTlr94x2iwfve8veI2933p5zm9Rvsme0Wzx4t+zagqp6ErwNGriVA2ef7W6xz811W5ZGxYIFRW1Rfv3VjUrPO8+tRz36aJtHNb5i4eoxEaFD0w50aNqBlG4pJY7tLtjNz7//XGaa4etfv+bNxW+yW3fvObfF/1oUXVArNtpNapVE84bNo/p3qF/fXTeKi4Mbb3RTBLfcEqGsW7PGhWl6umuLUr++a4syfrxbtmC79xufsnD1sdiYWLq16Ea3Ft0Y2nNoiWN5u/NYuXklWdlZTP9yOtpSydqYxcerPuaVb19BKWqI265Ju7ArGnq16kXj+o0jUmu9em5QGRfn3p3n5rrVTtUK2G3b3Nv99HT39r+gwG3N9cQTbg6ideuI1GxMNFm4BlRcbBz7td6P/VrvR5M1TUgptqvKjl07WL5pedhphue/eb7E63Rq1ins/G6Plj2Ii42rUk2xsfDcc24wec89bgT7wAOVDNjdu90FqfR0eOMNF7Bdu7qh8LhxsN9+VarFGK/5LlxFZDIwElivqmU2yBQ3sfgoMALYDpyvqv+r2Sr9rVH9Rhzc7mAObld2f9GtuVv3XFgrDN2l2Ut54/s3yN6Rvee8GImhW4tuYacaujbvWu6KhpgYePJJN4J96CEXsI8+WkHALlpU1BZl9Wq3kHbsWPe2/7jjbPd+E1i+C1fgBeBxIK2c48OBpNDHUcCToc+mEpo2aMoRiUdwROIRZY5t3LGxzGh3afZSPl31KVvztu45Ly42jh4te4QN3g5NOxAjMTz6aFHALlr0MsuW3cTPP6+iS5cu3HXddZyTl+fmEf73PzfkHTbMnTxqlO3eb2oF34Wrqs4VkW4VnDIaSFNVBT4XkRYikqiqa2ukwFqsVaNWHNXpKI7qVPJ3laqybtu6sCsaPlj+ATvzd+45t1G9RkXrd4ftR59lG/hoahpoLgA//fQTEy65BIBzjjjCtUU5+2xo167m/qLG1ADfhWsldAR+LvZ4deg5C9coERHax7enfXx7+nftX+JYgRaw+vfVJVY0ZG3M4tv13/LOknfIn51PsWtrgJvLubRNPPr8X92It2l9bKm/qW3EDQD9JTRy/U85c67vAveo6iehx7OAiapapp2eiEwAJgAkJCQkZ2RkVLqGnJwc4uODsSmL32pt9PPPJMycSesPP6D52nWls7XIbUVfNq/fnE6NOhV9NHafOzbqSKPYmp0m8NvPszxBqROCU2t16kxNTZ2vqn1LPx/EketqoHOxx52ANeFOVNVngGcA+vbtqylV6FMye/ZsqnK+l3xRa3a222w6Pd3t4h8TA4MH0yVvPj9lZ5c5vXPnLrx/yXtlphm+zf6WGetmlDi3Q9MOYVc09GzZkwb1Ir/O1Rc/z0oISp0QnFojWWcQw3UqcJmIZOAuZG2x+VaP5Oa6HVzS093nXbvgkEPc+quzz4YOHbjr5ZeZMGEC27dvL/YHG9Oq1d30aNqb3m17l3nZbXnbSqxoKAzft394u8TmODESQ5fmXcpcVCvcHKdeTBD/9za1he/+7xORKUAK0EZEVgO3AvUBVPUpYDpuGdYy3PTdBd5UWkeFaYtC+/Zw+eVuPWqfPiVOP+eccwC46aabWLXKrRYYNOguJk8+h5NOgrffdpvAFNckrgl92vehT/uSrwWwaccmlm1cxtLspSVGvC9/+zJbcrfsOa9wc5zSo939Wu9Hp2adiBFb4mWiy3fhqqpn7eW4Al7sv1S3rVjh2qKkpxe1RTnlFBeogwe7W7TKcc4553DOOeeUeMt13HHwxz+6O1mnTXPbGFZGy0Yt+UPHP/CHjn8o8Xzh5jjF998t/PzRjx+xI3/HnnMb1mtIr1a99ox4i2+E3j6+fZV/NMaE47twNT6yaZPblSU9HT75xN0JkJICN90Ep50GTZtW+6UvuMCtgx0/3m21+u677v6B6iq+OU6/Lv1KHCvQAtZsXVNifndp9lIWb1jMu1nvltgcJz4unsS4RI7YcESZEW+rRq2qX6CpcyxcTUmFbVHS0tyQMjcXDjwQ7r7btUXp0iVi3+qcc1zAnnWWa374/vuulUykxUgMnZp1olOzTqR2Ty1xbHfBblZtWVViH97Psz7nqzVf8fr3r1OgBXvObdWoVdgLa0mtkmjaoPq/aEztZOFqitqipKW57a02bHAtX//0Jze0POKIqG3nd/rpbqOrM85wswsffACtanCAGBsTS/eW3enesvuezXFmN3LTF3m781ixaUWZC2uZKzNJX5he4nUSmiSEnd+tre1+zN5ZuNZlhW1R0tNh8WK348pJJ7lAPeGEGmuLcvLJbhOs006DgQNdE8S2bWvkW1coLjaOA9ocwAFtDihzbPuu7a7dT6ngfTfrXSZ/M3nPeYXtfsKNeLu36B7odj+mYhaudc3vv7u2KOnpkJnpnjvuOHjmGTeMjMb78ko48UTXNmb0aEhNdTsNtvfxtaXG9RtzaMKhHJpwaJljv+f+HnaPhtLtfmLFjZrDXVjr0rxLINr9mPJZuNYF+fluOJie7tY+7dgBvXrB7be7tijdu3tdIeDmXd97D0aOdNfNZs2Cjh29rqrqmjVoRnKHZJI7JJd4XlXJ3pFdJnizsrOY+9PcEu1+4mLj6NmyZ9j53Q5NO/iq3Y8Jz8K1tlIt2RZl3To3mXn++e5t/1FH+bItSkqKu7A1fDgMGAAffRTRa2ieEhHaNG5Dm8ZtOKbzMSWOqSq/5vxaFLrFOgu/v+x9cnfn7jm3Sf0m9GrVixYFLei3u1+JEW+bxm0seH3CwrW2WbOmaB71228D2RbluOPcQHvYsKKA9cngOmpEhMSmiSQ2TWRAtwEljhW2+yk94l2wegGffPpJyXY/DVsUTTG02q9Ez7Vot/sxJVm41gIxO3a4Bf5pae69dEGBa9gX4LYoRx/t/ipDhkD//m56uFcvr6vyRvF2P0N6Dtnz/OzZs+l3fD/X7qdUS/dPV33KlG+nlGj307Zx26KVDK2KVjREst2PKWLhGlSFbVHS0uj32muwcyd06+YW+I8bB0lJXle4z5KTXagOHuwC9qOP4ICyF+7rtPqx9d3otHUSI5JGlDi2M38nyzcuL3PX2oxlM3gh54US53Zs2rFM6Ca1cu1+orE5Tl1g4Ro0ixa5EerLL8Mvv0Dz5qwbNIgOEyfWyrYohx3mfocMGuSmCGbNgoPLdq8xYTSs15CD2h3EQe0OKnNsa+5Wt5SsWOhmZWfx5g9vsmH7hj3nxUgMXZt3DbuGt6J2P8bCNRjWrXOL+9PS4OuvXVuU4cNh0iQ46SSWfv45Hfr33/vrBNRBB7mAHTjQLdOaObPM/jCmipo2aMrhiYdzeOLhZY5t2rGpzGg3KzuL//783xLtfurH1KdHyx4l5nULP3ds1rHOb45j4epXO3a4hZ9paTBjhpsGSE6GRx5x94vWsbYoBxwAc+e6gB040N3J1bfM9sQmElo2alluu5/129aHDd4PV3xYpt1Pr1a99kw17P5tN/VW1SOpVRLtmrSrEysaLFz9pKAAPv7YBeobb7gF/506wd/+5uZRe5fd+7Qu6dUL5sxx4TpokPudc/TRXldVd4gICfEJJMQncHzX40scK9ACfvn9lzIX1hatX8S0JdPYVbCLB5Y+ALh1wMVHunt6rrVOqlWb41i4+sGSJW7p1EsvwU8/uf33xoxxgZqSUuvmUfdF9+5FATtkCEyfDscfv/c/Z6IrRmLo3LwznZt3ZmD3gSWO5Rfk8+qMV2md1LrEPryfrf6MVxe9WmJznNaNWpe4oFY8gOPj/N8mpjgLV69s2ACvvupGqV9+6QJ0yBC46y53s32TJl5X6FtduhRNEQwbBv/5j5uLNf5UL6YeHRt1JKVXCsN6DStxLDc/122OU3jjRPZSlm5cyqwVs0hbkFbi3MT4xLAX1og6kLEAABlRSURBVHq07EHDeg1r8q9UKRauNSk31yVBYVuU/Hw49FB48EHXFiUx0esKA6NDBzeCHTTI3Rvxzjvu9lkTLA3qNeDAtgdyYNsDyxzblreN5ZuWlxjtLs1eyjtL3inR7kcQ1+4nzIW1bi26ebY5joVrtKnCZ5+5QH311aK2KFdd5d72H1p24w9TOQkJbh3skCEwahS8+aa7Gc3UDk3impS7Oc7mnZuLWrkX25nslW9fCdvuJ9wa3s7NO0d1RYOFa7QsX17UFmX58qK2KOPHu+FWBW1RTOW1betuLhg61P14X33VfTa1W4uGLSps91N6miErO4vZK2ezfVdRo8yG9RrSs2XPEsGrW5UUUiJSo/0Lj6RNm+C111ygfvqp2xglNRVuvhlOPXWf2qKY8rVq5bYoHD7c7Zr4yitu821T9xRv93Ns52NLHFNV1+6n2IqGpRuXsmTDEqZnTSdvdx6D2w3mYi6OSC0WrvsqL69kW5S8PNcW5d573Txq585eV1gntGjh1r6OGOGWAe/a5drIGFNIROjYrCMdm3UkpVtKiWOF7X4++/yziH0/C9fqUIV584raomRnu/enf/5z1NuimPI1bep+z40a5aaz8/JcI0Rj9qaw3c9PjX6K2Gv6LlxFZBjwKBALPKuq95Y6fj7wAPBL6KnHVfXZGinup5/cPf1paW5taoMGbuv88ePdpF8NtUUx5WvSxC3IOOUUuPBCF7B/+pPXVZm6yFfhKiKxwBPAEGA18JWITFXV70ud+qqqXlYjRRW2RUlLcze4g1u1fu21bqG/R21RTPkaN3ZLs8aMcW8m8vLg8su9rsrUNb4KV+BIYJmqrgAQkQxgNFA6XKMrP59WX3wBTz/t2qLs3Om28PNZWxRTvoYN3dKsM8+EK65wAXvNNV5XZeoSv4VrR+DnYo9XA0eFOe80EekPLAX+qqo/hzmn+m64gUMffNBdhr7wQjeB59O2KKZ8cXFu8ca557o3Grm5cOONXldl6gpR1b2fVUNE5HTgBFW9KPR4HHCkql5e7JzWQI6q5orIn4EzVHVgOa83AZgAkJCQkJyRkVGpOhqvXIlkZbE9JQUNwDxqTk4O8fH+v+/aqzp37xbuu29/Zs5sz/jxKzn//JUV/p60n2fkBaXW6tSZmpo6X1XL7tGmqr75AI4BZhR7fANwQwXnxwJbKvPaycnJWhWZmZlVOt9LQanVyzrz81UvuEAVVK+/XrWgoPxz7ecZeUGptTp1AvM0TOb4bVrgKyBJRLrjVgOMBc4ufoKIJKrq2tDDk4DFNVuiCaLYWHj2WTdVcO+9borgoYdspsdEj6/CVVXzReQyYAZuVDpZVReJyO243w5TgStE5CQgH9gInO9ZwSZQYmLgySddwD78sLvI9dhjtqOjiQ5fhSuAqk4Hppd67pZiX9+Amy4wpspE4NFH3RLlBx90AfvUUxawJvJ8F67GRJsI3H+/G8HefbcL2Oeec1MHxkTKXsNVRA5Q1R9qohhjaooI3HmnC9jbbnN7Ebz4om1WZiKnMv8rfS0izwC3qeqmaBdkTE0RgVtvdVMEN9zgRrCvvOJ1Vaa2qMxM05HAQUCWiFweukXVmFrj+utdl/I33nC3zObl2RICs+/2OnJV1W+BwSJyMm7DlL+IyDWq+l7UqzOmhvz1r26K4LLLYN26g+nf391Ca0x1Vfoaqaq+jRvBvghkiMh0ETkgapUZU8MuvRSeeQa+/LIVo0bB9u17/zPGlKeqC1AaA/NxAXsCsFBEHhOR5hGvzBgPXHwxTJz4Ax995Ppx5eR4XZEJqr2Gq4hcJSIvi8hSIBuYBvwBt+fqRcD+wPciEm6DFWMCZ9iwdaSnw8cfu9bdv//udUUmiCqzWuAa4DPgSeBzYL6q5hU7niYi1wGTcdMGxgTe2We7OdizznLdZd9/H1q29LoqEySVuaBVmSZQzwN373s5xvjHmDEuYMeMgcGDXY+u1q29rsoERaRu+vsNCLvtnzFBdtJJrqvBokUwcCCsX+91RSYoIhKuoZ235kTitYzxm+HDXV+urCzXKX3t2r3/GWNsuwpjKmHwYHjvPdejMiUFfvllr3/E1HEWrsZU0oABMGOGG7n27++C1pjyWLgaUwX9+sGHH0J2tgvbFSu8rsj4lYWrMVV05JHw0UewdasL2KwsrysyfmThakw1HHEEZGa6djH9+8NiazZkSrFwNaaaDj0UZs8GVTeC/fZbrysyfmLhasw+6N0b5syB+vXdMq2vv/a6IuMXFq7G7KP994e5c6FJE3ejwVdfeV2R8QMLV2MioGdPN4Jt2dKtif3sM68rMl6zcDUmQrp1cwHbrh0MHepGs6bu8l24isgwEVkiIstE5PowxxuIyKuh41+ISLear9KY8Dp3dqHaqZO7bXbWLK8rMl7xVbiG+nM9AQwHegNniUjvUqf9Edikqr2Ah4H7arZKYyqWmOhWEfToASNHuu0KTd3jq3DFNUNcpqorQnvGZgCjS50zGtcJAeANYJCIWEc54ysJCW4d7AEHwOjRMG2a1xWZmiaq6nUNe4jIGGCYql4UejwOOEpVLyt2znehc1aHHi8PnbMhzOtNACYAJCQkJGdkZFS6lpycHOLj4/flr1NjglJrXazz99/rMXHioSxbFs8tt3xP//5l/jettqD8PCE4tVanztTU1Pmq2rfMAVX1zQdwOvBsscfjgH+UOmcR0KnY4+VA6729dnJyslZFZmZmlc73UlBqrat1bt6seswxqrGxqhkZkXvdoPw8VYNTa3XqBOZpmMzx27TAaqB454NOwJryzhGRekBzYGONVGdMNTRv7nbT6tfPtY9JT/e6IlMT/BauXwFJItJdROKAscDUUudMBc4LfT0G+Cj028MY32raFKZPd3vBnncePPec1xWZaPNVuKpqPnAZMANYDLymqotE5HYROSl02nNAaxFZBlwNlFmuZYwfNWniOhoMHQoXXQRPPul1RSaaKtP9tUap6nRgeqnnbin29U7c3KwxgdOoEbz9Npx+OlxyCeTlwZVXel2ViQZfjVyNqQsaNoR//xtOPRWuugoeeMDrikw0WLga44G4OMjIgLFjYeJEuPNOrysykea7aQFj6or69eGll9znm292UwR//zvYLTG1g4WrMR6KjYXnn3cj2TvucJ0N7r3XArY2sHA1xmOxsfDMMy5g77/fjWAnTbKADToLV2N8ICYGnnjCBewjj7iA/cc/3PMmmCxcjfEJEXj4YRewDzzgAvbppy1gg8rC1RgfEYH77oMGDdwKgrw8mDzZTR2YYLFwNcZnRNzFrbg4uOUW2LUL0tKgnv1rDRT7z2WMT918swvY6693I9hXXnGPTTDYbI4xPnbddW4e9t//hjFj3FItEwwWrsb43FVXuZUE06bBySfDjh1eV2Qqw8LVmAC45BL417/cvrCjRsG2bV5XZPbGwtWYgLjoInjhBdeba8QI2L7dlhD4mYWrMQEyfjy8/DJ8+ilMnHgoW7Z4XZEpj4WrMQEzdiy89hr88ENThgyBTZu8rsiEY+FqTACdeircfvsiFiyAQYNgQ+SaypoIsXA1JqCOPTabd96BxYshNRXWr/e6IlOchasxATZsmOvLtWKFa364dq3XFZlCFq7GBNygQfDee7BqFQwYAD//7HVFBixcjakV+veHDz6AdetcwK5c6XVFxsLVmFri2GPhww/d6oEBA2D5cq8rqtt8E64i0kpEZopIVuhzy3LO2y0i34Q+ptZ0ncb42R/+AB995O7g6t8flizxuqK6yzfhClwPzFLVJGBW6HE4O1S1T+jjpJorz5hgOPxwdxdXfr67yPX9915XVDf5KVxHAy+Gvn4RONnDWowJtEMOgdmz3dcpKbBwoZfV1E2iql7XAICIbFbVFsUeb1LVMlMDIpIPfAPkA/eq6tsVvOYEYAJAQkJCckZGRqXrycnJIT4+vgp/A+8EpVarM7IqU+fq1Y3461/7kJcXwwMPLGC//XJqqLqSatPPtLTU1NT5qtq3zAFVrbEP4EPguzAfo4HNpc7dVM5rdAh97gGsBHpW5nsnJydrVWRmZlbpfC8FpVarM7IqW+fy5apduqi2aKH6xRfRrak8te1nWhwwT8NkTo12IlDVweUdE5F1IpKoqmtFJBEIe7+Jqq4JfV4hIrOBwwG7LmpMOXr0gLlzYeBAGDwY3n/frSww0eWnOdepwHmhr88D3il9goi0FJEGoa/bAP0Am643Zi+6doU5cyAxEYYOdV+b6PJTuN4LDBGRLGBI6DEi0ldEng2dcyAwT0QWAJm4OVcLV2MqoVMnd5GrSxcYPtytiTXR45sGhaqaDQwK8/w84KLQ1/8FDqnh0oypNRITXcAOHgwjR8Jbb7mgNZHnp5GrMaYGtGvn1sH27u16ck21W3GiwsLVmDqodWuYNQv69IHTTnPdZU1kWbgaU0e1bAkzZ8KRR8KZZ8KUKV5XVLtYuBpThzVr5jrKHnccnHsuvPji3v+MqRwLV2PquPh4mD7drYO94AJ49tm9/xmzdxauxhgaN3YXtk44AS6+GP75T68rCj4LV2MMAI0awdtvw0knwaWXwsMPe11RsFm4GmP2aNAAXn/drSC4+mq47z6vKwouC1djTAlxcZCRAWedBddfD7ff7nVFweSbO7SMMf5Rrx6kp7ugvfVWyMuDO+4AEa8rCw4LV2NMWLGxMHky1K8Pd90Fublw//0WsJVl4WqMKVdMDDz9tBvBPvigG8E+8ogFbGVYuBpjKhQTA48/7gL2kUdg1y73OMau2FTIwtUYs1ciMGmSW01w331uBPv0027qwIRn4WqMqRQRuOceN4K94w4XsJMnu4tfpiz7sRhjKk3ELc2Ki4Obb3YBm57uLnqZkixcjTFV9n//56YIJk50c7BTprjANUVsStoYUy1/+xs8+ii8+SaMGeOWapkiFq7GmGq74gp48kmYNg1Gj4YdO7yuyD8sXI0x++TPf4bnnoMPPnB9ubZt87oif7BwNcbsswsvdBttz57tGh5u3ep1Rd7zTbiKyOkiskhECkSkbwXnDRORJSKyTESur8kajTHlGzcOXnkF/vtfty/sli1eV+Qt34Qr8B1wKjC3vBNEJBZ4AhgO9AbOEpHeNVOeMWZvzjwTXnsN5s1z7bs3bvS6Iu/4JlxVdbGqLtnLaUcCy1R1harmARnA6OhXZ4yprFNPdSsIFi6EQYNgwwavK/KGb8K1kjoCPxd7vDr0nDHGR0aOdG1jfvgBUlJg48a6d5dBjd5EICIfAu3DHLpJVd+pzEuEeU4r+H4TgAkACQkJzJ49uzJlApCTk1Ol870UlFqtzsjye50NGsBdd7XgppsO4YorDqOg4L+0aZPndVkViujPVFV99QHMBvqWc+wYYEaxxzcAN1TmdZOTk7UqMjMzq3S+l4JSq9UZWUGpc+5c1UaNdmmvXqqrVnldTcWq8zMF5mmYzAnatMBXQJKIdBeROGAsMNXjmowxFTj+eHjwwYWsXw8DBsDKlV5XVDN8E64icoqIrMaNTt8VkRmh5zuIyHQAVc0HLgNmAIuB11R1kVc1G2Mqp3fv35k1CzZvhv79YdkyryuKPt+Eq6q+paqdVLWBqiao6gmh59eo6ohi501X1f1Utaeq3uVdxcaYqujbFz76yN0iO2AALNnb2qCA8024GmNqvz59IDMT8vNdwC6qxe87LVyNMTXq4INhzhzXJiYlBRYs8Lqi6LBwNcbUuAMOgLlzoVEjSE2F+fO9rijyLFyNMZ7o1cuNYJs1c3dyffGF1xVFloWrMcYz3bu7EWybNjBkCHzyidcVRY6FqzHGU126uBFsYiIMG+YueNUGFq7GGM917OgCtmtXOPFEmDnT64r2nYWrMcYX2rd3m20nJcGoUTB9utcV7RsLV2OMb7Rt6240OOggOPlkeKcy2zn5lIWrMcZXWreGWbPgiCNcV9nXX/e6ouqxcDXG+E6LFq7h4dFHw9ixrn1M0Fi4GmN8qVkzeO89t9HLuefCCy94XVHVWLgaY3wrPh7efdf147rgAnjmGa8rqjwLV2OMrzVu7FrGjBgBf/oTPP641xVVjoWrMcb3GjZ0TQ9Hj4bLL4dJk7yuaO8sXI0xgdCggVs5cPrpcM01cO+9XldUsRptUGiMMfuifn23cqB+fbjhBsjLg5tvBgnXutRjFq7GmECpVw/S0iAuDm69FXJz4c47/RewFq7GmMCJjYXnnnMBe/fdLmAfeMBfAWvhaowJpJgYeOopF7APPeSmCB591D8Ba+FqjAksEXjsMXexqzBg//lPF7xes3A1xgSaiJsSaNDATRHk5cG//uWmDrzkm3AVkdOB24ADgSNVdV45560EtgK7gXxV7VtTNRpj/EnEXdRq0MBd5MrLc7fL1vMw4XwTrsB3wKnA05U4N1VVN0S5HmNMgIjALbe4ZVo33gi7dsFLL7nHXvBNuKrqYgDxy2y0MSaQbrjBjWCvucYFbEaGu+hV03ww7VtlCnwgIvNFZILXxRhj/Ofqq+Ef/4C33oJTT4WdO2u+BlHVmvtmIh8C7cMcuklV3wmdMxu4toI51w6qukZE2gEzgctVdW45504AJgAkJCQkZ2RkVLrWnJwc4uPjK32+l4JSq9UZWUGpE7yrddq0RCZN2p++fTdyxx3f0bBhQYXnV6fO1NTU+WGv/aiqrz6A2UDfSp57Gy6I93pucnKyVkVmZmaVzvdSUGq1OiMrKHWqelvr5MmqIqqpqao5ORWfW506gXkaJnMCNS0gIk1EpGnh18BQ3IUwY4wJ64IL3O2yc+a41t1bt9bM9/VNuIrIKSKyGjgGeFdEZoSe7yAihX0gE4BPRGQB8CXwrqq+703FxpigOPdcmDIFPvsMhg6FzZuj/z39tFrgLeCtMM+vAUaEvl4BHFbDpRljaoEzznCrBs44w3U2+OADaNUqet/PNyNXY4yJtpNPdisIvvsOBg6E336L3veycDXG1CknnujaxixZAqmp8Ouv0fk+Fq7GmDpn6FDX+PDHHyElBX75JfLfw8LVGFMnDRwI77/vgnXAAFi1KrKvb+FqjKmzjj8eZs6EDRtcwK5d2zBir23haoyp044+GmbNgi1b4LHHkiL2ur5ZimWMMV5JTnY3GSxf/gPQLyKvaSNXY4wBDjkEWrTYFbHXs3A1xpgosHA1xpgosHA1xpgosHA1xpgosHA1xpgosHA1xpgosHA1xpgosHA1xpgosHA1xpgosHA1xpgoqNHW2l4Skd+An6rwR9oAG6JUTqQFpVarM7KCUicEp9bq1NlVVduWfrLOhGtVicg8DdeL3IeCUqvVGVlBqROCU2sk67RpAWOMiQILV2OMiQIL1/I943UBVRCUWq3OyApKnRCcWiNWp825GmNMFNjI1RhjosDCtQIicrqILBKRAhHx3ZVOERkmIktEZJmIXO91PeURkckisl5EvvO6loqISGcRyRSRxaH/7ld6XVM4ItJQRL4UkQWhOv/udU0VEZFYEflaRP7jdS0VEZGVIvKtiHwjIvP29fUsXCv2HXAqMNfrQkoTkVjgCWA40Bs4S0R6e1tVuV4AhnldRCXkA9eo6oHA0cClPv2Z5gIDVfUwoA8wTESO9rimilwJLPa6iEpKVdU+kViOZeFaAVVdrKpLvK6jHEcCy1R1harmARnAaI9rCktV5wIbva5jb1R1rar+L/T1VlwgdPS2qrLUyQk9rB/68OXFExHpBJwIPOt1LTXNwjW4OgI/F3u8Gh8GQVCJSDfgcOALbysJL/RW+xtgPTBTVX1ZJ/AIMBEo8LqQSlDgAxGZLyIT9vXF6nxrbRH5EGgf5tBNqvpOTddTBRLmOV+OXoJGROKBfwNXqervXtcTjqruBvqISAvgLRE5WFV9NactIiOB9ao6X0RSvK6nEvqp6hoRaQfMFJEfQu+6qqXOh6uqDva6hmpaDXQu9rgTsMajWmoNEamPC9aXVfVNr+vZG1XdLCKzcXPavgpXoB9wkoiMABoCzUTkJVU91+O6wlLVNaHP60XkLdzUW7XD1aYFgusrIElEuotIHDAWmOpxTYEmIgI8ByxW1Ule11MeEWkbGrEiIo2AwcAP3lZVlqreoKqdVLUb7v/Pj/warCLSRESaFn4NDGUff1lZuFZARE4RkdXAMcC7IjLD65oKqWo+cBkwA3fh5TVVXeRtVeGJyBTgM2B/EVktIn/0uqZy9APGAQNDy3G+CY26/CYRyBSRhbhfsjNV1dfLnAIgAfhERBYAXwLvqur7+/KCdoeWMcZEgY1cjTEmCixcjTEmCixcjTEmCixcjTEmCixcjTEmCixcjTEmCixcjTEmCixcjTEmCixcjWHPxui5ItK12HOPishyEUnwsjYTTHaHljHs2VfgK+BrVb1YRK7FbZXXT1WzvK3OBFGd3xXLGHAbUIvIjbg9JJYDN+F2+7dgNdViI1djihGR/+K2mhulqu95XY8JLptzNSZERAYCh+E2Il/ncTkm4GzkagwgIocBc4CrcT2f4lX1BG+rMkFm4WrqvNAKgf8CT6vq7SJyMLAQN+c629PiTGBZuJo6TURaAZ8Cc1X1T8WefxXooqrHeFacCTQLV2OMiQK7oGWMMVFg4WqMMVFg4WqMMVFg4WqMMVFg4WqMMVFg4WqMMVFg4WqMMVFg4WqMMVFg4WqMMVHw/2PRmhu1osmkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[2, 3], [1, -4], [1, 10]])\n",
    "\n",
    "\n",
    "# take b to be the sum of the two columns of A\n",
    "b = A[:,0] + A[:,1]\n",
    "\n",
    "# Form the matrix A.T @ A\n",
    "ATA = A.T @ A  \n",
    "\n",
    "# Form the RHS vector:\n",
    "rhs = A.T @ b\n",
    "\n",
    "# plot the three lines\n",
    "x = np.linspace(-1,5,100)\n",
    "\n",
    "y1 = -(2./3.)*x + (b[0]/3.)\n",
    "y2 = (1./4.)*x - (b[1]/4.)\n",
    "y3 = -(1./10.)*x + (b[2]/10.)\n",
    "\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "ax1.set_xlabel(\"$x$\", fontsize=14)\n",
    "ax1.set_ylabel(\"$y$\", fontsize=14)\n",
    "ax1.set_title('Two lines', fontsize=14)\n",
    "ax1.grid(True)\n",
    "\n",
    "ax1.plot(x,y1,'b')\n",
    "ax1.plot(x,y2,'r')\n",
    "ax1.plot(x,y3,'g')\n",
    "\n",
    "# plot what we hope should be the solution\n",
    "ax1.plot(1, 1, 'ko')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - Least squares solution as a compromise between all constraints\n",
    "\n",
    "At the end of the lecture we showed an example and noted that the least squares solution was attempting to satisfy all three constraint equations, and that the specific value found was the one which minimised $\\| A\\boldsymbol{x} - \\boldsymbol{b}\\|_2$.\n",
    "\n",
    "By perturbing the values of the obtained least squares solution, show that it is indeed the case that these lead to $\\| A\\boldsymbol{x} - \\boldsymbol{b}\\|_2$ growing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2  1  1]\n",
      " [ 3 -4 10]] [ 7  3 -1]\n",
      "array([16, -1])\n",
      "[[ 2  1]\n",
      " [ 3 -4]] [7 3]\n",
      "array([17,  9])\n",
      "[[ 2  1]\n",
      " [ 3 10]] [ 7 -1]\n",
      "array([13, 11])\n",
      "[[ 1  1]\n",
      " [-4 10]] [ 3 -1]\n",
      "array([  2, -22])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x24e10f15e48>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAFFCAYAAABYL66IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeVhV1frHP4tZZFBQQQXFGCRTroqWNImSlkNOabfC0rxeNc30mpVD2XCD6pdZVna7di2HcMzyZmqTgkN6S8sxTXEC5wEcIBkE1u+PxXTggIDncAbW53n2o2fvdfZ+F8OXddZ61/sVUko0Go1GY1ocLB2ARqPR2CNaXDUajcYMaHHVaDQaM6DFVaPRaMyAFleNRqMxA1pcNRqNxgxYlbgKIQKFEIlCiANCiN+FEBOMtIkWQlwRQuwqPGZYIlaNRqOpDCdLB1CGPOBZKeVvQghP4FchxA9Syv1l2m2WUva1QHwajUZTJaxq5CqlPCOl/K3w/xnAAaC5ZaPSaDSa6mNV4loaIUQQ0AH42cjlKCHEbiHEOiHEbbUamEaj0VQBa5sWAEAI4QGsBCZKKa+Wufwb0FJKmSmE6A2sAkKN3GMUMAqgXr16kYGBgdWKoaCgAAcHq/3bUy10X6wPe+kH1O2+HDp06KKUsrHRi1JKqzoAZ+A7YFIV2x8HGlXWJjIyUlaXxMTEar/HWtF9sT7spR9S1u2+ADtkBbpjVX9uhBACmAcckFLOqqCNf2E7hBC3o6Y20movSo1Go7kx1jYtcBfwOLBXCLGr8Nw0oAWAlPJjYDDwlBAiD8gCHin8C6LRaDRWg1WJq5RyCyBu0OZD4MPaiUij0WhqhlVNC2g0Go29oMVVo9FozIBVTQtoNJVx9epVvLy8OHDggKVDuWm8vb3toh9g332pX78+AQEBNUo10+KqsQmuXr3KuXPnaNGiBb6+vhQmjNgsGRkZeHp6WjoMk2CvfSkoKODUqVNcvHiRJk2aVPteelpAYxOcP3+e5s2b4+rqavPCqrENHBwc8PPz48qVKzV7v4nj0WjMwvXr16lXr56lw9DUMZydncnLy6vRe7W4GmHfPkhMNL6jTWM59IhVU9vczM+cnnM1wj//CStWtOGWW2DYMEtHo9FobBE9cjXCZ59Bx46XePJJ+OQTS0ej0WhsES2uRnB3h/j4fTzwAIwaBXPmWDoijUZja2hxrQAXlwK++gr69YOnn4Z337V0RBpb5Y033qBz5854eXnRuHFjHnzwQfbvL2uuYX7mzJlDREQEXl5eeHl5ERUVxZo1a2o9jqCgIIQQ5Y4+ffrUeizmRItrJbi6wooV8NBDMGkSvPWWpSPS2CJJSUmMHTuWrVu3smHDBpycnOjXrx/p6ek3fe/hw4fzyiuvVKltQEAAb731Fr/99hs7duyge/fuDBgwgD179tx0HNVh+/btnDlzpvj47bffEELw8MMP12oc5kaL6w1wcYGlS+HRR2HKFLXYpdFUh++++44nn3yStm3b0q5dOxYtWsTFixf56aefAFixYgWurq6kpKQUv2fChAkEBwdz7tw5k8XRv39/evXqRUhICGFhYcTFxeHp6cm2bduK2wQEBDBrlmG1z7179+Lm5may0Xbjxo3x9/cvPtauXYuXlxdDhgwxyf2tBS2uVcDJCRYtgieegBkz4KWXQBc51NSUjIwMCgoKaNiwIQCDBw+mXbt2vP766wDMnDmTJUuW8O233+Ln52eWGPLz81m6dCmZmZnceeedxeejoqLYvn27QduJEycycuRI2rRpY3A+Pj4eDw8PmjZtioeHh9Fj8+bNlcYhpWTevHkMHToUd3d303XQCtCpWFXE0VFlEbi4wOuvQ06OmibQqZeWY+JE2LXrxu1MSfv28N57N3ePCRMmEBERQVRUFKByKePj4+nTpw/BwcHExcWxYcMGQkPLuRfdNHv37iUqKors7Gw8PDz46quvaNeuXfH1qKgoPvroo+LXq1atYufOnSxfvrzcvcaMGcPDDz9MZmYmHh4eRp/XvHnl/qI//PADx44dY+TIkTXskfWixbUaODjAv/+tBPbtt5XAvveeFlhN1Zk0aRJbtmzh22+/xdHRsfh8z5496dy5My+++CKrV6+mc+fORt8fHx9PfHx88eucnByEEMycObP43Lp167jnnnuMvr9169bs2rWLy5cvs3LlSoYNG0ZSUhJt27YFoEuXLjz77LOkp6dTv359Jk+ezIwZM/D19S13Lx8fH3x8fG6qtsAnn3xC586dad++fY3eb81oca0mDg7w4YdqsevddyE3V6Vq2Yk/m01xsyPI2uYf//gHS5cuJTExsdyIbsOGDezevRspZaVTAUWjxSJeeOEFmjdvzjPPPFN8rrLRoouLCyEhIQB06tSJ7du38+677zJv3jwAIiMjcXFxYceOHezcuRMnJyfGjRtn9F5lhd4YlQn9+fPn+e9//8scO8111OJaA4SAd95RI9i33lICO3eumjrQaIwxYcIEli5dSlJSEuHh4WRkZBRf2717N4MGDeKDDz5gzZo1TJ06le+++87ofYpGi0V4enri4+NTLJjVpaCggJycnOLXrq6udOjQgdWrV7NgwQIWL16Ms7Oz0ffe7LTAZ599hqurK4888kiNYrd2tLjWECHgjTfUCPa115TAfvaZWvzSaEozbtw4Fi1axKpVq2jYsCFnz54lMzMTIQRpaWn07t2bSZMmMWLECG6//XYiIiJISkoiOjrapHFMmTKFPn36EBgYSEZGBosXLyYpKalcrmtUVBSzZ8+mR48e9O3bt8L73cy0gJSS//znPzzyyCN2U66wLFoKbgIh4NVX1Qj2xReVwH7+OVTwh15TRylaIIqJiTE4P378eH744Qf69u3LjBkzAGjbti1Dhgxh6tSpBilSpuDs2bMMHTqUs2fP4u3tTUREBOvWreP+++83aNe+fXscHBzKpWSZkqSkJA4fPkxCQoLZnmFptLiagOnT1Qj2uefg+nWVF+viYumoNNaCMXPiykZ7y5Ytq/K958+fb/K2CQkJjB49mttuu63K964u3bp1M/p1sSe0uJqIyZOVoE6YAIMGwRdfgJubpaPSaKpGQUEBFy5cYP78+ezdu7daAq8xjhZXE/LMM0pgn3oK+veHVatA13fW2AKbNm2ie/futG7dmpUrVxZvcNDUHC2uJmbMGCWwI0dCnz6wejXUr2/pqDSayomOjqagoMDSYdgVOjvTDIwYAQsXwsaN0KsXlMq60Wg0dQQtrmZi6FBYsgS2boX774caepxpNBobRYurGXn4YVWycMcOuO8+MEGFOY1GYyNocTUzAwfCl1/Cnj0QEwMXL1o6Io1GUxtYlbgKIQKFEIlCiANCiN+FEBOMtBFCiPeFEIeFEHuEEB0tEWt16NsXvv4a/vgDoqPBhCU6NRqNlWJV4grkAc9KKW8FugDjhBBtyrTpBYQWHqOAf9VuiDXj/vthzRo4dkwJ7OnTlo5Io9GYE6sSVynlGSnlb4X/zwAOAGUrP/QHFkrF/4AGQoimtRxqjejeHb79Fk6ehK5d4cQJS0ek0WjMhVWJa2mEEEFAB+DnMpeaA6Vl6STlBdhquece+P57OH9eCezx45aOSKPRmAOr3EQghPAAVgITpZRXy1428pZym5SFEKNQ0wb4+fmRlJRUrRgyMzOr/Z7q8H//58nkyRHccUc+s2btpnnzLLM9y9x9qQ28vb3JyMggPz/foFyfrWIv/QD770t2dnbNfn+klFZ1AM7Ad8CkCq7/G3i01OuDQNPK7hkZGSmrS2JiYrXfU1127pTS11fKZs2k/OMP8z2nNvpibvbv3y+llPLq1asWjsQ02Es/pLT/vhT97BkD2CEr0B2rmhYQQghgHnBASllRvbOvgScKswa6AFeklGdqLUgT0r49JCVBXp6aIvj9d0tHpLEH4uPjEULw9NNPW+T5c+bMISIiAi8vL7y8vIiKiipXM7YuYFXiCtwFPA50F0LsKjx6CyHGCCHGFLZZCxwFDgOfAGMtFKtJaNtWbZN1cFBZBLt3WzoijS3zv//9j08++YSIiAiLxRAQEMBbb73Fb7/9xo4dO+jevTsDBgxgz549FovJEliVuEopt0gphZQyQkrZvvBYK6X8WEr5cWEbKaUcJ6UMllK2k1LusHTcN0t4uBJYNzfo1g1+/dXSEWlMyYoVK3B1dSUlJaX43PPPP09wcDDnTJj0fOXKFWJjY5k3b57RqlYBAQHlCmDv3bsXNzc39u/fb7I4+vfvT69evQgJCSEsLIy4uDg8PT1NXvzb2rEqca3LhIbCpk3g5aV2cv1cNkdCY7MMHjyYdu3a8frrrwMwc+ZMvvjiC7799ttKzQiry6hRoxg8eDDdu3c3ej0qKort27cbnJs4cSIjR46kTRvDdPL4+Hg8PDwqPTZv3nzDmPLz81m6dCmZmZnceeedNe+cDWKV2QJ1lVatlMB26wY9esDatXD33ZaOyoqZOBF27ardZ7ZvX23bWSEE8fHx9OnTh+DgYOLi4li9ejWhoaEmC+uTTz7h8OHDLFq0qMI2UVFRxZYzAKtWrWLnzp0sX768XNuyLrPGqMx8cO/evURFRZGdnY2HhwdfffUV7dq1q0JP7Ac9crUyWrRQAtu0qdrVlZho6Yg0pqBnz5507tyZF198keXLlxMZGWm03YsvvogQotKjbFrQwYMHmTZtGgkJCbhU4i/UpUsXjhw5Qnp6Ojk5OUyePJkZM2bg6+tbrm2Ro2xlR71KKsG3bt2aXbt28b///Y+nnnqKYcOGsW/fvqp9sewEPXK1Qpo3V3OwMTGq4PZ//6tGspoyVHMEaUk2bNjA7t27kVJWOhUwceJEhg4dWum9WrRoYfB627ZtXLx4kbZt2xafy8/PZ9OmTXz88cf8+eefuLq6EhkZiYuLCzt27GDnzp04OTkxbtw4o8+Ij48nPj6+0jjWrVvHPffcY/Sai4tLsd13p06d2L59O++++y7z5s2r9J72hBZXK8XfX6Vp3XcfPPigqqzVu7elo9LUhN27dzNo0CA++OAD1qxZw9SpU/niiy+Mtm3UqBGNGjWq1v0HDBhAp06dDM49+eSThIaGMm3atOLRrKurKx06dGD16tUsWLCAxYsX41yBVfHNTguUpaCggJycnCq3twe0uFoxjRvDhg3QsycMGADLl6t/NbZDSkoKvXv3ZtKkSYwYMYLbb7+diIgINm/eTG8T/bVs0KABDRo0MDhXv359fHx8DEazoOZdZ8+eTY8ePejbt2+F9/Tx8cHHx6dG8UyZMoU+ffoQGBhIRkYGixcvJikpqc7luuo5VyvH1xfWr4eOHWHIEFV8W2MbpKen88ADD9C3b19mzJgBQNu2bRkyZAivvPKKRWJq3749Dg4O5VKyTMnZs2cZOnQorVu3JiYmhu3bt7Nu3Tp69epltmdaI3rkagM0aKCKvfTpA488Arm5EBtr6ag0N8LHx4cDBw6UO79s2TKz78WvaC98QkICo0eP5rbbbjPbs+fPn2+2e9sSWlxtBC8vWLcO+vWDxx+H69dh+HBLR6WxBQoKCrhw4QLz589n7969LFu2zNIh1Qm0uNoQHh7wzTdq3vXJJ9UIdtQoS0elsXY2bdpE9+7dad26NStXrjS6e0tjerS42hju7soyZvBgGD1aCayF6nNobITo6GgKCgosHUadQy9o2SBubrByJfTvD+PHgxnXJjQaTQ3R4mqjuLqqzIHBg+HZZ+GNNywdkUajKY2eFrBhnJ1hyRJwcYFp09QUwYwZIIx5NWg0mlpFi6uN4+QECxcqgX3lFcjJgbg4LbAajaXR4moHODrCvHlKYN94Q41g335bC6xGY0m0uNoJDg7w8cdKYN95Rwns7NlaYDUaS6HF1Y4QAt5/XwnsrFlqiuBf/7J0VBpN3USLq50hBMycqbIJiqYIblDBTqPRmAEtrnaIEGpRy9VVLXKdOHErXbuqxS+NRlM76DxXO0UIePllJbLr1/vx2GOqHoFGo6kdtLjaOdOmwVNPHWbFCnj4YTUPq7E/Nm3aRL9+/WjevDlCiAorU3300Ue0atUKNzc3IiMjq2QyaK1Utc+WQotrHeDhh0/y/vuwahUMGgTZ2ZaOSFNVhg8fXqXar5mZmbRt25bZs2dX6G21bNkyJkyYwLRp09i5cyd33nknvXr1IjU11cRR1w5V6bMl0eJaRxg/Hv79b+Uo268fXLtm6YjqDitWrMDV1ZWUlJTic88//zzBwcGcO3fOJM/o3bs38fHxDB48GAcH47/Ws2bNYvjw4fz973/n1ltv5YMPPqBp06b8q1RKSUBAQLlC2nv37sXNzY39+/ebJFZTUZU+WxLri0hjNkaNgk8/hR9/hL59ITPT0hHVDQYPHky7du14/fXXAZg5cyZffPEF3377baVmhaYkNzeXX3/9lZ49exqc79mzJ1u3bi1+HRUVxfbt2w3aTJw4kZEjR9KmTRuD8/Hx8Xh4eNC0aVM8PDyMHrY87XCz6PXjOsaTT6o82CeegF69YM0aVYjbFpn47UR2nd1Vq89s79+e9x6onuusEIL4+Hj69OlDcHAwcXFxrF69mtDQUDNFWZ6LFy+Sn59fTsz9/Pz48ccfi19HRUXx0UcfFb9etWoVO3fuZPny5eXuWWRimJmZiYeHh9HnVsfE0N7QI9c6SGwsLF0K//ufMj+8fNnSEdk/PXv2pHPnzrz44ossX76cyMhIo+2KRoNFR0JCQrlzNzMaFGW27EkpDc516dKFI0eOkJ6eTk5ODpMnT2bGjBn4+vqWu5ePjw8hISEEBwcTEhJi9LjRXOiLL76IEKLSoyLLGmtHj1zrKEOGqKpaDz+s7Lu//x5qaPZpMao7grQkGzZsYPfu3UgpK50KKGtp/cILL9C8eXOeeeaZ4nM1GQ02atQIR0dHzp49a3D+/PnzBvFERkbi4uLCjh072LlzJ05OTowbN87oPePj44mPj6/0uevWreOee+6p8PrEiRMZeoNdLi1atKj0urVideIqhPgU6Aucl1K2NXI9GvgvcKzw1JdSytdqL0L7YcAA+OoreOgh6NZNzcU2bmzpqOyP3bt3M2jQID744APWrFnD1KlT+eKLL4y2LWtp7enpWTxCvBlcXFyIjIzkhx9+YMiQIcXnf/jhBx566KHi166urnTo0IHVq1ezYMECFi9ejLOzs9F7mmJaoFGjRjRq1KgGPbJ+rE5cgfnAh8DCStpsllJWbLquqTJ9+ijbmP79ITpa2Xj7+1s6KvshJSWF3r17M2nSJEaMGMHtt99OREQEmzdvpnfv3iZ7TmZmJocPHwaUIWFqaiq7du3Cx8eneOQ3adIkHn/8cW6//XbuuusuPv74Y06fPs2YMWMM7hUVFcXs2bPp0aMHfftW/GtW9IcgIyMDT09Pk/WlqlSlzxZFSml1BxAE7KvgWjTwTXXuFxkZKatLYmJitd9jrVSlL4mJUrq7S9m6tZQnT5o9pGqzf/9+KaWUV69etXAkVSctLU2Gh4fLUaNGGZx/+OGHZefOnat0j2HDhsmXX375hu0SExMlUO4YNmyYQbs5c+bIli1bShcXF9mxY0e5cePGcveaP3++dHR0lPv27atSjJb6nlS1z9XBWF+KfvaMAeyQFelYRRcseVRBXNOA3cA64LYb3U+La2KV2m3eLKWHh5TBwVKmpJg3pupii+JaGdbcjx49esixY8dWub0196W6mFJchbpuXQghglCjU2Nzrl5AgZQyUwjRG5gtpSyX0yKEGAWMAvDz84tcunRptWKobB7J1qhOX/bv9+L55yPw9LzOrFm7adrUOrZzeXt7ExISQn5+Po6OjpYO56axtn4UFBRw8eJFEhIS+Oijj/jll1+qbMFtbX25GYz15fDhw1y5csVo+27duv0qpexk9GJFqmvJg0pGrkbaHgcaVdZGj1wTq9V++3YpGzaUMjBQyuRk88RUXfTI1bwkJiZKIYQMDw+XP/30U7Xea219uRlMOXK1xgWtShFC+APnpJRSCHE7Klc3zcJh2RWdOsGGDdCjB9x7r/p/eLilo9KYk+joaAoKCiwdhl1hdZsIhBBLgG1AayHESSHE34QQY4QQRUuag4F9QojdwPvAI4V/QTQmpH17SEyEggKVRbBvn6Uj0mhsC6sbuUopH73B9Q9RqVoaM9O2LSQlQffuJXmwf/mLpaPSaGwDqxu5aqyL8HDYtAnq1VMCu2OHpSPSaGwDLa6aGxISAhs3grc3xMSomgQajaZytLhqqkSrVkpgGzdWC111uJKcRlMltLhqqkyLFmqKoHlzeOABteCl0WiMo8VVUy2aNVMj2FatoHdvVU1Lo9GUR4urptr4+alRa+vW8OCD8M03lo5Io7E+tLhqakTjxmpzQbt2yvTwq68sHZF9Eh0dzdNPP33T90lKSkIIwcWLF00QlaYqaHHV1BgfH5X7Ghmpim8bcQLRWICgoCBmzpxpcO7OO+/kzJkzRh0FTMmZM2d47LHHCA8Px9HRkeHDh5v1edaMFlfNTdGggZp3jYqCRx+Fzz+3dEQaY7i4uODv71/O5sXU5OTk0KhRI6ZMmcIdd9xh1mdZO1pcNTeNpyd8+y107aqMDz/91NIRVUxCQgJBQUE4ODgQFBREQkKCWZ+3adMmunTpgoeHB97e3txxxx3sK7WX+Msvv6Rdu3a4uroSGBhIXFwcle3mNjYqLT11EB0dTUpKCs8991yxBxUYnxa40bODgoJ4/fXXGT16NF5eXgQEBPD2229X2t+goCDef/99hg8fbuCoUBfR4qoxCfXrKyfZnj3hb3+Djz+2dETlSUhIYNSoUaSkpCClJCUlhVGjRplNYPPy8ujfvz933303u3fv5ueff2bChAnFJe1+/fVXhgwZwqBBg9i7dy9vvvkmb7zxBh9+WPPd3V9++SUBAQHMmDGDM2fOcObMGaPtqvrsd999l3bt2vHbb7/xwgsv8Pzzz7Nt27Yax1eXsLraAhrbpV49WLVKzb8+9RTk5kIpXz2LM336dK5du2Zw7tq1a0yfPp3Y2FiTP+/q1atcvnyZBx98kODgYADCC8uLZWRkMGvWLLp27cqrr74KQFhYGMnJybz11luMHz++Rs/08fHB0dERT09P/Cvx66nqs3v27Fk8Kh4/fjzvv/8+69evJyoqqkbx1SX0yFVjUtzcYOVKGDgQJkyAMp9gLUpqamq1zt8sPj4+DB8+nPvvv58+ffowa9YsTpw4UXz9wIED3HXXXQbvufvuuzl16hRXr141S0zVfXZERIRBm2bNmnH+/HmzxmYvaHHVmBwXF1i2DP76V3juObiB+3KtUZFpnTnN7D777DN+/vln7r33Xr7++mvCwsL47rvvAFWovqIFporOOzg4lJuTvX79erXjquqzyzq/CiF03dcqosVVYxacnVXmwNChMH06vPIKWLrqblxcHO7u7gbn3N3diYuLM+tz//KXv/DCCy+QlJREdHQ0CxYsAKBNmzZs2bLFoO2WLVsICAio0E21cePGBvOo2dnZ/PHHHwZtXFxcyM/PrzSmmjxbUz20uGrMhpMTzJ8PTz4Jr74K06ZZVmBjY2OZO3cuLVu2RAhBy5YtmTt3rlnmWwGOHTvGlClT2Lp1KykpKSQmJrJnzx7atGkDwLPPPsvGjRt55ZVXOHToEAkJCbzzzjs8//zzFd6ze/fuJCQkkJSUxO+//86IESPKjVyDgoLYvHkzp06dqnDTQE2eXVV27drFrl27uHr1Kunp6ezatYv9+/ff9H1tjor8X+zp0B5aiRZ9fn6+lGPGSAlS/uMfUhYUVP8etuihdfbsWTlw4EDZrFkz6eLiIgMDA+Vzzz0nc3Nzi/uxcuVK2bZtW+ns7CwDAgLk66+/LgtKfYG6du0qx40bV/z6ypUr8pFHHpFeXl6yWbNmcs6cOeXabNu2TUZEREhXV1epfsVLbKgvXLhQ3O5Gz27ZsqV8++23DfpU9llSlv+eYMTuumXLljX8KtYudm+tbepDi2uipUOQBQVSPvOM+okbN04JbnWwRXGtDHvph5T235c6Y1CosU2EgPfeU4tdM2dCTg78+9/goCemNHaKFldNrSEE/N//gasrxMXB9eswbx7YieW9RmOAFldNrSIEvP66EtgZM9RGg4UL1eKXRmNP6B9pjUV46SU1RTBlihrBLl6s0rc0GntBi6vGYrzwghLYSZPUCHb5cjWi1WjsAb2coLEo//gHfPghfP212jKblWXpiDQa06DFVWNxxo2DuXNV2cJ+/aBMbRWNxibR4qqxCv7+d/jsM2Ud07s3ZGZaOiKN5ubQ4qqxGoYNU/UItmxR1t1mLgyl0ZgVLa4aq+LRR2HpUvj5Z+jRAy5dsnRElkUbFNouVieuQohPhRDnhRD7KrguhBDvCyEOCyH2CCE61naMGvMyeLCqCbtrF8TEQFqapSOyLSxpUPjll1/Ss2dPGjdujKenJ3fccQdff/21WZ9prViduALzgQcqud4LCC08RgH/qoWYNLVMv37K1WD/fujWDW5QQU9zA2rLoHDjxo10796dNWvWsHPnTnr37s3AgQPZvHmzWZ9rjViduEopNwHplTTpDywsrJvwP6CBEKJp7USnqU169YJvvoHDh+HcOZULe7Nog0LzGhTOnj2bKVOmcPvttxMSEsLLL79MZGQkq1atqvoXzU6wOnGtAs2BE6Venyw8p7FD7rsP1q2DvDw4dAiuX6/5yEsbFJZQmwaFGRkZNGzYsMZ9slVscYeWsd+ucn/qhRCjUNMG+Pn5kZSUVK2HZGZmVvs91oo99KVBg4bk5kpOnKgHZOLsXP2q21OnTjVqUDh16lT69etnokhLSE9P5/Lly8TExNCkSRMAmjdX44D8/Hzeeust7r77biZPngxAv3792LdvH2+++SbDhw8vbpebm0tGRgagSoTm5OQUvy7bxtnZGQcHB1xcXKhfvz6gxK2o35mZmbi6ulbp2VJKunXrxrBhwwAYPnw47733HmvXrqVt27YGzy8dT2nmzp3LyZMnGThwYIVtrAljfcnOzq7R748tiutJILDU6wDgdNlGUsq5wFyATp06yejo6Go9pMiSwx6wh74cOHCAwEDBoUMOnDzpQevW1d8qe+RVYFYAACAASURBVPLkyQrPm8PaxNPTk+HDhzNw4EBiYmKIiYlhyJAhBAYGkpGRweHDh+nTp4/Bs2NiYnjzzTeRUuLl5YWjoyMuLi7FbYQQuLq6GrynKm2K7G08PDzw9PSs0rOFEERGRhq0CQgI4MqVKwbnMjIyjH79Vq5cyUsvvcTSpUu57bbbbvbLWSsY64ubmxsdOnSo9r1scVrga+CJwqyBLsAVKaXxzz4au8LDAwIDr1FQAAcPQnZ29d6vDQpLqOqza2pQuHLlSh5//HEWLlxolk8FtoDViasQYgmwDWgthDgphPibEGKMEGJMYZO1wFHgMPAJMNZCoWosgJtbAWFhFAtsdWoRaIPCEsxpULh8+XKGDh3K/PnzGTx48E3dy5axOnGVUj4qpWwqpXSWUgZIKedJKT+WUn5ceF1KKcdJKYOllO2klDssHbOmdnF3h9at1f+rI7DaoLAEcxkULl26lNjYWN58803uvfdezp49y9mzZ0lPrywByE6pyP/Fng7toZVo6RBuGmMeWllZUu7aJeXOnVL++aelIquYumhQ2LVrV6MGhV27dr2Jr2TtYUoPLSEryamzFzp16iR37KjeANceFoGKsIe+HDhwgFtvvbXcgkN2tkrRys+HsDAoXCC3eipaBLJF7L0vRT97xhBC/Cql7GTsmtVNC2g01cHNTU0RODkpkdXVtDTWghZXjc3j6moosDaQTqmpA2hx1dgFLi5KYF1cIDlZlyvUWB4trhq7oUhgXV2VwF65YumINHUZLa4au8LZWS1s1aunCr5cvmzpiDR1FS2uGrujtMAeOaILbmssgxZXjV3i5KQE1t1dCWxdzGHXWBYtrhq7pUhgPTzg6FHtaKCpXbS4auwaR0cIDQVPTzh2DC5csHREmrqCFleN3VMksF5ekJIC589bOqKqow0KbRctrpo6gYMDhISAtzekpirbGHvFkgaFGzdu5M4778TX15d69eoRHh5eLpa6gi0Wy9ZoaoSDAwQHq/nXEydASvD3t3RUtUORQaG58fDw4JlnnqFdu3a4u7vz008/MXr0aNzd3Rk7tm5VB9UjV02dYsmSBGJigrj9dgc6dAjiww+1QaEpDQojIyN55JFHuO2222jVqhVDhw7l/vvv1+6vGo09U2RQmJqqDArPnk3huedG8eGHCZijOJw2KISdO3eydetWunbtWuM+2So3nBYQQoRLKf+4UTuNxtqZPn16OYPC7OxrxMVNZ8CAWJo3hwqcT2rE1atXuXz5Mg8++CDBwcEAhIeHA6q03axZs+jatSuvvvoqAGFhYSQnJ/PWW28xfvz4Gj3Tx8cHR0dHPD09K50GqOqze/bsWTwqHj9+PO+//z7r168nKiqq0jgCAgK4cOECeXl5vPzyy4wZM6bS9vZIVUauO4UQs4UQdc8bV2NXpKamGj1/7lwqZ8+WzMOaCh8fH4YPH879999Pnz59mDVrFidOlLjCHzhwgLvuusvgPXfffTenTp3iqpkrz1T12REREQZtmjVrxvkqpFts3ryZHTt28PHHH/Pee++xaNEi0wRuQ1RFXG8HbgOShRDjhRCOZo5JozELFRkRBga2oEkTlaKVmmpaga2rBoWtWrWiXbt2/P3vf2fSpEm88sor1Y7R1rmhuEop90op7wNGAs8Ae4UQvcwemUZjYioyKIyPjyMwUGUOXLigcmFNKbB1zaCwLAUFBeTk5Jj0nrZAlRe0pJSrUCPYBcBSIcRaIUS42SLTaExMZQaFQkDz5tC0KVy8CMeP37zA1kWDwg8++IBvvvmG5ORkkpOTmTdvHjNnzmTo0KE3dV9bpLp5ru7AryiBHQfsEUJ8DLwkpdTVMzVWT2xsbIVur0UCKwScPq3ENShI5cfWBHd3dw4dOsSQIUO4ePEifn5+xMbG8sILL5CdnU3Hjh1ZsWIFL7/8MvHx8fj5+TFlypRKd2RNnTqV48eP079/fzw8PJg+fTqnT582aPPaa68xevRogoODycnJMZraVZNnV4X8/HxeeOEFjh8/jpOTE8HBwbz55pt1ckHrhgaFQoiJQOfCIxjIBXYB2wr/jQXaAoOklD+bNdoaog0Kbb8vFRkUmouzZ+HkSWjQAG65peYCWxH2bupnq5jSoLAqI9dnUUL6L+B/wK9SytxS1xcKIV4APkVNG2g0No+/vxrBnjihShYGB5teYDX2zQ3FVUoZWIX7fAbE33w4Go314OenBDY1VbkahIRogdVUHVP9qFwAupvoXhqN1dCkCbRsqQwPk5PhBovwGk0xJhFXqdhointpNNZG48bQqpWy7NYCq6kq+kOORlMFfH3VwlZmJhw6BHl5lo5IY+1YnbgKIR4QQhwUQhwWQkwxcn24EOKCEGJX4THSEnFq6h4+Pmph69o1NYLVAqupDKsS18KttXOAXkAb4FEhRBsjTZdJKdsXHv+p1SA1dZqGDUsE9tAhqMHOU00dwarEFVXH4LCU8mhhutdSoL+FY9JoDGjQQGUOZGVpgdVUjLWJa3PgRKnXJwvPleUhIcQeIcQXQoiqpIppNCbF21v5cuXkwMGDkJt74/dorJjcXDh7FscyJSlvBmuzeTFWpqfsFrLVwBIpZY4QYgxqK265NDAhxChgFICfnx9JSUnVCiQzM7Pa77FW7KEv3t7eZGRkkJ+fT0ZGhqXDAYq2yzpy8mQ9/vijgICALJydq1aQwJL9GDNmDGlpaaxYseKm7pOSkkK7du3YsGEDnToZ3aRkUsaMGUNQUBBTppRbijFKUXxJSUl07Nix3Ovf9+5l0KBBHFizBi8pEYBDw4ZklCnuk52dXbPfHyml1RxAFPBdqddTgamVtHcErtzovpGRkbK6JCYmVvs91oo99GX//v1SSimvXr1q4UjKk5Eh5a+/Srlnj5TZ2YbXhg0bJvv06VPuPUX92LVrl+zXr5/08/OTrq6uMjAwUA4aNEgeP37cbPFWFFNldO3aVY4bN87gXF5enjxz5oxMT083ZXhG2bNnj2zQoIG8cuVKld9z7NgxCcjt27cbvk5MlPLoUSl//VUO6tZNvjZ2rJQnT0qZlWX056voZ88YwA5Zge5Y27TAdiBUCNFKCOECPAJ8XbqBEKJpqZf9gAO1GJ9GUw4PD2jdWmUPHDwI2dlVe9+FCxeIiYnBw8ODNWvW8Mcff7Bo0SKCg4PNXizbFDg6OuLv74+Tk/k/AH/wwQc89NBDeHl51ewGWVlQVKoxNRUuXQIfH54cM4Z/ffUVeX5+4OZmuoCxsjlXKWUe8DTwHUo0l0spfxdCvCaE6FfY7BkhxO9CiN2o+rLDLROtxhZJSCipdBUUpF6bgvr1ISwMCgqqLrA//fQTly5d4rPPPiMyMpKgoCC6du3K//3f/9GuXbsK37d3715iYmLw8vLC09OTv/zlLyQmJhZf37RpE3fccQdubm74+fnxj3/8g9xKJoVLGxwWMXz4cPr27Vv8/40bNzJnzpxi08Pjx49z/PhxhBD89ttvVX52dHQ0Y8eOZdq0aTRq1IgmTZowefLkSgtw5+fns3z5cvr162dw/vPPP6dz5854enrSpEkThgwZwqlTp0oaFK00Hj0Kv/+uivWCKn3Wvj0EBdFzwADS09PNMm1mVeIKIKVcK6UMk1IGSynjCs/NkFJ+Xfj/qVLK26SUf5FSdpPa30tTRRISYNSokmLYKSnqtakFVkolsFlZlbf39/enoKCAL774olLH17I89thjNG3alF9++YWdO3fyyiuv4FY46jp16hS9evWiQ4cO7Ny5k3nz5rFkyRKmTp1a437Nnj2bqKgonnzyyWLTw8DA8uvIVX12QkICTk5ObN26lQ8//JD33nuPZcuWVfj8PXv2cOXKlXLzurm5ubz66qvs3r2bb775hosXL/Loo49CerpKRD5Q6kNtYCAUVbby9i4uEuHi4kL79u3ZuNH0G0ytTlw1GnMxfbrKTy3NtWvqvKlwd1dTBKAEtrKNBl26dGHatGkMGzYMHx8fevbsSXx8PCkpKZU+IyUlhR49ehAeHk5ISAgDBw4sNgz86KOPaNq0KR999BG33norffv25c033+TDDz8sZ85YVby9vXFxccHd3R1/f3/8/f2LHWxLU9Vnt2nThtdee42wsDAefvhhunXrxvr16yvtrxCinOHiiBEj6N27N7e0asXtt97Kv158kc2bN3Ny2zb1jW3cWDW85RZVhaeMZU0RzZo14/jx49X/wtwALa6aOkMF/oQVnq8p9eopgRVCFXyprBZBXFwcZ8+eZe7cubRr14558+bRpk2bSsVm0qRJjBw5ku7duxMXF2dg83LgwAGioqJwKFW+6+677yY3N5fDhw+bpH8VUdVnV9f0MCsrC2dnZ4P7Avy2bRv977+fls2b4+nvT6cHHwQg1cUFIiKUrUQVqFevHlk3+phRA7S4auoMFfgTVnj+ZnBzg/BwJbCZmfDnnxW39fX1ZciQIbzzzjscOHCAoKAg/vnPf1bY/pVXXmH//v0MGDCArVu3EhERwaeffgrYp+lho0aNyM3NVaPfvDw4f54/f/2V+3v3xh1Y9PbbbF+/nm/XrQMg19m5Wh7p6enpNC4a5ZoQLa6aOkNcnPrYXhp3d3XeHLi6gpeXmt47eFCJ7I1wcXEhODiYzBs0Dg0N5ZlnnmHNmjX87W9/4z//UbvA27Rpw7Zt2wzEasuWLcX3NUZZ00OA3bt3l4urKqaH1X12VWhfONLd/913sHs3pKbyx5EjXLx8mfj33+fe2FjCu3ThfFpaje6/b98+OnbsWOP4KkKLq6bOEBsLc+eq+qxCqH/nzlXnzYWjIxQUXOXo0V2sXr2Ln37axa5du0hJSeGbb75h6NChfPPNNxw6dIiDBw8yc+ZM1q5dy8CBA43eLysri3HjxpGUlMTx48f5+eef2bJlS7Hp4dixYzl9+jRjx47lwIEDrFmzptgbq6zzbRHdu3dn3bp1fP311xw8eJBJkyZx4sQJgzZBQUH88ssvHD9+nIsXLxodadbk2RUipfprlJJC49On6RgezpatW1WB3TZtaBEdjaurKx/OncvRo0dZs2YNL730UvWeARw/fpxTp07Rs2fPar/3RljbDi2NxqzExppXTI2xZctmtmzpYHCuf//+zJo1Cw8PDyZPnsyJEydwcnKiVatWzJw5kwkTJhi9l6OjI5cuXWLYsGGcPXsWX19f+vbty8yZMwFo3rw569at47nnnqN9+/Y0aNCAxx57jPj4io1CRowYwZ49exgxYgSgRHLgwIEGzrGTJ09m2LBhtGnThqysLI4dO1buPjV5djlyciAtTR05OeqvYMOGjBo9mv8kJDDx7bcBaOzuzoIFC5g2bRpz5swhIiKCWbNm8cADD1T9WcCSJUvo2bMnLVu2rNb7qsINDQrtAW1QaPt9qW2DQnNw/boq9JKdDc2aXaNp02qO5qyUm/6e5OWppP60tJK5E09PVUS3YUNwdCQnJ4fw8HAWLlzIPffcY5K4c3JyCA0NZcmSJdx1111A7RsUajQaE+DsrLIIDh2C06frUa+eqrBVJykoUKkUaWlw+bKaBnBzUwn+Pj5qwroUrq6uLFiwgPT0dJOFkJKSwvTp04uF1dRocdVoahEnJ7XR4I8/CjhyxJFbblGDszqBlCr/NC1NJfrn5akvSOPGapTq7l7pKv+9995r0nDCwsIICwsz6T1Lo8VVo6llnJwgIOAaZ896cuSI8ufy9bV0VGYkN7dkHjU7Wwlogwaq00XpFHaIFleNxgI4Oqp6sIcPw7FjalDXqJGlozIh+fkl86hFpRU9PFSKRsOG6i+MnWP/PdTYDfa2+OroqBwNjhyB48eVwJohl732kNJwHrWgQM2dNmumRqll5lFtgZv5mdPiqrEJnJ2dzbJF0dKUFtiigjJNmlg6qurhkJ2txDQ9XaVEODoqMfX1VdVsqrFbytq4fv16jUsqanHV2ARNmjTh1KlTNGzYEA8Pjwq3WdoiDg7K9PDoUVXnQEpVZ8Sqyc1VYpqWRv2sLCWg3t5KUEtVnbJlCgoKOHfuHN7e3jV6vxZXjU1QVCT58OHDXCiqy2nDZGdnF5cJLKJoMX3nTrXeU8PfafNRUKDqKGZmlhSsdXHhuqsrzt7eatR69qw6bJSy35f69evTqIaT4VpcNTaDl5cXV69eNcs+8NomKSmJDh06lDuflwdPPAFLlsBrr0ENdnSalvx8SEqCRYtg5UolrEFBMHQoPP44hIWpTSpt21o4UNNQ0felJmhx1WisCCcnpWPOzjBjhtoB+s9/WmDa8vffVSAJCXDypEqZ+utflfLffbddfOw3N1pcNRorw9ERPvsMXFxUxa7cXHjrrVoQ2PPn1ZB54UL47TcVyAMPwMyZ0K+fKlSrqTJaXDUaK8TBAf79bzWCffttJbDvvmsGgc3Ohq+/VoL67bdqGqBjR/WwRx+1gZU160WLq0ZjpTg4wJw5Kj30vffUFMGcOSb4RF5QAFu2qI/9y5er3NSAAJg8Wc2j3nabSeKv62hx1WisGCFg1iwlsG+9pUawc+eqT+zV5tAhJaiff652LdSvD4MHK0GNjq7hTTUVocVVo7FyhIA33lAC+9pragQ7f34Vd5CmpcGyZepj/88/q2HvffepVbKBA5XAasyCFleNxgYQAl59VS1yvfiiSin9/PMKDE1zcmDtWiWoa9aoxu3aqcnbxx5T21E1ZkeLq0ZjQ0yfrkawzz0HR48mcP78dE6cSKVFYCBxTz5J7PnzsHSpKpri7w/jx6uP/e3bWzr0OocWV43Gxpg8GfbsSWDRolHANQBSUlMZ9eqr4OxM7JAhSlDvu69OVJ+yVvRXXqOxNS5dYtOaZygS1iKuAdP9/YlNSLBIWBpD9DYLjcYWuH4dVq+GIUOgaVNSK7A7ST15spYD01SEFleNxlqRErZvh2eeUYtQ/fqpff6jR9PC39/oWwICWtRujJoK0dMCGo21kZqq9vQvXAh//KFWsB58UO3rf+ABcHYm7vbbGTVqFNeulZ4acMfFJY7Ll+uw8aEVYXUjVyHEA0KIg0KIw0KIKUauuwohlhVe/1kIEVT7UWo0JiYjQyWvdu+urFCmTVO+L//+tyrht2KFEtjC3KvY2Fjmzp1Ly5YtEULQsmVLJk6cS2pqLPfdp0qtaiyLVY1chRCOwBygB3AS2C6E+FpKub9Us78Bl6SUIUKIR4C3gL/WfrQazU2Slwc//KB2Ta1apWqlhoSonQJDhyrnwkqIjY0lNjbW4FxMDDz0kNLoH3+0M18uG8PaRq63A4ellEellLnAUqB/mTb9gQWF//8CiBH2VJZeY99ICbt2EfzRRxAYCL17q4Ipw4fDtm1qi+pLL91QWCuib1+17nXwoNrReu6cSaPXVANhTaZvQojBwANSypGFrx8H7pBSPl2qzb7CNicLXx8pbHOxzL1GAaMA/Pz8IpcuXVqtWDIzM/Hw8LiZ7lgNui+Wx+XiRfzWr8fv++/xOHqUAicn0rp04VyPHqR16YJ0cTHp83bubMC0ae1o0iSbd97ZTaNGuSa9f2ls9XtijOr2pVu3br9KKTsZvSiltJoDGAL8p9Trx4EPyrT5HQgo9foI4FvZfSMjI2V1SUxMrPZ7rBXdFwuRmSnlokVS9ughpYODlCBlly5SzpkjN69aZfbHb9okpYeHlCEhUqammu85NvU9uQHV7QuwQ1agO9Y2LXASCCz1OgA4XVEbIYQT4A3o6XuNdZCfD+vXw7Bhqhbq449DcrJaoDp4UH30HzuWvFowyLrnHjWle/48dO2qCmFpag+rWtACtgOhQohWwCngEeCxMm2+BoYB24DBwIbCvyAajeX4/XeVOpWQAKdOKXfBRx9V4mpBW5QuXZTW9+gB994LiYnKaVZjfqxKXKWUeUKIp4HvAEfgUynl70KI11DD76+BecAiIcRh1Ij1EctFrKnTnDunbFEWLSqxRenVSxVg7dcPyri7WopOnZSo3nefEtgNG6B1a0tHZf9YlbgCSCnXAmvLnJtR6v/ZqLlZjab2ycpStiiLFpXYokRGwuzZ8Mgj0KSJpSM0Svv2anNXTIyaIli/XhsOmBurE1drIDM3k3yZb+kwNNZCQQFs3qwEdcWKEluU555TH/vbtLF0hFWibVvYuFHlwEZHqzzYv/zF0lHZL1pcjTD5+8l8+tunhB0II9Q3lFCfUMJ8w4r/9ffwR6fW1gEOHiyxRUlJAQ8PlaH/xBNKnWzQXjo8vERgu3VTC16RkZaOyj7R4mqEfq37cfncZbLrZ3Pw4kHWJq8lN78kT9DDxYNQn9BywhvqG4pvPV8tvLbMxYsltii//KIEtEcP5XE9YIBd2KKEhsKmTUpcY2Lgu+/gjjssHZX9ocXVCL1De+N+yp3o6GgA8gvySb2SSnJ6MslpyRxKO0RyejK/nv6VlftXGkwhNHRrSKiv4Ui3SHi9XL0s1CNNpeTkKDuUIluUvDyIiICZM5UtStOmlo7Q5LRqpQS2e3f1t2PtWpXUoDEdWlyrgKODI60atqJVw1b0DO5pcC03P5fjl48rwU1LJjldie+mlE0k7ElAUpIl5lffTwmvT5iBAIf4hFDPuV5td6tuI6XKOV20SI1Ui2xRJkxQ86h1YDKyRQs1RRATo4ptffONmu3QmAYtrjeJi6MLYb5hhPmGlbuWdT2LI5eOGIx2D6UdYk3yGs7tMtz0HegVaCC8RaPdWxregoujabdG1mmOHlWCumgRHDkC9eopF9QnnlAqU8dsUZo3L8ki6NVLJUL06GHpqOyDuvWTVMvUc65H2yZtadukbblrV3OuFo90k9OSOZSuRr7Lfl/GpexLxe0chANBDYLKTTOE+YbRwrsFjg7aa/6GXLoEy5crQf3pJ2Wl2q2bKpAyaBB4elo6Qovi768E9r77VFXDL79U9WQ0N4cWVwvh5epFZLNIIpuVX6pNu5ZWLLrJ6ckcTDtIcloyW1K3kJmbWdzOxdGF4IbBJSPdIvH1DaWZZzMchO2tZpuM3FyVh7pokRqO5eaqlKk33oDYWFWRSlNM48Zqc0HPnmrdbsUK6F+2Hp2mWmhxtUJ83X3xdfelS0AXg/NSSs79eY7ktBLBTU5Xx/dHvic7L7u4rbuzOyE+IYT6hOKa6cox72PF87yN3RvbZ0aDlLBjhxLUJUvUyn/jxjBmjPrY37GjGrVqjOLrqzYXPPAADB4Mixcryy5NzdDiakMIIfD38Mffw597Wt5jcK1AFnDy6sni+d2iOd695/dyJP0Ii08sLm7r7eptdLQb5htGAzcb9AdJSSmxRTl4UNmi9O+vFqbuv7+4er/mxjRoAN9/D336qA1n16+rhAlN9dHiaic4CAdaeLeghXcLYm6JMbi2PnE9Lf/SslxGw7aT21i6b6lBRkMj90Yl6WOlhDfUJ5T6LlaU43n1KnzxhRqlJiWpc/feC5Mnq2GXNpGqMV5esG6dKo8wdKiaURk+3NJR2R5aXOsAjsKREJ8QQnxCINTwWnZeNscuHTMY7R5KO8QPR39gwe4FBm2beTYzurB2S8NbcHVyNX9H8vLw+fln5Su1ahVkZ6uM+CraomiqjoeHSs0aMACefFIJ7KhRlo7KttDiWsdxc3Lj1sa3cmvjW8tdy8zN5HD64ZKphsKMhq/++IqL10qMH4pGzWVFN9Q3lKAGQTg53PjHLCEhgenTp5OamkqLFi2Ii4tT/lBSwu7d6iP/4sVEnDsHPj4wYoT62H/HHXoe1Uy4u6u1wIcegtGjlcA+/fSN36dRaHHVVIiHiwft/dvT3r99uWuXsi5xOP2wwWg3OT2ZRXsWcTXnanE7Jwcnbml4i9ERb3Ov5jgIBxISEgxsolNSUhj197/DN98Qu28f7Nun5k379mVfx460ff55MLEtisY4bm4qNeuvf4Xx45XATppk6ahsAy2umhrRsF5DOjfvTOfmnQ3OSym5cO1CuYW1Q2mHWH90PVl5WcVt6znVI9gnmCMzjpB1LcvgPteyspi+dCmxXbrARx/Bww+Dry8Xk5K0sNYyrq4qNSs2Fp59VgnslHKm95qyaHHVmBQhBE3qN6FJ/Sbc1eIug2sFsoBTV08Z5PAmpx1i38V9Ru+VAjzybEvCfM8QenItYVlhXL1+1WhbjXlxdlapWc7OMHWqKscwY4aekakMLa6aWsNBOBDoHUigdyDd/2wCPx6BhF8JQglpWdx83dh+ejsr9q+gQBYUn/fZ6WOQ0VC6ToOnq/3vtkpIgOnTITVV1QeIi1OjSnPj5KSmvl1c4JVX1Aj29dfN/1xbRYurpvY4e7bEFmXnTvXb+sADxD30EKPmzSuecwVwd3dn7uy5xMbGkpufy9FLR0lOS2bdL+soaFhAcnoySceTWLRnkcEj/D38jS6sBTcMtoviOAkJatW+6EuVklKyil8bAuvoCPPmKYGNj1cCq7fKGkeLq8a8ZGXBf/+rhjzff19ii/L++ypLvXFjYgHuuMN4tgBqm294o3DCG4XjecazuBQkwLXr14ozGkovrK0+tJrzf54vbicQBHoHGi0F2apBK5wdbWOjwfTpJcJaxLVr6nxtiCuoErcff6wEduZMOHIkhOhoPUVQFi2uGtNTUKCKhRbZomRkqL38zz+v0qduLZ/2FRsbWyym1cHd2Z0Ivwgi/CLKXbuSfcVoRsOSfUu4nH25uJ2jUCUlje1YC/QKtKriOKmp1TtvLoRQfx9dXGDWrADGjIF//csmzRnMhhZXjen44w8lqAkJJbYogwerff1du9b6b563m7fR4jhSStKy0sqNdovq8P55/c/iti6OLsU1GspONTT1aFrrNRpatFBfWmPnaxsh1Mj17NkU5s5tSW4u/Oc/aupAo8VVc7NcvAhLl6qP/du3l9iixMer7T3u7paOsBxCCBq5N6KReyOiAqMMrkkpOZN5xqjwrju8zsDup75zfaNWP6E+oTRyb2QW4Y2LM5xzBfUljosz+aOqhBAwcuQx1UffhQAAEFFJREFUwsJaFi9yLVhQ58riGkV/CTTVJydH7Y1cuFD5g+Tlqcr9dmCLIoSgmWczmnk2o2tQV4Nr+QX5nLh6olyNhp1nd/LlgS8N7H4auDUwOtoN9QnF2827xvEVzZxYIlugIoSAl19WaVrTp6tiLwkJul6OFldN1ZAStm4tsUW5fFmJ6IQJ6mN/RPk5T3vD0cGRoAZBBDUIKmf3cz3/eondT6kC6FtSt7B472KD4jhN6jehiWMTOl3pZGD5E+ITgrvzjUf6sbGWFdOKmDZNbTiYPFkJ7NKl6nVdRYurpnKOHCmxRTl6VH0GHTRILUzFxOgJtkKcHZ3VyNQ3tNy17LxsgxoNyenJbD+6nW8Pf8v8zPkGbZt7Ni832i0qjmMLdj/PPqsWuZ55Rv2YrFypttDWRbS4aspTZIuycKEarQqhbEJnzNC2KDXAzcmtnN1PUlIS0dHRZORkGGQ0FI16Vx5YSVpWWnH7IrsfY1MNLb1bWlVGw/jxasQ6erQqW7hqlVVOvZsdLa4aRW6uKuK5aBGsXq1tUWoJT1dPOjTtQIemHcpdS89KNxjtFs3xbj2xlYzcjOJ2zg7OBPsEG92xVlQcp7YZNUrNuf7tb6rw9urVKnmkLqHFtS4jpVrhL7JFSUuDJk1g7Fj1sb9DB50ZbkF86vlwR8Ad3BFwh8H5Iruf0gtrpevwlrb7qedUjxCfEKMLa03qNzFrKtmTT6opgieeUNYxa9eqQtx1BasRVyGED7AMCAKOAw9LKS8ZaZcP7C18mSql7FdbMdoNKSnw+efqY/+hQ+oz3IABSlB79tTLvFZOabufe1vea3CtqDhO2Y0T+87v478H/0teQV5xWy9XL6NpZGG+YTSs19AkscbGKoF97DH1o/Xtt3XHJMJqxBWYAqyXUr4phJhS+PoFI+2ypJTlC4xqKqfIFmXhQti4UZ279161a2rwYPCueXqQxnooXRynrN1PXkEeKZdTymU0GLP78a3nW7KgViajwcOlep/vhwxRAjtkiLLv/v57Ve/c3rEmce0PRBf+fwGQhHFx1VSVvDz4/ntufecdtTCVnQ1hYfDPfypblKAgS0eoqUWcHJwI9gkm2CeYXvQyuJaTl6OK46Qb1uH98eiPLNy90KBtM89mBiPenIs5ND7fmGCfYNycjKcG9O+vFrYGDVJroz/8oIx57Rkhpbxxq1pACHFZStmg1OtLUspyn02EEHnALiAPeFNKuaqC+40CRgH4+flFLl26tFrxZGZm4mGLM/BS4nH4MH7ff4/f+vW4XLpErqcnF2JiONujBxm33mrT86g2+30pgy31Iys/i1NZpziZdZKT106qf7NOcirrFJevl9RoEAiauDYhoF4AAe4B6t/Cw9/NHycHJ3bsaMiLL7bF03MBMI20tPM0adKEkSNHct9991muk4VU9/vSrVu3X6WUnYxdq1VxFUL8CPgbuTQdWFBFcW0mpTwthLgF2ADESCmPVPbcTp06yR07dlQr1qJUGZvh1KkSe+nff1fzpg8+CI8/zsb69enao4elIzQJNvd9qQB76cfl7Mss/X4pXq28iqcZihbaruRcKW7n5OBEqwatCPMN48KGfH75bAMUlGwldnd3Z+7cuTUq3mNKqvt9EUJUKK61Oi0gpazwT5MQ4pwQoqmU8owQoilw3lg7KeXpwn+PCiGSgA5ApeJqt2RmKoOjRYtg/Xq1+h8VpWxR/vrX4oktWWQ9rdGYmAZuDQj3Cie6XbTBeSklF69dLLewlpyWzJ4v9kCB4aDu2rVr/H3i39nTZI/Bwpq/h3+tF8cxFdY05/o1MAx4s/Df/5ZtIIRoCFyTUuYIIRoBdwH/V6tRWpr8fNiwQY1Qv/xSVfBo1Uol+A8dCiEhlo5Qo0EIQeP6jWlcv3E5ux+HscbzbrMuZvHez+8ZFMfxcPEoyd0ttbAW6hOKr7uvWftws1iTuL4JLBdC/A1IBYYACCE6AWOklCOBW4F/CyEKAAfUnOt+SwVcq+zdW1LO7/RptbofG6uSCO+6y6bnUTV1ixYtWpBipG5is2YtSZ12hNQrqSWj3cKphh2nd/DF/i8M7X7q+RhNJQv1DcXL1fIJtVYjrlLKNCDGyPkdwMjC/28F2tVyaJajyBZl4ULYtUvVcevVC2bPhr596+6mbY1NExcXZ2ClrnDnzz/jSD7kSHh4K1o1bFWuOE5ufi7HLh0zEN6K7H786vsZHe2G+ITUmt2P1YirppBr15QtyqJFJbYonTsb2KJoNLZM0aJVaVufp56K4913Y+naVS0ftG1b/n0uji60btSa1o1al7uWdT2LI5eOlIx2C+d41ySv4dyucwZtA70Cje5Ya9WwlUn7qcXVGiiyRVm4UCX6l7ZFeeIJCA+3dIQajUkxZuvTv7/KgY2Ohh9/hPbV2CpUz7leueI4RVzNuWq0RsOy35dxKbtkE6ijcOTxFo8TXZxuf3NocbUkRbYon3+uKh97eqrdUo8/bhFbFI3GkoSHqzFG9+7q+P576GQ0yal6eLl6GbX7AUi7lmaQyeB52XQV37S41jYXLqgqwosWldii9Oypqk9ZqS2KRlNbhIQoge3WTZUL/u476NLFfM/zdfclyj2q2O4nyYRpi3poVBtkZ6uP+/37Q7NmqpLw9eswa5ZK/l+3TlW20MKq0RAUpAS2SRNlx7Z5s6Ujqhl65GoupISfflIj1OXLlS1Ks2bwj3+oj/3t6k7Sg0ZTXQIDVX2hmBhVrnD1ajVVYEtocTU1hw+XzKMW2aI89JAS1O7dtS2KRlNFmjWDpCQlsH36qCSanj1v+DarQYurKUhPV6PTRYtKbFFiYpQl5qBBda8Eu0ZjIvz8lMD26KFKZXz5pRJaW0DPudaU3FxVQ+2hh5QL6lNPwZUr8Oab/9/eHcdaXdZxHH9/hIswpNsYjsgr2h/J5pxWuyONLRlaoUgFQ8fdSDcZ4EYbyTCou+Fg6pzTLGfTEHJmEjSN1YIyGt7QLEAETSMLGxXTRS2R7tzA8tsfzymge67ee8753eec3/28trvdc+45v/N9uLsfnvP8nvP7pjP/27enbVQOVrO6TJiQ9r5efDHMmQNbtuSuaGA8cx2MCNi9O81QN206vS3K9denjXn+GKpZw40fn/a+zpyZLrq9cSNcd13uqt6dw3UgDh1Ka6iPPpraoowefXpblJH+ZzQrWnt72vs6axZ0daU3jwsW5K6qf06Ffozo7YUNG9KnpnbuTHdefjmsXJmWAtwWxWzIjRuXdi7Onp3eLJ44ATfemLuq6hyu1XR384m7706/uQsugNtuS1egclsUs+zGjoWtW9Obx4UL05/pTTflrqovh2s1kyfz+qxZdKxalS6a4nVUs6YyZkzamjVvXjqXfOJE+mxOM3G4VrNkCQenTKFj6tTclZhZP0aPTluz5s+HZctSwK5Ykbuqk7wVy8xa1qhRsHlz2jlwyy1wxx25KzrJM1cza2ltbalBR1sbdHenGeytt+ZfzXO4mlnLGzkSHnkkzWTXrIHjx9MsNmfAOlzNrBRGjID161PA3nlnCth77skXsA5XMyuNM86ABx5IAXvvvWmJ4L778lx33uFqZqUipR6eZ54J/92u/uCDQx+wDlczKx0J7rorBeztt6eA3bBhaK/46XA1s1KS0ocrR41Kuwfefjud9BqqS4E4XM2s1FavTjPYVavSDHbjxrRtq2gOVzMrvZUr0wx2+fI0g928OQVukfwJLTMbFm6+Ge6/P12TYO7c1De0SA5XMxs2li6FdetOXrbwrbeKey2Hq5kNK4sWwcMPw44d6cLbvb3FvE7ThKukayW9LOkdSZ3v8riZkl6RdFDSqqGs0czK4YYbUnORp59OrWOOHWv8azRNuAIvAXOBnf09QNII4JvAVcCFQJekC4emPDMrk66u1Apv167UXfaNNxp7/KYJ14g4EBGvvMfDpgIHI+KPEXEC2AR8rvjqzKyM5s2Dxx+HffvgiivgzTcbt4Gq1bZinQP85ZTbh4GPV3ugpMXAYoCJEyfS09MzqBfq7e0d9HOalcfSfMoyDmj9sbS3w9q141m9+iIeeqiD9vaehhx3SMNV0s+BD1T5UXdE/HAgh6hyX1R7YESsA9YBdHZ2xvTp0wdaJgA9PT0M9jnNymNpPmUZB5RjLNOnp6+jR//csLEMabhGxJV1HuIwcO4ptzuA1+o8ppkZl14KPT3vNOx4TbPmOkB7gA9L+pCkUcB84EeZazIz66NpwlXSHEmHgcuArZKerNz/QUnbACLiX8AXgSeBA8D3I+LlXDWbmfWnaU5oRcQWYEuV+18Drj7l9jZg2xCWZmY2aE0zczUzKxOHq5lZARyuZmYFcLiamRXA4WpmVgCHq5lZARyuZmYFUETVj+aXiqS/AX8a5NMmAH8voJwcPJbmU5ZxwPAey3kRcXa1HwyLcK2FpOciot+LdrcSj6X5lGUc4LH0x8sCZmYFcLiamRXA4dq/dbkLaCCPpfmUZRzgsVTlNVczswJ45mpmVgCHaxVlad8t6duSjkh6KXct9ZB0rqSnJB2otF9flrumWkkaLWm3pBcqY1mTu6Z6SBohaZ+kH+eupR6SDkn6jaT9kp5ryDG9LHC6Svvu3wOfIrWV2QN0RcRvsxZWA0mfBHqB70TERbnrqZWkScCkiHhe0jhgL/D5Fv2dCBgbEb2S2oBngGUR8evMpdVE0nKgE3hfRFyTu55aSToEdEZEw/breubaV2nad0fETuAfueuoV0S8HhHPV77/J6kLxTl5q6pNJL2Vm22Vr5ac4UjqAGYB63PX0owcrn1Va9/dkn/IZSTpfOCjwK68ldSu8lZ6P3AE2B4RrTqWrwNfBhrX1S+fAH4maa+kxY04oMO1rwG377ahJeks4AngSxFxLHc9tYqIf0fER0jdi6dKarklG0nXAEciYm/uWhpkWkR8DLgKWFpZUquLw7Uvt+9uQpX1ySeAxyLiB7nraYSIOAr0ADMzl1KLacBnK2uVm4AZkr6bt6TaVXr1ERFHSL38ptZ7TIdrX27f3WQqJ4E2AAci4mu566mHpLMlvb/y/RjgSuB3easavIj4SkR0RMT5pL+RHRGxIHNZNZE0tnKiFEljgU8Dde+wcbj+nzK175b0PeBXwBRJhyUtzF1TjaYBXyDNjvZXvq5+ryc1qUnAU5JeJP1Hvj0iWnobUwlMBJ6R9AKwG9gaET+t96DeimVmVgDPXM3MCuBwNTMrgMPVzKwADlczswI4XM3MCuBwNTMrgMPVzKwADlczswI4XM0ASddKOi7pvFPu+4akVyVNzFmbtSZ/QsuM/12/YA+wLyIWSVpBupzetIj4Q97qrBWNzF2AWTOIiJD0VWCrpFeBbmCGg9Vq5Zmr2SkkPUu63NzsiPhJ7nqsdXnN1axC0gzgEtIF0/+auRxrcZ65mgGSLgF+ASwn9YU6KyI+k7cqa2UOVxv2KjsEngW+FRFrK21XXiStufZkLc5alsPVhjVJ44FfAjsjYskp928GJkfEZdmKs5bmcDUzK4BPaJmZFcDhamZWAIermVkBHK5mZgVwuJqZFcDhamZWAIermVkBHK5mZgVwuJqZFeA/sOjDSnEj0XcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[2, 3], [1, -4], [1, 10]])\n",
    "b = np.array([7,3,-1])\n",
    "\n",
    "def ls_solution(A,b):\n",
    "    ATA = A.T @ A  \n",
    "    # Form the RHS vector:\n",
    "    print(A.T,b)\n",
    "    rhs = A.T @ b.T\n",
    "    pprint(rhs)\n",
    "    # solve the system\n",
    "    return sl.solve(ATA, rhs)\n",
    "\n",
    "# solve our 3x2 problem using LS\n",
    "ls_sol = ls_solution(A, b)\n",
    "# solve the three individual 2x2 problems - can still use the LS code\n",
    "ls_sol1 = ls_solution(A[ [0,1] ,: ], b[ [0,1] ])\n",
    "ls_sol2 = ls_solution(A[ [0,2] ,: ], b[ [0,2] ])\n",
    "ls_sol3 = ls_solution(A[ [1,2] ,: ], b[ [1,2] ])\n",
    "\n",
    "# plot this solution to see where it lies\n",
    "x = np.linspace(0,5,100)\n",
    "\n",
    "y1 = -(2./3.)*x + (7./3.)\n",
    "y2 = (1./4.)*x - (3./4.)\n",
    "y3 = -(1./10.)*x - (1./10.)\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "ax1.set_xlabel(\"$x$\", fontsize=14)\n",
    "ax1.set_ylabel(\"$y$\", fontsize=14)\n",
    "ax1.grid(True)\n",
    "\n",
    "ax1.plot(x,y1,'b', label='$2x+3y=7$')\n",
    "ax1.plot(x,y2,'r', label='$x-4y=3$')\n",
    "ax1.plot(x,y3,'g', label='$x+10y=-1$')\n",
    "ax1.plot(ls_sol1[0], ls_sol1[1], 'ko', label='solution 1')\n",
    "ax1.plot(ls_sol2[0], ls_sol2[1], 'ko', label='solution 2')\n",
    "ax1.plot(ls_sol3[0], ls_sol3[1], 'ko', label='solution 3')\n",
    "ax1.plot(ls_sol[0], ls_sol[1], 'bo', label='LS solution (all)')\n",
    "\n",
    "ax1.legend(loc='best', fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2  1  1]\n",
      " [ 3 -4 10]] [ 7  3 -1]\n",
      "array([16, -1])\n",
      "[[ 2  1  1]\n",
      " [ 3 -4 10]] [ 7  3 -1]\n",
      "array([16, -1])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x24e10fb0d48>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAFFCAYAAABL+O58AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5xcdX3v8dc7C7sxhpBsEiDkt5CHEVCjzU31crGUHwqUG7BiQ3hYoYWmtFK0xlYo99IrrX1A+6BWr1aMgGKVEIoXjCaURhS1VagRIxB+lAQSswbkR8KPdM0uyX7uH3MWJpOZ7Myc2TlnZt7Px2MfO+ec7/me7+wmnzn7Pd/P96uIwMzMmmtM1g0wM+tEDr5mZhlw8DUzy4CDr5lZBhx8zcwy4OBrZpYBB18zM0DSjZKekfRQ0b6/kvSApA2S/lXSkRXO3ZuU2SBpdVXX8zhfMzOQ9C5gF/CViDgu2TchIl5KXl8KHBMRF5c5d1dEjK/ler7zNTMDIuL7wI6SfS8Vbb4eaNjd6kGNqsjMrB1J+iTwQeBF4DcrFBsraT2wB7g6Iu4Ysd5O73Z4/aTumHjkOF7aM7buOgYHX/sMmzqmm2eHBmuuQ4Oq+/rDumq/bOW6dpf/dzFxUjcv7GzghXLm0Mnd7PivVxpS197u9HVEd7r/n93de159PWloLDvH7K7p/AkH1Va+1KSu/orHNj74ynMRMbXWOk84cWzs3DFUc1s2PvjKRqD4Da2IiBXFZSTNAb413O1QcuxyYGxE/GWZY0dGxHZJbwC+A5wcEZsP1J6Ov/OdeOQ4/vjW/8G6p+enqmdLX+Hf0PLxs7l219aaz+/5efr/qYdsbdwH6aGbB8ruX7J0LqtWPtmw6+TNORe8gRt+8IuG1PXy7PQfqAOz0n3QzZnx7Kuvz+tfwM3jNtR0/qlHPJrq+udMuL/isfmznqr9Pwqwc8cQX18zpebz5s96andELKznmombgTXAfsE3IrYn35+QdA/wNuCAwdd9vmZmFUiaV7S5GNjv00jSJEk9yespwPHAwyPV3fF3vmbtbEvf1H3ufq0ySSuBE4Epkvoo3OGeIemNwBCwFbg4KbsQuDgiLgLeBHxB0hCFG9qrI8LB18zqt+7p+am7HlpFRCwts/uGCmXXAxclr38IvLnW67nbwaxI10Dj+s0b0QffiGcBlk8OvmZmGXDwNTPLgIOvmY2a2156e9ZNyC0HXzOzDDj4mpWolGBi1kgOvmZtbjj70vLFwddsFDUy5dvai4OvWc5lPdY37bwnVp6Dr5lZBhx8zcwy4OBrVoZHPNhoc/A1s1HlRIvyHHzNzDLg4Gs2yjzczMrJXfCVdJqkxyRtknRZmeMXSHpW0obk66KiY+dLejz5Or+W63bKnKXWmZxokT+5mkxdUhfwOeBUoA/4saTVZWaFXxURl5Sc20th5vmFFJZ3/kly7s4mNN1sVPX8vDv1em6WL3m7810EbIqIJyJiELgFOKvKc98DrIuIHUnAXQecNkrttFHU/Wjfq19Z8oiH1zjRovHyFnynA9uKtvuSfaXeJ+kBSbdJmlnjuZZjpQE36wBsNlpy1e0AlFtru/RpxTeBlRExIOli4CbgpCrPLVxEWgYsA5h82CSO2nY2AFNfGVtns2FgfOFHeXhXD8vHz669gmNgzGC6z8Ixsxr3YKdrUfm6ent7WLJ0bsOuU0y7B6n0eRljm5NiW/r+9vakX/4dYKg7fT1D3UOpzu8d6ua8/gV1nz9hW/13v1u63l1m74frrq8d5C349gEzi7ZnANuLC0TE80WbXwSuKTr3xJJz7yl3kYhYAawAmH7sxNg88w4g5Z9W4woPNZaPn821u7bWVUXaHP5GPlWv9Cf3kqVzWbXyyYZdp9hId7mD82eMynWLlb6/F4/qaUi9L89OH3zT9vn+Re90bh63oe7z0zyUPmfC/XWf267y1u3wY2CepLmSuoFzgdXFBSRNK9pcDDySvL4LeLekSZImAe9O9lmNsujrrKZ7oZW7IDp9uJkTLfaXqzvfiNgj6RIKQbMLuDEiNkq6ClgfEauBSyUtBvYAO4ALknN3SPorCgEc4KqI2NH0N2FmVoVcBV+AiFgLrC3Zd2XR68uByyuceyNw46g20Bquljva7kf7mtL9MOzQzQMN63owK5a3bgczqyDtM4GBwdzda3U0B19rOa3c92v5JelGSc9Ieqho399JejQZ2nq7pIkVzj1gZm45Dr6WKQfS1tEBiRZfZv/ErHXAcRHxFuA/KdPlWZSZezpwDLBU0jEjXczB11qSg7Y1WkR8n8JD/OJ9/xoRe5LNeykMYS1VV2aug6+ZWXV+H7izzP66smvdA28tq1kjHxo14uGQrdGQZItOt3PvuDrHDa+ZIml90Y4VScLViCRdQWF469fKHS6zb8SB3Q6+beTl2WqpwfzuOrAmey4iFtZ6UjI97ZnAyRFR7j/YiJm55bjbwcyaohWz3CSdBnwcWBwR/RWKjZiZW46Dr7W0Trt7TjvW1yqTtBL4EfBGSX2SLgQ+CxwCrEsWb7guKXukpLVQyMwFhjNzHwFujYiNI13P3Q5mHWRL31TmzHg262bkUkQsLbP7hgpltwNnFG3vl5k7Et/5FvFSQgVOpzUbfQ6+lolGdhd0WteDtQcHX7MqeEmhgg7IcmsaB1+zJmqloYA2uhx8rS2468FajYOvmVkGHHzNzDLg4NtAHj9ZHXcRpNPKiRatmOU2Whx8rW2MdlBvlxEPW/qmZt0Ew8E3V9IuDW5mrcPB18wsAw6+1lZaoT/ZY30NHHzNrEbOcmsMB18zsww4+JqZZSB3wVfSaZIek7RJ0mVljn9U0sOSHpB0t6TZRcf2JhMeb5A04kzy1p5Gs9+3XYabWfZyFXwldQGfA04HjgGWSjqmpNhPgYUR8RbgNuBvi479KiIWJF+Lm9Joswy0cqKFFeQq+AKLgE0R8UREDAK3AGcVF4iI7xatpXQvhcXqrEW0wmiETpBlooWz3AryFnynA9uKtvuSfZVcCNxZtD1W0npJ90o6u54GeDULM2uGvK3hpjL7yg6KlPQBYCHwG0W7Z0XEdklvAL4j6cGI2Fzm3GXAMoDJh03iqG37xumpr4yts/kwsauH5eNnj1ywgjFHp/s8HDOrMWNIuxaVr6e3t4clS+fWXa92H+iztLFibO1/mlfz/vb2lPtnWruh7nT1DHUP1VT+8JJ/mz399f8uJmxLN9xsS9e7gQ+nqqPV5S349gEzi7ZnANtLC0k6BbgC+I2IePUJSLKoHRHxhKR7gLcB+wXfiFgBrACYfuzE2Dzzjn2OpxnH+L4XFnLtrq11n5+2L69RA/grPVhasnQuq1Y+WXe9zex2GJxfe49UNe+vUWvcvTw7XfCtNR19+fjZ+/zbTDMRVNq/EM+ZcH+q89tB3rodfgzMkzRXUjdwLrDPqAVJbwO+ACyOiGeK9k+S1JO8ngIcDzzctJabmdUgV8E3IvYAlwB3AY8At0bERklXSRoevfB3wHjgn0uGlL0JWC/pZ8B3gasjwsHXbBQ4yy29vHU7EBFrgbUl+64sen1KhfN+CLx5dFtnVuiSaVTXg3WuXN35Wnpp+xHbiYe1WZ45+Jq1KCdaNJakGyU9I+mhon29ktZJejz5PqnCuTVn1zr4NlhP956sm2Bm9fkycFrJvsuAuyNiHnB3sl1Ozdm1Dr5mHcrLCe0rIr4P7CjZfRZwU/L6JqCu5K1yHHzNMtLJk6q3UIrx4RHxFEDy/bAK5WrOrs3daIdONzBr0H15DdT9aF9dyRaWXy/tGVvnULc1UyStL9qxIkm4aoSqsmuLOfia1cHDzVrScxGxsMZzfilpWkQ8JWka8Ey5QtVm1xZzt4OZWWWrgfOT1+cD3ygtUG92rYOvmdWl3bLcJK0EfgS8UVKfpAuBq4FTJT0OnJpsI2mhpOuTU+vKrnW3g5kZEBFLKxw6uUzZ9cBFyeu6smt951uG5/RtL850szxy8DVrYR4Z07ocfM3MMuDga5ahrBMtnOWWHQdfM7MMOPia1anSUktm1XDwtabJctSBRzxY3jj4mpllwMF3FKRZFTYvPG+B2ehy8DWzurVbinEzOfi2Ia/jZpZ/Dr7WMdr1oZuz3FqTg69ZxrJOtLBsOPjm0MCswaybYFXyWF+rl4OvWYdzinE2chl8JZ0m6TFJmyTtt1SzpB5Jq5Lj90maU3Ts8mT/Y5Le08x2m5lVK3fBV1IX8DngdOAYYKmkY0qKXQjsjIijgU8B1yTnHgOcCxwLnAb8Y1JfzTynr5mNptwFX2ARsCkinoiIQeAW4KySMmcBNyWvbwNOlqRk/y0RMRARTwKbkvrMzHIlj8sITQe2FW33Ab9eqUxE7JH0IjA52X9vybnTSy8gaRmwDGDyYZM4atvZZRsy9ZWxNTe+d2gc5/UvYGB8uh/tmKPTfS6OmZX+CXrXov3r6O3tYcnSuXXVp937/SoyEWMrD82q5/3t7Uk/rnqoO10dQ91DI5Y5vKuH5eNnlz3W01//72bCtnoTLdbUfc12kMfgW+5fYWkUqFSmmnOJiBXACoDpx06MzTPvKNuQerJ3zutfwM3jNrBlR7qHGGnHbjZi+FK5J/lLls5l1con66ovL+NsB+fPqHisnvfXiFTstIkx1YyQWT5+Ntfu2lr2WNqUeHfT1S6P3Q59wMyi7RnA9kplJB0EHArsqPJcy8iBgp5Zp8lj8P0xME/SXEndFB6grS4psxo4P3l9DvCdiIhk/7nJaIi5wDzgP5rUbrPMOMut9eQu+EbEHuAS4C7gEeDWiNgo6SpJi5NiNwCTJW0CPgpclpy7EbgVeBj4F+BDEbG32e/BOksjEi2c5dZ58tjnS0SsBdaW7Luy6PVu4P0Vzv0k8MlRbaCZWUq5u/NtF+0wp6+ZjR4HX+s4eRl1kSdOMQZJH5b0kKSNkj5S5rgkfSbJoH1A0tvTXM/B18w6nqTjgD+gkJT1VuBMSfNKip1O4SH+PAp5Ap9Pc00HXzMzeBNwb0T0Jw/9vwe8t6TMWcBXouBeYKKkafVe0MG3TXk1C7OaPAS8S9JkSeOAM9g3ZwDKZ9/WnRqYy9EOVshYynrs5otH9Xi+WsudwcGD6u2jniJpfdH2iiTblYh4RNI1wDpgF/AzYE/J+VVl0FbLwfcATj3iUS8QaNY+nouIhZUORsQNFHIIkPQ3FO5sizU0g9bdDmaWWjvcpEg6LPk+C/htYGVJkdXAB5NRD+8AXoyIp+q9nu98zRrg0M0DDZlgJ42en3d7Cap0vi5pMvAKhezYnZIuBoiI6ygkfp1BYarafuD30lzMwdc6Uvejfbmb6OeQreEHpRmKiBPK7Luu6HUAH2rU9dztYGaWAQdfM7MMOPiaGeAU42Zz8DUzy4CD7yjyzGZmVomDr5lZBhx8zcwy4OBrHcvz+lqWHHzNGqTTJyFqhxTjZnLwzbG0qaLOluo8Wc+EZ9Vz8B3BqUc8mnUTrIN4FePO4eBrB5T1ZDFm7crB18wsAw6+ZmYZcPC1psrbNI62L8/v0Dy5Cb6SeiWtk/R48n1SmTILJP1I0kZJD0haUnTsy5KelLQh+VrQ3HdgZla93ARf4DLg7oiYB9ydbJfqBz4YEccCpwH/IGli0fE/i4gFydeG0W+ymVl98hR8zwJuSl7fBJxdWiAi/jMiHk9ebweeAfx3ktXNWW6WlTwtI3T48GJ0EfHU8GJ2lUhaBHQDm4t2f1LSlSR3zhFRNuVI0jJgGcDkwyZx1Lb94vw+pr4ytuo30Ts0jvP6i3o8emFgsP4f85ij030+jpmVftxo16LX6ujt7WHJ0rmp6tPu6Wmb1FAx9rXEhLTvb29P+sSWoe50dQx1D5Xdf3hXD8vHzx7x/J7++n8/E7bVkuW2pu7rtIOmBl9J3waOKHPoihrrmQb8E3B+RAz/S7sceJpCQF4BfBy4qtz5EbEiKcP0YyfG5pl3HPB6taRNnte/gJvH7dvjsWVH/TfnaTOWGjFovzhtdsnSuaxa+WSq+vJ2t1n8ELAR7y/t2Oi0mYmVMiOXj5/Ntbu2jnh+mqlQnZRUvaYG34g4pdIxSb+UNC25651GoUuhXLkJFD4y/1dE3FtU9/ASzgOSvgR8rIFNNzNrqDz1+a4Gzk9enw98o7SApG7gduArEfHPJcemJd9Fob/4oVFtrdkoSfvXSpbzO3hynerlKfheDZwq6XHg1GQbSQslXZ+U+R3gXcAFZYaUfU3Sg8CDwBTgr5vbfDOz6uXmgVtEPA+cXGb/euCi5PVXga9WOP+k0WrbqUc8mtkn+sCswVR3Mi/PlidrMcuhPN35mpl1jBGDryR34nQ4z2xmnUDSG4u6MzdIeknSR0rKnCjpxaIyV9Z7vWrufH8q6dPl0n3NsjRp739xzfO3M2lvf9ZNaSudOr9DRDw2nCEL/BqFjNrbyxT9QVEmbdnhrNWoJvguAo4FHpf0J5K66r2YWSOdt2s9x73yFOft+nHWTbH2czKwOSJGHhhdpxEfuEXEg8Apks4G/g74I0nLI+LO0WqUtbfB+TNSJVrc8fQX6GHvq9tn/mojZ/5qIwN0cfYRf1hzfd2P9nm2NSt1LrCywrF3SvoZsB34WERsrOcCVY92iIg7JK0F/hS4RdK/Ax+NCKe0WFP93tQPcNHLP+S/736SsexhNwfxw7Fzuf6Q47NumjWBBlXvCKApktYXba9Isl33rb+QT7CYQtZsqfuB2RGxS9IZwB3AvHoaU+toh3HATyhMfPMe4AFJn5F0aD0XN6vHzq7X069uutnDAF10s4d+dbOza1zWTbN8ey4iFhZ97Rd4E6cD90fEL0sPRMRLEbEreb0WOFjSlHoaU81oh49I+pqk/wSeB74J/Dfg0xTG374ReFjSr9fTgE6QJlfeyps01M+a1x3Ln05+H2tedyyThvLz0K3Tl5BvA0up0OUg6Ygki3Z4cq8xFOJizarpdlgO/Aj4PHAv8JOIKJ654yuSPg7cSOHBnNmo++tJp7/6+h8P/Y0MW2LtRNI4Chm2f1i072KAiLgOOIfCc689wK+AcyOiriymah64zayini8Bf1NPA8ysvax7en7Lzm4WEf3A5JJ91xW9/izw2UZcq1EZbs8Co5bea9ZpWnlyHatOQ4JvFHyvEXXZ/irNz1qttPPDmlnjeW6HKrXqn1Fmlk8OvlYVz+9g1lgOvmZmGXDwNbP9dOrkOs3k4GuZ8FwK1ukcfM3MMuDga0b+lrO39ufg2ySe38HMijn4mpllwMHXbBR4ZjMbiYOvmTXcuqe97u5IHHxbhOd36DyeXKe9OfjWwPM7mFmj5Cr4SuqVtE7S48n3ssvVS9oraUPytbpo/1xJ9yXnr0rWYrIG8fwOZo2Tq+ALXAbcHRHzgLuT7XJ+FRELkq/FRfuvAT6VnL8TuHB0m2tpOMvNOlnegu9ZFBbnJPl+drUnJusqnQTcVs/5ZmbNVPXS8U1yeEQ8BRART0k6rEK5sckS0HuAqyPiDgpLf7wQEXuSMn3A9HInS1oGLAOYfNgkjtpWfYye+srYAx7vHRrHef0Lyh4bGJ/uxz3m6HSflWNmpXuAA9D7+m6WLJ2bup5h2l32V5SJGNtNb29Pw97f3p70DzmHutPVMdQ99Orrw7t6WD5+dvUnvzCbnu49I5erYMK2kUY8rKm77nbQ9OAr6dvAEWUOXVFDNbMiYrukNwDfkfQg8FKZcmWjTbJk9AqA6cdOjM0z76j6wiMNoTmvfwE3j9tQ9tiWHelmikr79Drt03OAZYuOZNXKJ1PXMyxPab2D82ewZOnchr2/RvSRpx2lUjxKZvn42Vy7a2tN56fJzPQD6gNrevCNiFMqHZP0S0nTkrveacAzFerYnnx/QtI9wNuArwMTJR2U3P3OALY3/A2YmTVA3vp8VwPnJ6/PB75RWkDSJEk9yespwPHAw8nyzd+lsLRzxfOz5PkdzGxY3oLv1cCpkh4HTk22kbRQ0vVJmTcB6yX9jEKwvToiHk6OfRz4qKRNFPqAb2hq683MqpSrB24R8Txwcpn964GLktc/BN5c4fwngEWj2UYzs0bI252vjaI8phh7rK91KgffFpJ2fodGaMTwKbM8kjRR0m2SHpX0iKR3lhyXpM9I2iTpAUlvT3M9B98aefiMWXVacGazTwP/EhHzgbcCj5QcPx2Yl3wtAz6f5mIOvmajxHP6tg5JE4B3kTykj4jBiHihpNhZwFei4F4KQ1un1XtNB1+zNuZpJav2BuBZ4EuSfirpekmvLykzHdhWtF0xi7YauRrtYGb7OmRr5PJBaZa6BuvO1pySTEswbEWS7QqFWPh24E8i4j5Jn6Ywsdf/Lipf7hdRd9qog6+ZdYrnImJhhWN9QF9E3Jds38b+syr2ATOLtlNl0brbwcw6XkQ8DWyT9MZk18nAwyXFVgMfTEY9vAN4cXgisHo4+DaZU4z3l5exvnma5Mcy8SfA1yQ9ACwA/kbSxZIuTo6vBZ4ANgFfBP44zcXc7WBmFW3pm9oxNwwRsQEo7Za4ruh4AB9q1PV852tmlgEH3xaTh1WMvZabWXoOvmZmGXDwrYNTjM0sLQdfM7MMOPhaLuRluJlZszj4mtmoacGZzZrGwdfMLAMOvmZmGXDwzUCnZAxZPnhayXxy8LW6ONHCLB0H3w7k+WHNsufg24LysJDmaPBws/LqnDjccs7B12wUeR03q8TBt05OMbZOsaVvatZNaEu5Cb6SeiWtk/R48n1SmTK/KWlD0dduSWcnx74s6cmiYwua/y7MzKqTm+BLYb2kuyNiHnA3+6+fRER8NyIWRMQC4CSgH/jXoiJ/Nnw8mRjZzCyX8hR8zwJuSl7fBJw9QvlzgDsjon9UW2VmNgrytIzQ4cOL0UXEU5IOG6H8ucDfl+z7pKQrSe6cI6Ls0w5Jy4BlAJMPm8RR20aK8+VNfWXsfvt6h8ZxXv/IPR4D49P96Mccne5zc8ys+p6gTxnfzYUnTC9snABdA41+Ej8X7c5uNMek3h6WLJ3b0Dr39qQf2jfUna6Ooe4hDu/qYfn42XWd39M/ve5rT9hWaX6HNXXX2Q6aGnwlfRs4osyhK2qsZxrwZuCuot2XA08D3cAK4OPAVeXOj4gVSRmmHzsxNs+8o5bLv6rcpCHn9S/g5nFV9HiMS/cgI23WUr3Dly48YTo3/OAXr26PxtP8LBeyfO8nTmDVyicbWmcjElLSjs0emDXI8vGzuXbX1rrOT5OV6YfT5TU1+EbEKZWOSfqlpGnJXe804JkDVPU7wO0R8UpR3cNLOA9I+hLwsYY0uk29PFseP2qWoTz1+a4Gzk9enw984wBllwIri3ckARtJotBf/NAotNHMrCHyFHyvBk6V9DhwarKNpIWSrh8uJGkOMBP4Xsn5X5P0IPAgMAX46ya0OTPtmuVm7cdz+paXmwduEfE8cHKZ/euBi4q2twD79f5HxEmj2T5rnsH5MzLp93V6szVTnu58zcw6hoNvCn6Ka9ZeJHVJ+qmkb5U5doGkZ4uyaC8qV0e1ctPtYGaWAx8GHgEmVDi+KiIuacSFfOdrqXhS9dbg1SxGJmkG8FvA9SOVbQQHXzOzgn8A/hwYOkCZ90l6QNJtkmamuVjHdztM6urcqSHynGiR1YgHy7+u3VFvZuUUSeuLtlck2a5IOhN4JiJ+IunECud/E1gZEQOSLqYwB03do6w6Pvhmac6MZz1XqlnzPBcRCyscOx5YLOkMYCwwQdJXI+IDwwWS4bDDvghck6Yx7nYwsxG1+01CRFweETMiYg6FSbu+Uxx44bUs2sRiCg/m6ubg28Kc5dY58to91O4kXSVpcbJ5qaSNkn4GXApckKZudztYbrnf17IQEfcA9ySvryzafzmF2RMbwne+lpqHm5nVzsE3JWe5mVk9HHzNzDLg4GtmlgEH3w6Xdnma0eZpHtuD5/Tdn4MvcM6E+zO7dpq1scysdTn4mpllwMG3xTnRwqw1OfhaQ4zmWN9m9Pu6b9mazcHXzCwDDr5mZhlw8G0AZ7mZWa0cfC33Y33BfbKNMGbQ/93zxL8NM7MMOPiajbI6l7yxNper4Cvp/clkxUOSKi33gaTTJD0maZOky4r2z5V0n6THJa2S1BJLtrZLlpunlmxv7b6aRbPlKvgCDwG/DXy/UgFJXcDngNOBY4Clko5JDl8DfCoi5gE7gQurvXCWKcZpdUqihft9rZ3kKvhGxCMR8dgIxRYBmyLiiYgYBG4BzpIkCiuJ3paUuwk4e/Raa+3CQd2y0IrLCE0HthVt9wG/DkwGXoiIPUX7p5erQNIyYBnAYYcfypZHCz0XR+0dV3ejpr4yFoDeoXGc17+g5vMHxqf7VYw5Ot3n6JhZ1a0RNmV8NxeeUPbHCkDXotFda0y7K1+7XjH2td6p3t4eliyd2/Br7O1JP6JkqDtdHYf1HMzy7tmp6ujpr//nP2Fb6cxma1K1pdU1PfhK+jZwRJlDV0TEN6qposy+OMD+/XdGrABWABz3lu6YM/9qANa/9PYqLl/e8JR55/Uv4OZxG2o+f8uOdP1pPT9P171d7QKNF54wnRt+8IuKx5vxcKnR67oV3/kuWTqXVSufbGj90Jj+8LRDAi89ejrXDm5NVUea5xMeD7+vpgffiDglZRV9wMyi7RnAduA5YKKkg5K73+H9VoWXZ8sr5Jo1Ua76fKv0Y2BeMrKhGzgXWB0RAXwXOCcpdz5QzZ10Q/hT3cxqkavgK+m9kvqAdwJrJN2V7D9S0lqA5K72EuAu4BHg1ojYmFTxceCjkjZR6AO+odnvodM1Y7iZH5C1Jq9msa9cPXCLiNuB28vs3w6cUbS9FlhbptwTFEZDmJnlWq7ufDtZuyRatJpWuot2n3x7cfBtE52SaDGslYKm5Z+ksZL+Q9LPkizbT5Qp05Nkzm5KMmnnpLmmg2+RVs5yM7NUBoCTIuKtwALgNEnvKClzIbAzIo4GPkUho7ZuDr5m1vGiYFeyeXDyVdrPcxaFzFkoZNKenGTW1sXB117VqHl9mzXBjrserJEkdUnaADwDrIuI+0qKvJpdm4y6epHCqKq65Gq0g1kzOXi3Ju0erDfLcYqk9UXbK5JsVwAiYi+wQNJE4HZJx0XEQ8WXLlNn3fNU8BwAAAgeSURBVE9BHXzNrFM8FxEVp6odFhEvSLoHOI3CTIvDhrNr+yQdBBwK7Ki3Me52aCBnuTWf716tESRNTe54kfQ64BSg9D/0agqZs1DIpP1OkllbF9/5Wkdy0LYS04CbkvnCx1DInP2WpKuA9RGxmkLG7D8lGbQ7KExtUDcH3xyZM+PZVKsFDMwaTD27WSsanD+j4TOdWWeJiAeAt5XZf2XR693A+xt1TXc7mJllwMG3hBMtGqPZ67nV0o3gLof6eR23xnHwtX00aqxvFhxUR9aJ3VJ55eBrHcUB2vLCwdfaioOrtQoHX2s7lQKwA7PliYNvg004eHfWTTD2D7QOvJY3Dr5tJk/z+jZ7xEOpwfkzXv2yfPBSQq9x8M0Zr2hh1hkcfG0/rTzczKxVOPiamWXAwbcMZ7mZ2Whz8LVRlfVDtzzwz8DKcfA1M8uAg6+ZWQZyE3wlvV/SRklDksou9SFppqTvSnokKfvhomP/R9IvJG1Ivs5oXuv35RUtzGwkuQm+FNZK+m3g+wcoswdYHhFvAt4BfEjSMUXHPxURC5KvtaPY1lxrRKKFh5uZja7cBN+IeCQiHhuhzFMRcX/y+mXgEQrLObcVJ1qYtb+WXUZI0hwKy37cV7T7EkkfBNZTuEPeWeHcZcAygMMOP5Qtj162X5mj9o6rq109gxM5atvZTH1lbF3nDxsYn+5XM+bo9J+rY2btvzbglPHdXHhCjZ93J0DXQN3rDDZVb28PS5bObWide3sa91fEUHf9dR3WczCXHj2doe6hVG3o6U93vzNh23CK8ZpU9bS6pgZfSd8Gjihz6IqI+EYN9YwHvg58JCJeSnZ/HvgrIJLv1wK/X+78iFgBrAA47i3dMWf+1fuVWf/S26ttzj6O2nY2m2fekTqHfcuOdCsGNGLS7EO27h8wLzxhOjf84Bc113Xo5oHU7WmGJUvnsmrlkw2ts5FDzdJ0B1169HQ+s+kXqbul0v5l5mciBU0NvhFxSto6JB1MIfB+LSL+X1Hdvywq80XgW2muc86E+7mtzgBsZjaS3PT5VkOSKCzf/EhE/H3JsWlFm++l8ADPzCyXchN8Jb1XUh/wTmCNpLuS/UdKGh65cDzwu8BJZYaU/a2kByU9APwm8KfNfg9mZtXKzQO3iLgduL3M/u3AGcnrfwPKdnpFxO+OagNbzMCswdT9vi/PVtl+XzNLLzd3vu3GDxX25fkNzPbl4JtTHutrebWlL91InLySdKOkZySVfV4k6URJLxZ1eV6Z5nq56XYwM8vYl4HPAl85QJkfRMSZjbiY73zNzICI+D6wo1nXc/A1M6veOyX9TNKdko5NU5EiOvtptqRnga0NrHIK8FwD68sbv7/Wlbf3Njsiau5AlvQvFN5LrcYCu4u2VyTZrsV1zwG+FRHHlbnuBGAoInYlQ1w/HRHz6mhHob5OD76NJml9RJSdErMd+P21rnZ+b41yoOBbpuwWYGFE1PWB5m4HM7MqSDoiybJF0iIK8fP5euvzaAczM0DSSuBEYEqSbfuXwMEAEXEdcA7wR5L2AL8Czo0UXQcOvo23YuQiLc3vr3W183tLLSKWjnD8sxSGojWE+3zNzDLgPl8zsww4+NbhQAt5FpVpaCpiM0kaK+k/kvGMGyV9okyZHkmrJG2SdF/ylDj3qnxvF0h6tuh3d1EWbU1DUpekn0rab17rVv3dtRv3+dZneCHP+yUdAvxE0rqIeLikXMNSEZtsADgpGc94MPBvku6MiHuLylwI7IyIoyWdC1wDLMmisTWq5r0BrIqISzJoX6N8mMIahxPKHGvV311b8Z1vHdp9Ic8o2JVsHpx8lT4cOAu4KXl9G3Dy8DCcPKvyvbU0STOA3wKur1CkJX937cbBN6UKC3kOa1gqYrMlf7ZuAJ4B1kVE6fubDmwDiIg9wIvA5Oa2sj5VvDeA90l6QNJtkmY2uYlp/QPw50CllTJb9nfXThx8U6iwkOew+ymkT74V+L/AHc1uXxoRsTciFgAzgEWSSjN+yt0ptcQdZBXv7ZvAnIh4C/BtXrtLzD1JZwLPRMRPDlSszL6W+N21EwffOlVayHNYRLw0/OdtRKwFDpZUTz56piLiBeAe4LSSQ33ATABJBwGH0sQZoRqh0nuLiOcjYni55S8Cv9bkpqVxPLA4SX29hcKSW18tKdPyv7t24OBbhwMt5FlUpqGpiM0kaaqkicnr1wGnAKVLc6wGzk9enwN8J022T7NU895KFmNdTKFPvyVExOURMSMi5gDnUvi9fKCkWEv+7tqNRzvUZ3ghzweTvkOAvwBmweikIjbZNOAmSV0UPjRujYhvSboKWB8Rqyl8+PyTpE0U7prOza65NanmvV0qaTGFUS07gAsya22DtMnvrq04w83MLAPudjAzy4CDr5lZBhx8zcwy4OBrZpYBB18zsww4+JqZZcDB18wsAw6+ZmYZcPC1liLp/ZIGJM0u2vdpSZslHZ5l28xq4Qw3aynJfBk/Bn4aEX8g6WMUpk88PiIez7Z1ZtXz3A7WUiIiJP0FsEbSZuAKCitTOPBaS/Gdr7UkST8EFgH/MyLuzLo9ZrVyn6+1HEknAW+lMCn4LzNujlldfOdrLUXSW4HvAR+lsE7Z+Ih4T7atMqudg6+1jGSEww+BL0TEVcnyPw9Q6PO9J9PGmdXIwddagqRe4N+B70fEHxbtXwXMioh3ZtY4szo4+JqZZcAP3MzMMuDga2aWAQdfM7MMOPiamWXAwdfMLAMOvmZmGXDwNTPLgIOvmVkGHHzNzDLw/wEhhyw6b29bjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[2, 3], [1, -4], [1, 10]])\n",
    "b = np.array([7,3,-1])\n",
    "\n",
    "# solve our 3x2 problem using LS\n",
    "ls_sol = ls_solution(A, b)\n",
    "\n",
    "# construct a grid of points in (x,y) around this solution location \n",
    "pert = 1.\n",
    "# construct a mesh\n",
    "# 100 points spread either side of ls_sol, of extent pert\n",
    "x = np.linspace(ls_sol[0] - pert, ls_sol[0] + pert, 100)\n",
    "y = np.linspace(ls_sol[1] - pert, ls_sol[1] + pert, 100)\n",
    "errors = np.zeros([len(x), len(y)])\n",
    "for i, xi in enumerate(x):\n",
    "    for j, yj in enumerate(y):\n",
    "        errors[i,j] = np.linalg.norm( A@np.array([xi, yj]) - b)\n",
    "\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "ax1.set_xlabel(\"$x$\", fontsize=14)\n",
    "ax1.set_ylabel(\"$y$\", fontsize=14)\n",
    "ax1.grid(True)\n",
    "\n",
    "cs = ax1.contourf(x,y,errors)\n",
    "fig.colorbar(cs, ax=ax1)\n",
    "# add our LS solution\n",
    "ls_sol = ls_solution(A, b)\n",
    "\n",
    "ax1.plot(ls_sol[0], ls_sol[1], 'r*')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - An even simpler over-determined case\n",
    "\n",
    "In the lecture, and above, we considered the simple case of three equations, two unknowns.\n",
    "\n",
    "Of course there is an even simpler case - two-equations, one unknown.\n",
    "\n",
    "An example might be\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "  2x &= 8 \\\\[5pt]\n",
    "   3x &= 9\n",
    "\\end{align*}\n",
    "   \\quad \\iff \\quad\n",
    "  \\begin{pmatrix}\n",
    "    2 \\\\[5pt]\n",
    "    3\n",
    "  \\end{pmatrix}\n",
    "  \\begin{pmatrix}\n",
    "    x \n",
    "  \\end{pmatrix}=\n",
    "  \\begin{pmatrix}\n",
    "    8 \\\\\n",
    "    9 \n",
    "  \\end{pmatrix}   \n",
    "$$\n",
    "\n",
    "Does this have a solution?\n",
    "\n",
    "No clearly not. What solution does the least square approach return?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "Let's do it by hand and then with code.\n",
    "\n",
    "The least squares solution will minimise the quantity $\\| A\\boldsymbol{x} - \\boldsymbol{b}\\|_2$, \n",
    "\n",
    "i.e. here (squaring for simplicity), that minimises\n",
    "\n",
    "$$\n",
    "\\left\\| \n",
    "\\begin{align*}\n",
    "  2x - 8 \\\\[5pt]\n",
    "   3x - 9\n",
    "\\end{align*}\n",
    "\\right\\|_2^2\n",
    " = { (2x - 8)^2 + (3x - 9)^2 }\n",
    "$$\n",
    "\n",
    "This is equal to \n",
    "\n",
    "$$ (4 x^2  -32 x + 64) + (9 x^2  -54 x + 81) = 13 x^2  -86 x +  145 $$\n",
    "\n",
    "to find the $x$ that minimises this let's differentiate and set to zero:\n",
    "\n",
    "$$26 x - 86 = 0 \\Rightarrow x = 86/26 = 43/13$$\n",
    "\n",
    "<br>\n",
    "\n",
    "Now let's compute by hand what the least squares implementation from the lecture does, i.e. solve \n",
    "\n",
    "$$A^TA\\boldsymbol{x} = A^T\\boldsymbol{b}$$\n",
    "\n",
    "Here \n",
    "\n",
    "$$ A^TA = (2, 3) \n",
    "  \\begin{pmatrix}\n",
    "  2 \\\\[5pt]\n",
    "   3\n",
    "  \\end{pmatrix}\n",
    " = 13\n",
    " $$\n",
    " \n",
    " and\n",
    " \n",
    "$$A^T\\boldsymbol{b} = (2,3)\n",
    "  \\begin{pmatrix}\n",
    "  8 \\\\[5pt]\n",
    "   9\n",
    "  \\end{pmatrix}\n",
    "= 43$$\n",
    "\n",
    "and so this also produces $ x = 43/13$.\n",
    "\n",
    "Let's implement in code and check we get the same answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3.30769231]] 3.3076923076923075\n"
     ]
    }
   ],
   "source": [
    "A = np.array([ [2], [3] ])\n",
    "b = np.array([ [8], [9] ])\n",
    "# Form the matrix A.T @ A\n",
    "ATA = A.T @ A  \n",
    "\n",
    "# Form the RHS vector:\n",
    "rhs = A.T @ b\n",
    "\n",
    "# solve the system\n",
    "ls_sol = sl.solve(ATA, rhs)\n",
    "\n",
    "print(ls_sol, 43/13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - Outer-product\n",
    "\n",
    "Compute the outer-product ($\\boldsymbol{a}\\boldsymbol{b}^T$) of the column vectors\n",
    "\n",
    "$$\\boldsymbol{a} = \n",
    "\\begin{pmatrix}\n",
    "1 \\\\\n",
    "2\\\\\n",
    "3 \n",
    "\\end{pmatrix}, \n",
    "\\qquad\n",
    "\\boldsymbol{b} = \n",
    "\\begin{pmatrix}\n",
    "4 \\\\\n",
    "5\\\\\n",
    "6 \n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "What is the resulting matrices rank?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution\n",
    "$$\n",
    "\\boldsymbol{a}\\boldsymbol{b}^T\n",
    "= \n",
    "\\begin{pmatrix}\n",
    "1 \\\\\n",
    "2\\\\\n",
    "3 \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "4, 5, 6\n",
    "\\end{pmatrix}=\n",
    "\\begin{pmatrix}\n",
    "4 & 5 & 6  \\\\\n",
    "8 & 10 & 12  \\\\\n",
    "12 & 15 & 18  \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "which with row operations we could trivially arrive at\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "4 & 5 & 6  \\\\\n",
    "0 & 0 & 0  \\\\\n",
    "0 & 0 & 0  \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "so there is only a single linearly independent row (equivalently only a single linearly independent column).  The rank is therefore indeed 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - Matrix rank and RREF (a non-square case)\n",
    "\n",
    "Consider the rectangular matrix\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "3 & 1 & 9 & 4 \\\\\n",
    "2 & 1 & 7 & 3 \\\\\n",
    "5 & 2 & 16 & 7 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "convert to REF and RREF. \n",
    "\n",
    "From these what is the rank of this matrix?\n",
    "\n",
    "What is the null space?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution \n",
    "\n",
    "\n",
    "The RREF is\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "1 & 0 & 2 & 1 \\\\\n",
    "0 & 1 & 3 & 1 \\\\\n",
    "0 & 0 & 0 & 0 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "This is formed of two linearly independent columns  (and equivalently  of  only  two independent rows) so its rank is 2. This is less that the minimum of $m$ and $n$ and so the matrix is not full rank.\n",
    "\n",
    "Now consider the augmented matrix representing three linear equations in four unknowns (the RHS vector is all zero)\n",
    "\n",
    "$$\n",
    "\\left(\n",
    "  \\begin{array}{cccc|c}\n",
    "3 & 1 & 9 & 4 & 0 \\\\\n",
    "2 & 1 & 7 & 3 & 0 \\\\\n",
    "5 & 2 & 16 & 7 & 0 \n",
    "  \\end{array}\n",
    "\\right)$$\n",
    "\n",
    "The RREF is\n",
    "\n",
    "$$\n",
    "\\left(\n",
    "  \\begin{array}{cccc|c}\n",
    "1 & 0 & 2 & 1 & 0 \\\\\n",
    "0 & 1 & 3 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 0 & 0 \n",
    "  \\end{array}\n",
    "\\right)$$\n",
    "\n",
    "so any vector for which $v_1 + 2v_3 + v_4=0$ and $v_2 + 3v_3 + v_4=0$\n",
    "will be a solution of $A\\boldsymbol{v} = \\boldsymbol{0}$ \n",
    "and thus will lie in the null space of $A$.\n",
    "\n",
    "If we chose arbitrary values for the variables that appear more than once, say $v_3 = \\alpha$ and $v_4 = \\beta$, then we obtain $v_1 = -2\\alpha - \\beta$ and $v_2=-3\\alpha-\\beta$.\n",
    "\n",
    "Therefore note that the solution for $v$ in the null space can be written as\n",
    "\n",
    "$$\\boldsymbol{v} = \\alpha \n",
    "\\begin{pmatrix}\n",
    "-2 \\\\\n",
    "-3\\\\\n",
    "1\\\\\n",
    "0\n",
    "\\end{pmatrix} \n",
    "+\\beta\n",
    "\\begin{pmatrix}\n",
    "-1 \\\\\n",
    "-1\\\\\n",
    "0\\\\\n",
    "1\n",
    "\\end{pmatrix} \n",
    "$$\n",
    "\n",
    "That is, any vector in the null space of $A$ can be written as a linear combination of the two \n",
    "vectors above.  This null space is a two-dimensional plane within $\\mathbb{R}^4$.\n",
    "The null space thus  forms a sub-space of  $\\mathbb{R}^4$.\n",
    "\n",
    "Note that the number of independent vectors that must be linearly combined to form the null \n",
    "space is equal to the number of non-pivot columns in the RREF.   \n",
    "\n",
    "Now consider the problem $A\\boldsymbol{x}=\\boldsymbol{b}$ where\n",
    "\n",
    "$$\\boldsymbol{b} = \\begin{pmatrix}\n",
    "22 \\\\\n",
    "17\\\\\n",
    "39\n",
    "\\end{pmatrix} $$\n",
    "\n",
    "One particular solution to this is \n",
    "\n",
    "$$\\boldsymbol{x}_{\\text{part}} = \\begin{pmatrix}\n",
    "1 \\\\\n",
    "2\\\\\n",
    "1 \\\\\n",
    "2\n",
    "\\end{pmatrix} $$\n",
    "\n",
    "We can add to this solution any vector from the null space and by linearity it will be another solution, e.g.\n",
    "\n",
    "\n",
    "$$\\boldsymbol{x} = \\begin{pmatrix}\n",
    "1 \\\\\n",
    "2\\\\\n",
    "1 \\\\\n",
    "2\n",
    "\\end{pmatrix}\n",
    "+\n",
    "2\n",
    "\\begin{pmatrix}\n",
    "-2 \\\\\n",
    "-3\\\\\n",
    "1\\\\\n",
    "0\n",
    "\\end{pmatrix} \n",
    "+3\n",
    "\\begin{pmatrix}\n",
    "-1 \\\\\n",
    "-1\\\\\n",
    "0\\\\\n",
    "1\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "-6 \\\\\n",
    "-7\\\\\n",
    "3\\\\\n",
    "5\n",
    "\\end{pmatrix}$$\n",
    "\n",
    "So the presence of a null space leads to non-uniqueness of solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework - Minimal-norm solution to under-determined problem\n",
    "\n",
    "Let's begin this example from a problem already in RREF, the augmented form of the matrix with zero RHS being\n",
    "\n",
    "$$\n",
    "\\left(\n",
    "  \\begin{array}{cccc|c}\n",
    "1 & 0 & 4 & 0 & 0 \\\\\n",
    "0 & 1 & -2 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 1 & 0 \n",
    "  \\end{array}\n",
    "\\right)$$\n",
    "\n",
    "Show that the null space is given by any multiple of the following vector.\n",
    "\n",
    "$$\\boldsymbol{v}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "-4\\\\\n",
    "2\\\\\n",
    "1\\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Now consider the solution to the problem with RHS vector\n",
    "\n",
    "$$\\boldsymbol{b} = \\begin{pmatrix}\n",
    "1 \\\\\n",
    "-2\\\\\n",
    "3\n",
    "\\end{pmatrix} $$\n",
    "\n",
    "Use the minimum norm solution formula from the lecture to compute the solution.\n",
    "\n",
    "Establish that it is indeed the minimal-norm solution, i.e. that all other possible solutions you obtain by adding multiples of the null vector have a larger norm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "From the RREF the null space is described by the vector whose components satisfy\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "v_1 + 4v_3&=0\\\\\n",
    "v_2 - 2v_3&=0\\\\\n",
    "v_4&=0\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "As above, let's encode all solutions to this via an arbitrary value: $v_3:=\\alpha$, then $v_1=-4\\alpha$, $v_2=2\\alpha$, \n",
    "and $v_4=0$. So the null space is a multiple of the vector\n",
    "\n",
    "$$\\boldsymbol{v}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "-4\\\\\n",
    "2\\\\\n",
    "1\\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[ 0.23809524,  0.38095238,  0.        ],\n",
      "       [ 0.38095238,  0.80952381,  0.        ],\n",
      "       [ 0.19047619, -0.0952381 ,  0.        ],\n",
      "       [ 0.        ,  0.        ,  1.        ]])\n",
      "array([[ 1.0000000e+00, -4.4408921e-16,  0.0000000e+00],\n",
      "       [ 0.0000000e+00,  1.0000000e+00,  0.0000000e+00],\n",
      "       [ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00]])\n",
      "Min norm solution:  [-0.52380952 -1.23809524  0.38095238  3.        ]\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# come up with an example with indpt equations\n",
    "G = np.array([\n",
    "    [1, 0, 4, 0],\n",
    "    [0, 1, -2, 0],\n",
    "    [0, 0, 0, 1]])\n",
    "d = np.array([1, -2, 3])\n",
    "\n",
    "# construct the right inverse:\n",
    "\n",
    "G_ri = G.T @ sl.inv(G@G.T)\n",
    "\n",
    "pprint(G_ri)\n",
    "\n",
    "# print it to check its the identity (to round off error)\n",
    "pprint(G@G_ri)\n",
    "\n",
    "x_m = G_ri@d\n",
    "\n",
    "print('Min norm solution: ',x_m)\n",
    "\n",
    "# check that this is a solution:   Gx = d?\n",
    "print(np.allclose(d, G@x_m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should be the case that the addition of any multiple of the null space vector is also a solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x24e10f67208>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAEvCAYAAAA+brZ3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXST55328e/tDWPjBRsb8IYXmS0QCHsgYGwSkpJmaZq0WSdJaSgBknbani4z75npvHNm2r7TTjstWzaapUmaNktD07QJCTb7EkOAQCBYXrEN2Nh4Ae/2/f4h0XGJDZYt6ZYe/T7ncGxLsp7rkeyLW49k/ZTWGiGEsKog0wGEEMKTpOSEEJYmJSeEsDQpOSGEpUnJCSEsTUpOCGFpId7c2KhRo3R6ero3NymECAAHDhw4p7VO6Os8r5Zceno6hYWF3tykECIAKKXK+ztPHq4KISxNSk4IYWlSckIIS5OSE0JYmpScEMLSpOSEEJYmJSeEsDQpOSGEpUnJCSEszWdLrrG1k2e2l9DW2W06ihDCSw6U1/P+sTP09LjvHcu9+mddrjhW3ch/vHuc4WHBPDhvnOk4Qggv+MlfTlB1vpXFExIJC1JuuU6fXcldnxnPdWmxbNxWTGd3j+k4QggP219az0dl51mxKJOwEPdVk8+WnFKK1YttVJ5v5U+Hq03HEUJ42Lp8O/GRYXx1dppbr9dnSw4gb2IiE8dEsb6g2K2P0YUQvuWTyka2naxl+cIMhocFu/W6fbrkgoIUq3Jt2Gsu8P6nZ0zHEUJ4yPoCO1HhIR45/u7TJQdw69SxpMdHsDbfjsyIFcJ67DXN/PXYGR6Zn050eKjbr9/nSy44SLEyJ4ujVU1sLzpnOo4Qws3WFxQTHhLMowsyPHL9Pl9yAHfNSGFsTDjrttpNRxFCuNGp+hbePlTN/XPTiIsM88g2/KLkwkKCeGxhJvvL6tlfWm86jhDCTZ7aXkywUjy2MNNj2/CLkgO4b04a8ZFhrMuX1ZwQVlDT1MbvCyv58swUxsSEe2w7flNyw8OC+doNGWw7WcvRqkbTcYQQQ/TszlK6unt4PCfLo9vxm5IDeOj6cUSFh8hqTgg/d/5iB7/dW87t05JIi4/w6Lb8quSiw0N5+Pp0/nrsDPaaZtNxhBCD9PzuMlo6ulmVa/P4tvyq5AAeXZBOeEgw6wuKTUcRQgzChfYunt9dxtLJoxk/Osrj2/O7kosfMYz75qTx9qFqTtW3mI4jhHDRy3vLaWztZLUXVnHghyUH8NiiDIKU4+lnIYT/aOvs5pkdpSzMHsW01FivbNMvS25szHDunpnC7wsrqWlqMx1HCDFAvy88xbkL7axa7J1VHPhpyQF8Y1EWXd09PLuz1HQUIcQAdHb38NS2EmaOG8m8zDivbddvSy59VCS3TUvit3vLaWjpMB1HCHEVf/y4iqqGVtbk2lDKPe/6OxB+W3IAqxbbaOno5je7ykxHEUJcQXePZsO2YiaPjWbxhASvbtuvS27CmChumjya53eXcaG9y3QcIUQ//nr0DCW1F1nt5VUc+HnJAazOtdHY2snLe8tNRxFC9EFrzbp8O5kJkdwyZYzXt+/3JTc9NZYbbKN4dmepjC8UwgcVfFbLp6ebeDwni2A3TeByhd+XHDhWc7XN7fzhQKXpKEKIXrTWrM23kxw7nDuvSzaSwRIlNy8zjhlpsWwskPGFQviSfaX1HCg/zzdyMgkNNlM3lig5pRRr8mxUNbSy+ZCMLxTCV6zLtzNqxDC+MivVWAZLlBxA7oREJo2NZn2BnW4ZXyiEcUcqG9hRdI6vL8wgPNS9YwZdYZmSU0qxOjeL4tqLvHdMxhcKYdrarXaiPTRm0BWWKTmAL0wZS+aoSNbJ+EIhjDp5tpn3Pz3LIwsyGDEsxGgWS5VccJBi5eIsjlU3UXCy1nQcIQLW+nw7EWHBPDo/3XQUa5UcwJ3Tk0mKCWe9vEW6EEZU1LWw+XA1D8xNY6SHxgy6YsAlp5QKVkp9rJR6x/l1hlJqn1KqSCn1mlLK/N7gGF/4jZwsPio7z76SOtNxhAg4G7YVExIUxNc9OGbQFa6s5L4JHO/19U+BX2its4HzwHJ3BhuKr85OZdSIMNbKak4IrzrT2MYbByq5e1YKo6M9N2bQFQMqOaVUCnAr8KzzawXkAa87L/ICcKcnAg5GeGgwy2/IZEfROY5UNpiOI0TAeGZHCd1ae3zMoCsGupL7JfA94NKfE8QDDVrrS2/9UQmY+ZuNfjw4L41oGV8ohNfUX+zglX0V3DEtidQ4z44ZdMVVS04p9UWgRmt9oPfJfVy0z9dsKKVWKKUKlVKFtbXee8YzKjyUR+an896xsxSdlfGFQnjab3aV0trZzeOLfWcVBwNbyS0AbldKlQG/w/Ew9ZdArFLq0gtgUoA+/55Ka/201nqW1npWQoJ33yzv0QUZRITJ+EIhPK2prZPnd5dxyzVjyPbCmEFXXLXktNY/1FqnaK3TgXuBrVrrB4B84G7nxR4G3vZYykEaGRnGA3PT2Hy4moo6GV8ohKf8dm85zW1dXhsz6IqhvE7u+8C3lVJ2HMfonnNPJPd6bGEmwUqxUcYXCuERrR3dPLejlEXjE5iaEmM6zue4VHJa6wKt9Redn5doredorW1a63u01u2eiTg0idHh3DMrhdcLKznTKOMLhXC31z6qoO5iB2t8cBUHFvyLh76szMmiW2ue3VFiOooQltLR1cNT20uYnT6SORneGzPoioAoudS4CO6YlsTL+yqovyjjC4Vwlz9+XMXpxjafPBZ3SUCUHMCq3Czaurr5zS4ZRi2EO1waMzglOZqc8d595YQrAqbkbIlR3Dx5DM/vLqO5rdN0HCH83rufnKb03EVWL/b+mEFXBEzJgWPgTXNbFy/J+EIhhuTSmEFb4ghuvsb7YwZdEVAlNzUlhkXjE3huRymtHTK+UIjB+vB4DSfONLNqcRZBBsYMuiKgSg5gTa6NuosdvPZRhekoQvilS2MGU0YO57ZpSabjXFXAldycjDhmp4/kqe0ldHTJ+EIhXLWnuI5DpxpYmZNlbMygK3w/oQeszrVxurGNtz6WYdRCuGptvp3EqGHcPTPFdJQBCciSyxmfwJTkaDYUFMv4QiFccLDiPLuL63hsYabRMYOuCMiSU0qxerGNsroW/vzJadNxhPAb6/PtxEaEcv/cNNNRBiwgSw7g5mvGkJUQyfp8Oz2ymhPiqo6fbuKD4zV8bUEGkYbHDLoiYEsuKEixarGNE2ea2XqixnQcIXze+oJiRgwL4eHr001HcUnAlhzA7dOTSBk5nLUyjFqIKyo9d5E/H6nmwXnjiIkINR3HJQFdcqHBQazMyeLQqQb2FMv4QiH6s6HATmhwEMtvyDAdxWUBXXIAd89MITFqmIwvFKIfVQ2tvHmwintnp5IQNcx0HJcFfMmFhwbz2MJMdhfXcbDivOk4QvicZ7Y73odxhQ+NGXRFwJccwP1z04iNCGW9rOaE+Du1ze28ur+Cu2Ykkxw73HScQZGSAyKHhfDo/Aw+OF7D8dNNpuMI4TM27Sqls7uHlX66igMpub95ZH46I4aFyPhCIZwaWzp5aU85y6aOJTNhhOk4gyYl5xQTEcqD88bx5yPVlJ67aDqOEMa9sKeMC+2+OWbQFVJyvSy/IYPQ4CA2ympOBLiL7V1s2lXKkomJTBobbTrOkEjJ9ZIQNYx7Z6fy5seVVDe0mo4jhDGv7q+goaWT1Xn+vYoDKbnPWZGThdbw9HYZXygCU1tnN09vL2F+Vjwz0kaajjNkUnKXSY4dzpeuS+bV/RXUNvvkvGwhPOr1A5XUNLf7/bG4S6Tk+vD44iw6u3vYJOMLRYDp6u5h47ZipqfGMj8r3nQct5CS60NmwgiWTR3LS3vKaWyR8YUicGw+XE3l+VZW5/r2mEFXSMn1Y9ViGxfau3hhT5npKEJ4RU+PZn1BMRPHRLFkYqLpOG4jJdePyUnRLJmYyKZdpVxs7zIdRwiPe+/YGew1F1iVa/P5MYOukJK7gtV5NhpaOnlln4wvFNamtWZdgZ30+AhunTrWdBy3kpK7ghlpI5mfFc8zO0po65Rh1MK6tp2s5WhVEytzsgi20CoOpOSuanWujZrmdl4/IOMLhXWty7czNiacu2b4x5hBV0jJXcX8rHimp8aycVsxnd0yjFpYz76SOj4qO8+KRZmEhVivEqy3R26mlGJNro3K86386XC16ThCuN26gmLiI8O4d7b/jBl0hZTcACyZlMjEMVGsLyiW8YXCUo5UNrD9ZC3LF2YwPMw/hkW7SkpuAJRSrM61Ya+5wHvHzpiOI4TbrMu3ExUewkPzxpmO4jFScgO0bOpYMkZFyvhCYRlFZ5t579hZHpmfTlS4f40ZdIWU3AAFBykez8niWHUTBSdrTccRYsjWFxQzPDSYRxf435hBV0jJueDO65JJigln3VZZzQn/VlHXwubD1TwwN424yDDTcTxKSs4FYSFBfCMni8Ly8+wvrTcdR4hB27i9mGCleGxRpukoHicl56Kvzk5l1AgZRi3815nGNl4vrOSeWSmMjg43HcfjpORcFB4azNcXZrCj6ByHTzWYjiOEy57ZUUK31n49ZtAVUnKD8MDcNKLDQ1gnqznhZ+ovdvDKvgrumJZEalyE6TheISU3CFHhoTyyIIP3Pz3LZ2eaTccRYsA27SylraubVbmBsYoDKblBe3R+OhFhwawvkNWc8A9NbZ28sKeMW64Zgy0xynQcr5GSG6SRkWE8OG8cfzpcTZkMoxZ+4KU95TS3+f+waFdJyQ3B1xdmEBIcxFPbZRi18G2tHd1s2lnK4gkJTEmOMR3Hq6TkhiAxKpyvzkrl9QOVnG6UYdTCd726v4K6ix2sCbBVHEjJDdmKRZn0yDBq4cPauxzDoudmxDErPc50HK+7askppcKVUvuVUoeVUseUUv/mPD1DKbVPKVWklHpNKWXtvw3pR2pcBHdOdwyjPndBhlEL3/PmwSrONLUF3LG4SwaykmsH8rTW04DpwC1KqXnAT4FfaK2zgfPAcs/F9G2PL86ivauHTTtlGLXwLV3dPWwoKObalBgWZo8yHceIq5acdrjg/DLU+U8DecDrztNfAO70SEI/YEscwRemjHEMo26VYdTCd7xz5DQV9S2WGhbtqgEdk1NKBSulDgE1wBagGGjQWl8aSFoJJHsmon9YtdhGc3sXL+0pMx1FCMAxLHpdvp3xo0dw06TRpuMYM6CS01p3a62nAynAHGBSXxfr63uVUiuUUoVKqcLaWuu+D9uU5BhyJyTw3M5SWjpkGLUw7/1Pz1JUc4HVFhsW7SqXnl3VWjcABcA8IFYpFeI8KwXoc8qL1vpprfUsrfWshISEoWT1eWvybJyXYdTCB2jtWMWNs+CwaFcN5NnVBKVUrPPz4cCNwHEgH7jbebGHgbc9FdJfzBwXx7zMOJ7ZUUJ7lwyjFubsKDrHJ1WNPJ6TRUhwYL9SbCB7PxbIV0odAT4Ctmit3wG+D3xbKWUH4oHnPBfTf6zJzeZsUztvHKgyHUUEsLUWHhbtqpCrXUBrfQS4ro/TS3AcnxO9LLDFM805jPors1IC/n9R4X0fldWzv7Sef71tsiWHRbtKbgE3uzSMuqK+hT8dkWHUwvvWbrVbeli0q6TkPGDJRMcw6nX5MoxaeNeRyga2WXxYtKuk5DwgKEixSoZRCwPW5duJtviwaFdJyXnIrTKMWnhZoAyLdpWUnIfIMGrhbevy7USEBfOIxYdFu0pKzoNkGLXwlvK6iwEzLNpVUnIe1HsY9T4ZRi08aOO2YkKCg3hsofWHRbtKSs7DLg2jlvGFwlNON7by+oFKvjIrhcQAGBbtKik5D+s9jPqQDKMWHvD09hJ6NHxjUeCMGXSFlJwXPDhvHDHDQ2U1J9zu3IV2Xt1fwZ3TkwNmWLSrpOS8YMSwEB6Zn86WT89y4kyT6TjCQjbtLKW9qyeghkW7SkrOSx5dkE5kWDDr82V8oXCPxtZOXtpTzrIpY8lKGGE6js+SkvOS2AjHMOp3jsgwauEeL+4uo7m9S1ZxVyEl50XLF2YQGhzEhgJZzYmhudjexaZdpSyZmMg1SYE1LNpVUnJelBgVzr2zU3nz40qqGmQYtRi8V/dXcL6lk1UBOmbQFVJyXrYiJwut4eltspoTg9PW6RgWPT8rnpnjRpqO4/Ok5LwsOXY4d81I5ncfnaK2WYZRC9e9fqCSmuZ21sgqbkCk5Ax4fLGNzu4ent1ZYjqK8DOd3T1s3FbM9NRYrs+KNx3HL0jJGZAxKpJbr03it3vKaWjpMB1H+JHNh6qpPN/KmgAeFu0qKTlDVudmcbGjm9/sKjMdRfiJ7h7N+gI7E8dEsWRSouk4fkNKzpCJY6K5cdJont9dxoV2GUYtru6vR89QXHuR1bKKc4mUnEFr8mw0tnby273lpqMIH6e1Zm2+ncxRkSwL8GHRrpKSM2h6aiwLs0fx7I5S2jplGLXoX/5nNRw/3cTji7MIDpJVnCuk5AxbnWvj3IV2XvvolOkowkdprVm71U5y7HDuvC7ZdBy/IyVn2NyMOGaNG8lT24rp6OoxHUf4oD0ldRysaGBlTiahMqzcZXKLGaaUYk2ejerGNt76uNJ0HOGD1m61kxA1jHtmpZqO4pek5HxAzvgEpibHsKGgmK5uWc2J/3Ww4jy7i+tYsTCT8FAZFj0YUnI+QCnF6twsyupa+PMnp03HET5k3VY7sRGh3D83zXQUvyUl5yOWTh5DduII1uXb6emR8YUCjlU38uGJGpYvyCByWIjpOH5LSs5HBAUpVufaOHn2AluOnzUdR/iA9fnFRA0L4R/mp5uO4tek5HzIF68dy7j4CNblyzDqQGevucC7R0/z0PWOIUhi8KTkfEhIcBCP52RxpLKRHUXnTMcRBq0vsBMeEszyGzJMR/F7UnI+5q4ZKYyNCWftVhlfGKhO1bfw9qFq7puTRvyIYabj+D0pOR8TFhLEikWZ7C+rZ39pvek4woAN24oJVooVizJNR7EEKTkfdO/sNOIjw1grw6gDzpnGNl4vrOTuWSmMiQk3HccSpOR80PCwYL6+MJPtJ2s5fKrBdBzhRU9vL6Fbax7PkTGD7iIl56MenJdGdHgI62Q1FzDqLrTzyv5y7piWRGpchOk4liEl56OiwkN5ZEEG7396ls/ONJuOI7xg065S2rt6ZMygm0nJ+bBH56cTGRYsq7kA0NjayYu7y1k2ZSy2xBGm41iKlJwPGxkZxoPzxvHOkWpKz100HUd40Iu7y2hu72JVrhyLczcpOR+3fGEGocFBbCiQ1ZxVXWzv4rldpSyZmMg1STGm41iOlJyPS4wK597Zqbx5sIqqhlbTcYQHvLKvgoaWTlbnybE4T5CS8wMrcrJQCp7aVmw6inCzts5unt5RwgJbPDPSRpqOY0lScn4gOXY4X56Rwu8+OkVNU5vpOMKNfl94itrmdtbkZpuOYllScn5iZU4WXd09PLuz1HQU4SYdXT08ta2EmeNGMi8zznQcy5KS8xPpoyK5fVoSv91bTv3FDtNxhBv88WPHcdY1eTIs2pOk5PzI6lwbLR3dbJLVnN/r7tGsL7AzJTmaxeMTTMexNCk5P5I9OoovTBnDC3vKaGrrNB1HDME7R6opq2thTW62rOI8TErOz6zOtdHc1sVLe8pNRxGD1NOjWZdvZ/zoESydPNp0HMu7askppVKVUvlKqeNKqWNKqW86T49TSm1RShU5P8rz314wJTmG3AkJPLujhJaOLtNxxCBsOX6Wk2cvsDrXRlCQrOI8bSAruS7gO1rrScA8YLVSajLwA+BDrXU28KHza+EFa/KyOd/SySv7KkxHES7SWrN2q530+AhunTrWdJyAcNWS01qf1lofdH7eDBwHkoE7gBecF3sBuNNTIcXfmzluJPOz4nlqewltnd2m4wgXbDtZyydVjaxabCMkWI4WeYNLt7JSKh24DtgHjNZanwZHEQKJ7g4n+rcm10Ztczt/OFBpOooYoEuruKSYcO68Ltl0nIAx4JJTSo0A3gC+pbVucuH7ViilCpVShbW1tYPJKPpwfVY8M9Ji2VhQTGd3j+k4YgD2ldZTWH6elYuzCAuRVZy3DOiWVkqF4ii4l7XWbzpPPquUGus8fyxQ09f3aq2f1lrP0lrPSkiQ1wO5i1KKJ/KyqWpo5a2Pq0zHEQOwdqudhKhhfGVWqukoAWUgz64q4DnguNb6v3udtRl42Pn5w8Db7o8nrmTxhASuSYpmfb6d7h4ZRu3LDlacZ6f9HCsWZhIeGmw6TkAZyEpuAfAQkKeUOuT8twz4CXCTUqoIuMn5tfAix2rORlldC+8cqTYdR1zBuq12RkaEcv/cNNNRAk7I1S6gtd4J9PdiniXujSNctXTyGLITR7Au385t1ybJ66580NGqRj48UcN3bhpP5LCr/soJN5Ojn34uKEixJs/GybMXeP/Ts6bjiD6sL7ATFR7CwwvSTUcJSFJyFnDr1LGkx0ewNr8IreXYnC8pOtvMX46e4ZH56USHh5qOE5Ck5CwgJDiIVbk2jlY1UfCZvEzHl6zLtzM8NJhHF2SYjhKwpOQs4kvXJZMcO5xfbZXVnK8oO3eRzYereWjeOOIiw0zHCVhSchYRGhzEysVZfFzRwO7iOtNxBLChoJjQ4CCWL5RVnElSchZyz8wUEqOGsXarjC80raqhlTcOVnLfnDQSo8JNxwloUnIWEh4azIpFmewpqaOwrN50nIC2saAYpWDFokzTUQKelJzF3D83jfjIMH4tqzljzja18VrhKe6emUpS7HDTcQKelJzFRISFsHxhBttO1nKkssF0nID09PYSuns0j+dkmY4ikJKzpIfmjSNmeKis5gyou9DOK/squGNaEmnxEabjCKTkLCkqPJRH5qez5dOzHD894HfFEm7w7M5S2rq6WZVrMx1FOEnJWdTXFmQwYlgI6/JlNectDS0dvLi7jFunjsWWOMJ0HOEkJWdRMRGh/MP14/jzJ6ex11wwHScg/GZXGRc7ulmTJ6s4XyIlZ2HLb8ggPCSY9bKa87jmtk5+s6uUpZNHM3FMtOk4ohcpOQuLHzGMB+am8fbhasrrLpqOY2kv7imnqa1LVnE+SErO4lYsyiQ4SLGhoNh0FMtq6ejiuZ2l5IxP4NqUWNNxxGWk5CwuMTqce2en8sbBSqoaWk3HsaSX91ZQf7GDJ5fIKs4XSckFgJXOF6VulNWc27V1dvPU9hLmZ8Uzc1yc6TiiD1JyASApdjh3z0zhtcJTnG1qMx3HUn63v4JzF9p5Ii/bdBTRDym5APF4jo3uHs1T20pMR7GM9i7HKm7WuJHMy5RVnK+SkgsQafER3DE9iVf2l3PuQrvpOJbw+oFKTje28cSSbByTO4UvkpILIKtzbbR39fDMDlnNDVVndw8bCoqZlhLDouxRpuOIK5CSCyBZCSO47dokXtpTzvmLHabj+LW3Pq6i8nwrT8oqzudJyQWYNXk2Wjq62bSr1HQUv9XV3cP6fDvXJEWTNzHRdBxxFVJyAWb86Ci+MGUMz+8qo7G103Qcv/TOkdOU1bXwRJ5NVnF+QEouAK3Js9Hc3sXzu8pMR/E7PT2atfl2xo8ewdLJY0zHEQMgJReArkmK4cZJo9m0q5TmNlnNueIvR89gr7nAE3nZBAXJKs4fSMkFqCeX2Ghs7eTFPeWmo/iNnh7Nr7cWkZkQybKpY03HEQMkJRegrk2JZfGEBJ7bWUpLR5fpOH7h/U/PcuJMM0/k2QiWVZzfkJILYE/kZVN/sYOX91aYjuLztHas4sbFR3DbtUmm4wgXSMkFsJnjRnKDbRRPbS+htaPbdByftvVEDceqm1i92EZIsPza+BO5twLck0uyOXehnVf3y2quP1prfvVhESkjh/OlGcmm4wgXSckFuDkZcczLjGPjtmLaOmU115dtJ2s5XNnI6lwbobKK8ztyjwmeXJJNTXM7vy88ZTqKz9Fa8z8fFpEUE86XZ6SYjiMGQUpOcH1mPLPTR7KhoJj2LlnN9bbLXsfHFQ08nmsjLER+XfyR3GsCpRRPLsnmdGMbfyisNB3HZ1w6FjcmOpyvzJJVnL+SkhMA3GAbxYy0WDYUFNPR1WM6jk/YU1LH/rJ6VuZkMiwk2HQcMUhScgL439VcVUMrbxyU1RzArz4sIjFqGPfOSTMdRQyBlJz4m5zxCUxLjWVdvp3O7sBeze0rqWNvST0rc7IID5VVnD+TkhN/o5Tim0tsVJ5v5a2DVabjGPWrrUWMGjGM++fKKs7fScmJv5M7IZFrU2JYG8CrucKyenbZ61iZkymrOAuQkhN/RynFk3nZVNS38MePA3M19z8fFjFqRBgPzB1nOopwAyk58TlLJiUyJTmatfl2ugJsNVdYVs+OonN8Y1EWw8NkFWcFUnLicxzH5sZTXtfCHw9Vm47jVf/zYRHxkWE8ME+OxVmFlJzo042TErkmKZpfby0KmNXcgXLnKi4nk4iwENNxhJtIyYk+OVZz2ZTXtfB2gKzmfvmBYxX34Dw5FmclUnKiXzdNHs3ksYGxmru0iluxSFZxViMlJ/qllOJbN2ZTFgDH5i6t4h66XlZxViMlJ64oEFZzf3tGVY7FWdJVS04ptUkpVaOUOtrrtDil1BalVJHz40jPxhSmXFrNlde18JZFXzf3yw8cr4uTY3HWNJCV3PPALZed9gPgQ611NvCh82thUTdNHu18ptV6fwVRWFbPTrvjdXGyirOmq5ac1no7UH/ZyXcALzg/fwG40825hA9xrObGU1FvvdXcLz446fjrBnldnGUN9pjcaK31aQDnx0T3RRK+6MZJiUxNjuHXW4sss5rbV1Ln/BtVWcVZmcefeFBKrVBKFSqlCmtraz29OeEhSim+fdN4TtW38sYBa7zf3C8+OElC1DA5Fmdxgy25s0qpsQDOjzX9XVBr/bTWepbWelZCQsIgNyd8weIJCUxPjeXXW+1+/+7Bu4vPsbeknlWL5f3irG6wJbcZeNj5+cPA2+6JI3yZUop/vGk8VQ2tvObHk7201vxiy0lGRw/jPnnXX8sbyEtIXgX2ABOUUpVKqeXAT4CblFJFwE3Or0UAWJQ9ipnjRrJuq6uEVNEAAAwBSURBVN1v57TuKDrHR2XnWZ1rk1VcABjIs6v3aa3Haq1DtdYpWuvntNZ1WuslWuts58fLn30VFqWU4js3jedMUxuv7q8wHcdlWmt+vuUkybHD+ersVNNxhBfIXzwIl823jWJeZhzr8otp7fCv1dzWEzUcPtXAE3k2mcAVIKTkxKB8Z+kEzl1o58U9ZaajDFhPj+bn758kLS6CL8+UOaqBQkpODMrs9DhyxiewcVsxzW2dpuMMyHvHzvDp6Sa+dWM2ocHyox8o5J4Wg/adpeM539LJpp1lpqNcVXeP41hcVkIkd0xPNh1HeJGUnBi0a1Niufma0Ty7o4TzFztMx7mitw9VYa+5wHeWTiA4SJmOI7xISk4MyXeWTuBCRxcbtxebjtKvjq4efvHBSaYkR3PLNWNMxxFeJiUnhmT86CjunJ7MC7vLONvUZjpOn14rPMWp+la+u3QCQbKKCzhScmLIvnVjNl3dmrVb7aajfE5rRze//rCIOc4nSkTgkZITQzYuPpKvzk7l1f0VlNddNB3n7zy/u4ya5na+e/MElJJVXCCSkhNu8eSSbEKCFT9//6TpKH/T2NLJhgI7eRMTmZMRZzqOMERKTrjF6Ohwlt+QwebD1RytajQdB4D12+w0t3fxvVsmmI4iDJKSE27zjZwsYiNC+elfT5iOwunGVp7fVcaXrktm4pho03GEQVJywm2iw0NZk2tjR9E5dtnPGc3yyy1FaA3fvmm80RzCPCk54VYPzhtHcuxw/vPd4/T0aCMZTpxp4g8HTvHQ9eNIGRlhJIPwHVJywq3CQ4P53i0TOFbdZGzozX++e4Ko8FCeyLMZ2b7wLVJywu1uuzaJa1Ni+Nn7n3n9rZi2naxl+8lansizERsR5tVtC98kJSfcLihI8c/LJnG6sY1Nu0q9tt3uHs2P3z1OWlwED10vw2mEg5Sc8Ii5mfEsnTya9fl2arz0516/LzzFiTPNfP+WifKGmOJvpOSEx/zTskl0dmt+4oWXlDS2dPJf733GnPQ4lk2VP8IX/0tKTnhM+qhIli/M4M2DVRwoP+/Rbf3ig5M0tHTwo9uvkT/fEn9HSk541JpcG6Ojh/Gjzcc89pKSE2eaeGlvOQ/MHcfkJHnhr/h7UnLCoyKHhfBPyybxSVWjR2a1aq350eZjRIWHyAt/RZ+k5ITH3T4tibkZcfz43eNufxLiDwcq2VtSz/dunsjISHnJiPg8KTnhcUopfnzXVNq6evjRn4657Xprmtv4jz8fZ05GHPfKDFXRDyk54RWZCSP45pJs3v3kDO8dO+OW6/y3P31Ka2c3P75rqrzjr+iXlJzwmhWLMpk4Jop/efsoDS1DG3zz3rEz/PnIaZ7Ms5GVMMJNCYUVSckJrwkNDuJn90yj/mIH3/3DEbQe3LOtVQ2tfO/1I0xJjmbFoiw3pxRWIyUnvGpKcgw//MIkPjh+lk27ylz+/s7uHp545SDdPZq1980gLER+hMWVyU+I8LpHF6Rz0+TR/OQvxzl0qsGl7/3Z+59xsKKBH981lfRRkR5KKKxESk54nVKK/7r7WhKjwvn6Cx9hr2ke0Pe9sLuMp7aVcP/cNG6bluThlMIqpOSEEbERYby4fA6guO+ZfZTUXrji5X+7t5x/3XyMpZNH82+3X+OdkMISpOSEMVkJI3j1sbn09Gjue2Yv207Wfu4ybZ3d/PKDk/yfPx7lxkmJrL1/BqHB8mMrBk4N9hmuwZg1a5YuLCz02vaEf/jsTDPfeKmQsroWFo1P4L7ZqYQEB3G2qY31+XaqG9u4bVoSP7vnWnkLJdEnpdQBrfWsvs4L8XYYIS43YUwU7/9jDi/uKeNXHxaxvdeKbkpyNP/91enMy4w3F1D4NSk54RPCQoL4+sJMvjI7lYq6FgBCghXjE6PkrxnEkEjJCZ8SHR7KlOQY0zGEhcgRXCGEpUnJCSEsTUpOCGFpUnJCCEuTkhNCWJqUnBDC0qTkhBCWJiUnhLA0KTkhhKVJyQkhLM2r70KilKoFyl38tlHAOQ/E8fVtB/r2A3nfTW/fH/d9nNY6oa8zvFpyg6GUKuzvLVSsvO1A334g77vp7Vtt3+XhqhDC0qTkhBCW5g8l93SAbjvQtx/I+256+5bad58/JieEEEPhDys5IYQYNOMlp5S6Ryl1TCnVo5Saddl5P1RK2ZVSnymlbu7n+zOUUvuUUkVKqdeUUmFDyPKaUuqQ81+ZUupQP5crU0p94ryc2ybzKKV+pJSq6pVhWT+Xu8V5m9iVUj9w4/b/Syl1Qil1RCn1llIqtp/LuW3/r7YvSqlhzvvF7ryf04eyvcuuO1Upla+UOu78GfxmH5dZrJRq7HWf/Iu7tu+8/ivelsrhV879P6KUmuGm7U7otU+HlFJNSqlvXXYZt+67UmqTUqpGKXW012lxSqktzt/fLUqpkf1878POyxQppR52acNaa6P/gEnABKAAmNXr9MnAYWAYkAEUA8F9fP/vgXudn28EHndTrp8D/9LPeWXAKA/cFj8CvnuVywQ7b4tMIMx5G0120/aXAiHOz38K/NST+z+QfQFWARudn98LvObG23ssMMP5eRRwso/tLwbecfd9PdDbElgG/AVQwDxgnwcyBANncLzWzGP7DiwCZgBHe532/4AfOD//QV8/c0AcUOL8ONL5+ciBbtf4Sk5rfVxr/VkfZ90B/E5r3a61LgXswJzeF1BKKSAPeN150gvAnUPN5LzerwCvDvW6PGAOYNdal2itO4Df4bithkxr/b7Wusv55V4gxR3XewUD2Zc7cNyv4LiflzjvnyHTWp/WWh90ft4MHAeS3XHdbnQH8KJ22AvEKqXGunkbS4BirbWrL9R3idZ6O1B/2cm979/+fn9vBrZoreu11ueBLcAtA92u8ZK7gmTgVK+vK/n8D2A80NDrF7OvywzGQuCs1rqon/M18L5S6oBSaoUbttfbGufDkk39LN0Hcru4w9dwrCD64q79H8i+/O0yzvu5Ecf97lbOh8HXAfv6OPt6pdRhpdRflFLXuHnTV7stvXF/30v//6F7ct8BRmutT4PjPx0gsY/LDOk28Mq0LqXUB8CYPs76Z6312/19Wx+nXf5U8EAuM5gs93HlVdwCrXW1UioR2KKUOuH8X+qqrrR9YAPw7zj24d9xPGT+2uVX0cf3Dvgp8oHsv1Lqn4Eu4OV+rmbQ+395nD5OG/J97HIIpUYAbwDf0lo3XXb2QRwP4y44j5H+Ech24+avdlt6dP+dx7BvB37Yx9me3veBGtJt4JWS01rfOIhvqwRSe32dAlRfdplzOJbvIc7/5fu6jEtZlFIhwF3AzCtcR7XzY41S6i0cD7sG9Es+0NtCKfUM8E4fZw3kdhn09p0Hdb8ILNHOAyJ9XMeg9/8yA9mXS5epdN43MXz+Ic+gKaVCcRTcy1rrNy8/v3fpaa3fVUqtV0qN0lq75W87B3BbDun+HoAvAAe11mf7yObRfXc6q5Qaq7U+7XwYXtPHZSpxHB+8JAXHMfwB8eWHq5uBe53PrmXg+B9kf+8LOH8J84G7nSc9DPS3MhyoG4ETWuvKvs5USkUqpaIufY7jYP3Rvi7rqsuOtXypn+v9CMhWjmeVw3A81Njspu3fAnwfuF1r3dLPZdy5/wPZl8047ldw3M9b+ytfVzmP7T0HHNda/3c/lxlz6RigUmoOjt+ZOjdtfyC35WbgH5zPss4DGi89vHOTfh+1eHLfe+l9//b3+/sesFQpNdJ5CGep87SBcdczJ0N4xuVLOJq6HTgLvNfrvH/G8ezbZ8AXep3+LpDk/DwTR/nZgT8Aw4aY53lg5WWnJQHv9treYee/Yzge5rnrtngJ+AQ44rzzx16+fefXy3A8E1js5u3bcRz7OOT8t/Hy7bt7//vaF+D/4ihagHDn/Wp33s+ZbtzfG3A87DnSa5+XASsv/QwAa5z7eRjHkzHz3bj9Pm/Ly7avgHXO2+cTer0CwQ3bj8BRWjG9TvPYvuMo09NAp/N3fjmO46sfAkXOj3HOy84Cnu31vV9z/gzYgUdd2a78xYMQwtJ8+eGqEEIMmZScEMLSpOSEEJYmJSeEsDQpOSGEpUnJCSEsTUpOCGFpUnJCCEv7/1yavnX9UsQnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from numpy import linalg as nla\n",
    "\n",
    "# the null vector\n",
    "n = np.array([-4, 2, 1, 0])\n",
    "\n",
    "# add on a multiple and check it's still a solution\n",
    "mult = 1.\n",
    "\n",
    "x_p = x_m + mult*n\n",
    "\n",
    "print(np.allclose(d, G@x_p))\n",
    "\n",
    "\n",
    "# is x_m the minimum norm solution?\n",
    "\n",
    "# plot the norm of the vectors we get by adding on multiples of the numm vector\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "mult = np.linspace(-10,10,100)\n",
    "norms = []\n",
    "for m in mult:\n",
    "    print\n",
    "    norms.append(nla.norm(x_m + m*n))\n",
    "\n",
    "ax1.plot(mult,norms)\n",
    "\n",
    "# as we hoped for the norm is at a minimum when the multiplier is zero!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From which we see that if we add any multiple of the null vector to the min norm solution, the value of the norm does indeed increase."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.13"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}