{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Table of Contents](table_of_contents.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topic 6.  Projection Operators\n",
    "\n",
    "Authors: Sequoia Ploeg & Alexander Petrie     \n",
    "Email: sequoiap4@gmail.com, alexander.petrie@gmail.com"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Introduction\n",
    "\n",
    "A **projection** is a mapping of some space onto a subspace of the original space. A **projection operator** is a linear operator P that maps a space to the subspace. A **projection matrix** is an $n \\times n$ square matrix that maps from $\\mathbb{R}^n$ to a subspace $W$.\n",
    "\n",
    "In order for an operator to be a projection operator, it must be idempotent, i.e. $P^2 = P$. This means that whenever P is applied twice to any value, it gives the same result as if it were applied just once.\n",
    "\n",
    "There are two kinds of projections:\n",
    "* **Orthogonal:** A projection is an orthogonal projection if its range and nullspace are orthogonal. Another test is if $P^T = P$ for real matrices; $P^H = P$ for complex ones. Calculating least squares regressions require orthogonal projections.\n",
    "* **Oblique:** Non-orthogonal projections are still valid projections, i.e. $P^2 = P$, as will be shown later.\n",
    "\n",
    "One example of a projection is mapping a 3D object onto a 2D subspace. This is done in video games. The graphics are often generated in a 3D environment and they are then projected onto a 2D screen so we can see them. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Explanation of the theory\n",
    "Let $W$ be a finite dimentional vector space and $P$ be a projection on $W$. Suppose the subspaces $U$ avnd $V$ are the range and null space of $P$, respectively. \n",
    "\n",
    "### Properties of a projection operator\n",
    "* $P^2 = P$ ($P$ is *idempotent*)\n",
    "    * This means that if you use $P$ to try to operate on a value that has already been operated on by $P$, the value will not change.      \n",
    "    \n",
    "    \n",
    "* $\\forall \\, x \\in U \\, \\colon \\, Px = x$\n",
    "    * This means that if $P$ operates on some value in the range space of $P$, that value does not change.\n",
    "    \n",
    "    \n",
    "* There is a direct sum $W = U \\bigoplus V$. Every vector x in $W$ my be decomposed as $x = u + v$ where $u$ is in range of $P$ and $v$ is in the null space of $P$\n",
    "    * In other words: if $P\\,\\colon\\;W\\rightarrow W$, then $W = \\mathcal{R}(P) + \\mathcal{N}(P)$\n",
    "    \n",
    "    \n",
    "* The range and null space of a projection are complementary, meaning they are disjoint.\n",
    "\n",
    "\n",
    "* If $P$ is a projection, $(I - P)$ is also a projection.\n",
    "    * Proof: \n",
    "    $$\n",
    "    \\begin{split}\n",
    "        (I-P)^2 &= (I-P)(I-P) \\\\\n",
    "        &= I-P-P+P^2 \\\\\n",
    "        &= I - P - P + P \\\\\n",
    "        &= I-P \\\\\n",
    "    \\end{split}\n",
    "    $$\n",
    "    \n",
    "    \n",
    "* $P$ projects onto the range space $U$.  $(I-P)$ projects onto the null space $V$.\n",
    "\n",
    "### Orthogonal projections\n",
    "When $U$ (range space of $P$) and $V$ (null space of $P$) are orthogonal subspaces (i.e., $\\mathcal{R}(P)\\, \\bot\\,\\mathcal{N}(P)$), $P$ is an **orthogonal projection**. Orthogonal projections are very useful because if you can orthgonally project a value onto a subspace, that projection is the closest representation of the value on that subspace (see engineering example at the end). \n",
    "\n",
    "Two subspaces are orthogonal if for $\\forall u \\in U$, $\\forall v \\in V \\implies \\langle u, v \\rangle = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Numerical Example\n",
    "\n",
    "### Example 1\n",
    "\n",
    "This first example is an orthogonal projection of a vector in $\\mathbb{R}^2$ to the line y = 0. It orthogonally projects all points in $\\mathbb{R}^2$ to this line, meaning the result of the projection is the closest point on the line $y=0$ to the initial point. The \"physical\" intuition would be that, for any point in $\\mathbb{R}^2$, the closest point to the line $y=0$ is found by dropping a vertical line from that point to the x-axis. Therefore, in this example, the projection matrix simply \"selects\" out the x-component of the vector representing that point, discarding the y-component.\n",
    "\n",
    "The projection matrix is:\n",
    "$$\n",
    "   P=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & 0 \\\\\n",
    "   0 & 0 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$\n",
    "\n",
    "We can verify it is a projection by computing $P^2$ and showing it is equal to $P$:\n",
    "\n",
    "$$\n",
    "  P^2 = P \\cdot P = \\left[ {\\begin{array}{cc}\n",
    "   1 & 0 \\\\\n",
    "   0 & 0 \\\\\n",
    "  \\end{array} } \\right] \\cdot \\left[ {\\begin{array}{cc}\n",
    "   1 & 0 \\\\\n",
    "   0 & 0 \\\\\n",
    "  \\end{array} } \\right] = \\left[ {\\begin{array}{cc}\n",
    "   1 \\cdot 1 + 0 \\cdot 0 & 1 \\cdot 0 + 0 \\cdot 0 \\\\\n",
    "   0 \\cdot 1 + 0 \\cdot 0 & 0 \\cdot 0 + 0 \\cdot 0 \\\\\n",
    "  \\end{array} } \\right] = \\left[ {\\begin{array}{cc}\n",
    "   1 & 0 \\\\\n",
    "   0 & 0 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$\n",
    "\n",
    "The code below plots the projection of ten randomly generated $(x,y)$ coordinate pairs onto the line $y=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAEWCAYAAAD7MitWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucFOWd7/HPV3BUUC6KiigYjZdE3YXEicmYk3VcNCqr\nMTmHBJNdNsScBWU9G/aYK5qNL3XVZUN0TxI1JGtMvHOSYNBFDWBYd5dxA7rjBQVDUDNkRBgFFBG5\n/faP5xko2u6eC9Nd1V2/9+vVr+m6dNWvqp/6zVNPVfUjM8M55/Jgn7QDcM65avGE55zLDU94zrnc\n8ITnnMsNT3jOudzwhOecy42aTXiSpkv6Ua0stxvr/ZSkNkmbJH2g2uuvFEmTJP172nFUgySTdFx8\nf6ukb6Ydk9tTVROepJckvR0P6lcl/VjSgb1ZlpldZ2b/ey/jaZa0uq+X20vfBi4zswPN7L8KJ0q6\nRtIzkrZLuqrI9M9JelnSW5Lul3RwqRXFA/Ot+D38QdJ3JPXr282pPkkNkq6S9Nu4fS9Juk3Se6od\ni5ldYmbX7O1yJL0nfl/9+yKuvhL37Vkpx3CwpDnxu35Z0ue6+kwaNbwLzOxA4IPAh4ArC2dQULO1\nz146GlhWZvpK4KvAvxROkHQy8ANgInA4sBm4uYv1jY7fwxnABODiXsScNT8DPgF8DhgMjAaeAMb2\n5Uqylnxy7PvAVkKZ/3PglngslGZmVXsBLwFnJYb/EXgwvl8E/D3wH8DbwHHACGAu8DrhgP+rxGev\nAu5MDH8EWAxsAJ4CmhPTDgZ+DLQD64H7gYFxPTuBTfE1oshyP0FIRBtijO8v2J4vA08DG4H7gP1L\nbPs+hOT+MrAW+CnhoNwvrtuAt4DfdbEP7wSuKhh3HXB3Yvi9sSAcVGIZBhyXGJ4NfD8x/AXgeeBN\nYBUwJTGtGVgNXB634xXgC4nph8Tv7A3gN8A1wL8npp8OLIn7awlwemLaIuDa+D1uAh6Iy7srLm8J\n8J4S23RW/D5Hltl35crTfsBNsYy0x/f7FWzz14A1wB1x/Ffi9rcT/mHs2q/A7cC13dxnfwb8V9zG\ntuT3C/w+LrezjDbF8RfH72g98AhwdBwv4Ma4no2EsnlKL/bHVbFc/DSWg2VAY5x2B+G4eTvG9NWu\njpWC9X4fmFkw7gFgWg9yyUBCGT8hMe4O4Iayn6tEYisT5EvEhAeMjDvnmkRh/z1wMtAf2Bf4V0JN\nZX9gDLAOGJv4Qu6M748EXgPGERLL2XH40Dj9XwjJaGhc7hnJglgQY3K5JxCS0Nnxc1+NBaMhsT2/\niQXn4FgALymx7RfHzx4LHAj8gnjgFEtCZfZhsYT3S+BrBeM2AaeWWEbywHwf4QD824ID8L2Eg+cM\nQo3xg4l9th24Ou6TcXH60Dj9XsKBMhA4BfgDMeHFfbSeUBPtD3w2Dh+SKAMr47oHA88BLxCSWX/C\nwffjEtt0A/CvXey7cuXpauBx4DDgUELSvaZgm/+BkBgPAM4FXo3bOBC4m/IJr9w+awb+iFB2/zgu\n95Nx2nvicvsntuOTcT+9P+6XK4HFcdo5hFrtkPj9vR84ohf74ypgS4y1H3A98HixY7k7x0rBek8j\n/JPYJw4Pi/vj8Dj8ICFpFnt1VpA+ALxdsNwvAw9kLeFtioG/HHf2AYnCfnVi3pHADhK1lLjTby+S\nmL5GInnEcY8AnweOIPw3GloknmbKJ7xvArMT0/YhHMDNie35i8T0GcCtJbZ9ITA1MXwisK2zILN3\nCW8hBYk2GWeRZRihNvFWfH8PsTZTYv77gS8l9tnb7HkAriXUsPvFbXpfYtp17E54E4HfFCy7BZiU\nKANXJKbNBB5KDF8AtJaI8YfAvWW2oavy9DtgXGLaOcBLiW3eSqL2DtxGojZBOODLJbyi+6xErDcB\nN8b37+HdCe8h4IsF5XIzoVnkTwn/JD5CTCi93B9XAQsS004ikWB4d8Ire6wUWf/zwNnx/WXAvK7K\nfsHnPwasKRj3V8Cicp9Lo53sk2Y2xMyONrOpZvZ2Ylpb4v0I4HUzezMx7mVCba7Q0cCnJW3ofAH/\ng5DsRsblrO9FrCPiOgEws50xxmQMaxLvNxNqb10uK77vT2h/2FubgEEF4wYRTkVK+SAh1gnAhwm1\nFAAknSfpcUmvx305jvBfuNNrZrY9Mdy53YcStin5PSa3uXAfdE5P7s9XE+/fLjJcav++Rvi+S+mq\nPBX7fkYkhteZ2ZaC5ZXazqLxldhnSPqwpF9LWidpI3AJe+7vQkcD/5Qo668TanNHmtmjwPcIp42v\nSpolqbBsdMbf1fFVWLb3L9N+2Z1jJeknwF/E939BOB3tid6U+czdlmKJ9+3AwZIOSowbRfivUaiN\nUMMbkngNNLMb4rSDJQ3pYn3FtBMKFxAuphASaLEYurLHsgjbsp09D+jeWkZooAdA0rGEU68Xyn3I\ngtmEWtbfxc/uB/yccNX4cDMbAswjHFBdWUfYppGJcaMS7wv3Qef03uzPQguA0yQdVWJ6V+Wp2PfT\nnhguLCuvUHo7e+puQlvaSDMbDNzK7v1drIy2EdpVk+X9ADNbDGBm/8/MTiU0D51AaGss1JPjq5jC\nuHp6rNwJXChpNOG0+/7EZx+KdxAUez0UZ3sB6C/p+MQyR1P+wl/mEt4uZtZGaEe5XtL+kv4Y+CKh\nAbvQncAFks6R1C/O3yzpKDN7hXAKcLOkoZL2lfQn8XOvAodIGlwijNnAn0kaK2lfQqPzOzGunroH\n+FtJx8Rbca4D7iv4r19SjHt/wnfWP25j560kd8Xt/5ikgYS2ol8U/Pcu5wZgsqThQAMhWa4Dtks6\nD/h4dxZiZjsIbZNXSRog6SRCs0KnecAJ8Raa/pImEE6VHuxmnOXWvQCYD8yRdGpc/kGSLpF0cTfK\n0z3AlZIOlTSM8A/gzjKrnA1MknSSpAHAt/Yi/IMIta0tkk4jXGXutI7QJHNsYtytwDc6r0hKGizp\n0/H9h2KNcV9Ck8UWwqnrHnp4fBXzakFMPTpWzGw14SLUHcDPk2d6Znaehduzir3Oi/O8RShrV0sa\nKOmjwIV0UVPMbMKLPktow2gH5gDfMrP5hTPFL+9CYDqhgLQR/qt1bt9EQtvSckLbybT4ueWEgr4q\nnh6MKFjuCkJ1+7tAB6EN6QIz29qLbbmN8GU8BrxIKIj/pwef/yHhlO6zwBXx/cQY5zLCadBdcfsO\nAqZ2d8Fm9gyhAfsrMUn+DaEAryccfHN7EOdlhFO1NYR2rB8n1vMacD7hYHiN0LB9vpl19GD55Ywn\nJNX7CFconwUaCbU/KF+ergWWEq5qPgM8GccVZWYPEdraHiU0zj+6F3FPJRy4bxIS7ezEejYT716I\nZfQjZjaHcAHlXklvxO08L35kEKGsrCecYr5GqK0X063jq4TrCf8gNkj6ci+PlZ8QLtb09HS201TC\nBaS1hOP40ngslKTY2FdzJF0NHGVm9XD/mHO5E8+07iTcarSzGuvMeg2vqNg+cBKhpuScqzHxtPdL\nwI+qlewg5YQXH/tZK+nZxLiDJc1XeDxovqShRT76JOG+oS/E+T5fZB7nXAZJej/h1rQjCM0C1Vt3\nmqe0sUq7CfipmZ0Sx80gNODeIOnrhPvnvlbwuYMJ7S2NhKtFTxBusu3NrSfOuZxItYZnZo8R7iFK\nupDQmEn8+8kiHz0HmG9mnffXzSfc+e6ccyVl8SHow+OtJJjZK5IOKzLPkex50+dqStzgKGkyMBlg\n4MCBp77vfe/r43Cdc52eeOKJDjM7NO04SsliwuuOYjfBFj03N7NZwCyAxsZGW7p0aSXjci7XJHX1\nxEmqsniV9lVJRwDEv2uLzLOaPe9yP4o974p3zrl3yWLCm8vuu/M/T/glkEKPAB+PT04MJTwJ8EiV\n4nPO1ai0b0u5h/Ac54mSVkv6IuExp7Ml/ZbwUzM3xHkbFX963cxeJ/zO2pL4ujqOc865kmr2SYve\n8DY85ypL0hNm1ph2HKVk8ZTWOecqwhOecy43POE553LDE55zLjc84TnncsMTnnMuNzzhOedywxOe\ncy43POE553LDE55zLjc84TnncsMTnnMuNzzhOedywxOecy43POE553LDE55zLjc84TnncsMTnnMu\nNzKZ8CSdKKk18XpD0rSCeZolbUzM83dpxeucqw2Z7JfWzFYAYwAk9QP+AMwpMuu/mdn51YzNOVe7\nMlnDKzAW+J2ZZbqDX+dc9tVCwrsIuKfEtCZJT0l6SNLJ1QzKOVd7Mp3wJDUAnwD+f5HJTwJHm9lo\n4LvA/SWWMVnSUklL161bV7lgnXOZl+mEB5wHPGlmrxZOMLM3zGxTfD8P2FfSsCLzzTKzRjNrPPTQ\nQysfsXMus7Ke8D5LidNZScMlKb4/jbAtr1UxNudcjcnkVVoASQOAs4EpiXGXAJjZrcB44FJJ24G3\ngYvMzNKI1TlXGzKb8MxsM3BIwbhbE++/B3yv2nE552pX1k9pnXOuz3jCc87lhic851xueMJzmdPS\nAtdfH/4615cye9HC5VNLC4wdC1u3QkMDLFwITU1pR+XqhdfwXKYsWgRvT2hmx8Rmtm4Nw871FU94\nLlOam2GfWCobGsKwc33FE57LlKYmGD0ajjnGT2f7kreLBt6GlxEtLeH0rbnZD/JBg8Ir7/uhr3i7\n6G6e8DLAC6SrpM52UYCtdyxi0aL8li9PeCVUs8blBXJPnzn5M2mHUFeam2GflbBzp7eLesIroto1\nLi+Qe5r6oalph1BXmppg9BLYsAHuyvnZg1+0KKLat0Z4Q/2eNm/bzOZtm9MOo64MGgSjRnnZ8hpe\nEWnUuLyhfrdxd40DYNGkRekG4uqOJ7wi/BTA1RtvFw084ZVQ7RqXF0hXSd4uGnjCywgvkK6SOttE\nB+w7IOVI0uUJr4Rq17i8QLpK8nbRwBNeCdWucXmB3G3SmElph+DqVGYTnqSXgDeBHcB2M2ssmC7g\nn4BxwGZgkpk92Vfr9xpXejzhuUrJbMKLzjSzjhLTzgOOj68PA7fEv33Ca1zp6dgcvvJhA97VzbBz\ne6WWbzy+EPipBY8DQyQd0Zcr2LjyJP+FiRSMnz2e8bPHpx2Gq0NZruEZ8CtJBvzAzGYVTD8SaEsM\nr47jXknOJGkyMBlg1KhR3V75xpUn8fQ/fodndvgD/a72eTNBkOWE91Eza5d0GDBf0nIzeywxXUU+\n866OuGOinAXQ2NjY7Y66Ny4fw87PnQOoKg/0e4F0leTlK8hswjOz9vh3raQ5wGlAMuGtBkYmho8C\n2vtq/YPf1wrtBqaqPF7mBdJVkreLBplsw5M0UNJBne+BjwPPFsw2F/hLBR8BNprZK/SRL40/jVEn\nrK/aA/0dmzt2FUrn+pq3iwZZreEdDswJd57QH7jbzB6WdAmAmd0KzCPckrKScFvKF/oygEljJnF7\n6+1AddruOgujXxWGSxsvTTsEV6cymfDMbBUwusj4WxPvDfjrSsXQsbmDbTu3se8++1ZqFa6ECadM\nSDsEV6cyeUqbBeNnj2ftW2v9of4UtG1so21jW9czOtdDmazhZcWRBx3pD/WnYOKciYCf3ru+5wmv\njB07d7B522Z/vMzVPG8XDTzhlfHM2mcYd9e4qtQ0vEC6SvJ20cATXkZ4gXSV1NkmOnLwyC7mrG9+\n0aKESxsvZcRBI6q2Pm+od5U0cc7EXW2jeeY1vBImnDKBW5beUrX1eUP9bpc3XZ52CK5OecIroW1j\nG+9sf4f9+u+Xdii5c8GJF6QdgqtTfkpbwsQ5E1m/Zb0/45qCFR0rWNGxIu0wXB3yGl4Zww8c7gkv\nBVMenAL46b3re57wyti2Yxsdmzty/wsTrvZ5u2jgCa+MZeuWMX72+KrUNLxAukrydtHAE15GeIF0\nldTZJnrisBNTjiRdftGihMubLmfkoOrdpOkN9a6Spjw4ZVfbaJ55Da+EC068gJktM6u2Pm+o3+3K\nP7ky7RBcnfKEV8KKjhX+wwEpOevYs9IOwdUpP6UtYcqDU9i0dZM/1J+C1jWttK5pTTsMV4e8hlfG\nYQMP84f6UzDt4WmAn967vpe5hCdpJPBTYDiwE5hlZv9UME8z8EvgxTjqF2Z2dV/HsmX7Fto2tuX+\nFyZc7fN20SBzCQ/YDlxuZk/GnsuekDTfzJ4rmO/fzOz8SgayvGM5E+dMrEpNwwukqyRvFw0yl/Bi\nV4uvxPdvSnoeOBIoTHh1xQukq6TONtExw8ekHEm6Mn3RQtJ7gA8A/1lkcpOkpyQ9JOnkMsuYLGmp\npKXr1q3r9rqv/JMrOXrI0T0Nude8od5V0rSHp+1qG82zzNXwOkk6EPg5MM3M3iiY/CRwtJltkjQO\nuB84vthyzGwWMAugsbHRurv+s449i6H7D+1V7L3hDfW7XTf2urRDcHUqkwlP0r6EZHeXmf2icHoy\nAZrZPEk3SxpmZh19FUPrmlY2bd3EgQ0H9tUiXTedPvL0tENwdSpzp7SSBPwz8LyZfafEPMPjfEg6\njbAdr/VlHNMensaW7Vv8of4ULG5bzOK2xWmH4epQFmt4HwUmAs9I6mzUmg6MAjCzW4HxwKWStgNv\nAxeZWbdPV7tr2IBh/lB/CqYvnA746b3re5lLeGb274C6mOd7wPcqHcvmbZtZ0bEi978w4Wqft4sG\nmUt4WfLCay8w5cEpValpeIF0leTtooEnvIzwAukqqbNNNO/lLHMXLbLiurHXcezQY6u2Pm+od5U0\nfeH0XW2jeeY1vBJOH3k6g/YbVLX1eUP9bjede1PaIbg65QmvhMVti3njnTeqmvRckPfHn1zl+Clt\nCdMXTmfHzh3+UH8KFqxawIJVC9IOw9Uhr+GVMfSAof5QfwqufexaIB8/qNDSAosWQXMzNDWlHU39\n84RXxqatm2hd0+qnWK4iWlpg7FjYuhUaGmDhwsolPW8XDTzhlbHy9ZVMe3haVS4keIHMn0WL4O0J\nzQBsvWMRixZVLuH5P+3AE15GeIHMn+Zm2Gcl7NwZanjNzZVbV2ebaB6aCcrxixYl3HTuTRx38HFV\nW5831OdPUxOMHg3HHFPZ01kI7aKdbaN55jW8EsYMH1PVn4bKU0N9V35w/g/SDqFqBg0KL79gUR2e\n8EpYsGoB67esr+qPgLrAf6zBVYonvBKufexa9mEff6g/BQ+seAAgFz/N9ZmTP5N2CLniCa+MwfsP\nzv3D1mmY2TITyEfCm/qhqWmHkCue8MrYuGUji9sWe9JzFbN522YABuw7oKLryVO7aDme8Mp4ccOL\nTF84vSr34XmBzKdxd40DKv+jEd4uGnjCywgvkK6S8tQuWk5m78OTdK6kFZJWSvp6ken7SbovTv/P\n2IdteWvWwK9/vee4X/8aZszYPTxjBkyZwk+HfIETDjlh9zxTpuw5X2/NmLFnDDNmwHe+w3NfmbSr\nUL4rphLLaGmB668Pjyh1+Zli6+7Ouqop6/H1pWpta1zPzJaZu9pG63afdoeZZe4F9AN+BxwLNABP\nAScVzDMVuDW+vwi4r6vlnnrCCWbDhpk9+qiZWfibHO4cN2iQ2eDBdsaNo+2MG0fvGt5jvt4qXOfM\nmWaSfWnqe+2MH59RPKYiy9g6ZJid0/Co9etndk5DGO4yvsJld2dd1RTj+f28e+33G36fvfj6Uty2\nM24c3f3vPcvriYClloEcUuol66KzL0mXEfqHXV+ppFtknU3AVWZ2Thz+BoCZXZ+Y55E4T4uk/sAa\n4FArs0EHH/1+O3vqd+C552DECGhvh5NOgiFD9pxxwwbeeOEZVg3eCcCYjn5w8invnq+3NmzYM4ZR\no2jdvAoaGhjTvrN4TAXWvrCBP9jTbKWBD67ZyfojTuKwE7oRX+G6u7Guqsp6fH1pwwZa1z7do+89\n0+sBZl9y+hNm1liRhfeB7pzSDgeWSJodTzPL9ijWR44E2hLDq+O4ovOY2XZgI3BI4YIkTZa0VNLS\nbdu2wZAhtI7Yh9Z3Xg4HVbEvfsgQBh02kjFrYMwa4Mij+raADBkS1v1yjOGooxjTMIoxL79TOqYC\n+w8fwiGvjuLDa97hFY1g/+HdjK8725+mwn2Ttfj60pAhPf7eM72eWtCdaiCh28RzgHuBlcB1wHsr\nVe0EPg38KDE8EfhuwTzLgKMSw78DDim33FNPPdXs0UftjL/qb2d86+jSVfvO09oBA8wOOKDvTmcT\ny19+3BBb/s1LQwwzZ4a/3/xm90834mntv535ze6dziY+1+X2p6nzlKsn+6JWVWtbq7hPyfgpbU+S\n0GjgJmA5cAvwX8CMigQFTcAjieFvAN8omOcRoCm+7w90QDhFL/XqbMMr256RaMOzRx999/DeKmxT\niW14NnPmHtO7asPrVVtcldtzeizrbYx9qVrbWuV9mvWE1+UpraS/kfQEMAP4D+CPzOxS4FTgf/Wu\nXtmlJcDxko6R1EC4KDG3YJ65wOfj+/HAo3GHl/bWWzB79u4q/ZlnhuElSxJrXgIXXQRz5oTpZ54J\n998PEybsOV+vt2zJnjFs3w7f/nb4WyqmUss488zuf6bYurv7uWrp7XbVompta572aTd056LF1cA/\nm9nLRaa938yer0hg0jhCjbIfcJuZ/X2MZamZzZW0P3AH8AHgdeAiM1tVbpmNjY22dOlSmm9vBtLt\nISzNGDo2dwAwbMCwqq/b1TdJmb5o0eWNx2b2d2WmVSTZxWXPA+aVisXMthDa+nps0phJexVbrfNE\n5/Iql09a5D3h3d56O+D7weVPLhNeFk7p0uz+0ROey6tcJrzxs8cD6bbh+S8bO1d9mX2Wtt61rmml\ndU1r2mE4lyu5rOFlwbSHpwHp1jKdyxuv4TnncsNreDk078/ndT2Tc3Uolwnv0sZL0w4hVZX+OXHn\nsiqXCW/CKRPSDiFVNy+5GfAOZFz+5DLhtW0Mvzw1cvDI1GJIs/vH2ctmA57wXP7kMuFNnDMRSPcK\nqfeE5lz1+VXalCxuW8zitsVph+FcruSyhpcF0xdOB9KpZb7xRvgl9ZYWaGqq+uqdS43X8HKmpQWe\negpefBHGjo09njmXE17Dy5lFi0A/WQQ7YGu/MOy1PJcXuUx4lzddnnYIqWluhoYG2Lo1/G1uTjsi\n56onlwkvz72vNzXBwoWhZtfc7LU7ly+5THgrOlYAcOKwE1OL4aZzb0pt3U1NnuhcPmUq4Un6R+AC\nYCuh28UvmNmGIvO9BLwJ7AC29/Q39Kc8OAVI9z68McPHpLZu5/Iqa1dp5wOnmNkfAy8Qumcs5Uwz\nG5PlDkPKWbBqAQtWLUg7DOdyJVM1PDP7VWLwcUL3i3Xp2seuBfyXj52rpqzV8JIuBh4qMc2AX0l6\nQtLkcguRNFnSUklL161b1+dBOudqR9VreJIWAMOLTLrCzH4Z57kC2A7cVWIxHzWzdkmHAfMlLTez\nx4rNaGazgFkQ+qXd6w1wztWsqic8Myt7Difp88D5wFgr0Uu4mbXHv2slzQFOA4omvGLS7DHMOZee\nTLXhSToX+BpwhpltLjHPQGAfM3szvv84cHVP1uPtZs7lU6YSHvA9YD/CaSrA42Z2iaQRwI/MbBxw\nODAnTu8P3G1mD/dkJZ29haV5a8gPzv9Baut2Lq8ylfDM7LgS49uBcfH9KmD03qwnCz2GpXnTs3N5\nleWrtHXtgRUP8MCKB9IOw7lcyVQNL09mtswE8v1cr3PV5jU851xueMJzzuVGLk9p0+wxzDmXnlwm\nPO8xzLl8ymXC6+wtLM3Ed8en7kht3c7lVS4TXpo9hnVKsxNw5/LKL1qk5L5n7+O+Z+9LOwznciWX\nNbwsuGXpLQBMOGVCypE4lx9ew3PO5YYnPOdcbuTylDbNHsOcc+nJZcLzHsOcy6dcJrzO3sLS/CHQ\nn33mZ6mt27m8ymXCy0KPYcMGDEtt3c7llV+0SMntrbdze+vtaYfhXK54wkuJJzznqi9zCU/SVZL+\nIKk1vsaVmO9cSSskrZT09WrH6ZyrPVltw7vRzL5daqKkfsD3gbOB1cASSXPN7LlqBeicqz1ZTXhd\nOQ1YGTv0QdK9wIVAtxKe9xjmXD5l7pQ2ukzS05JukzS0yPQjgbbE8Oo47l0kTZa0VNLSdevWAaHH\nMO81zLn8SaWGJ2kBMLzIpCuAW4BrAIt/ZwIXFy6iyGet2LrMbBYwC6CxsdGAXb2FpdmBzrw/n5fa\nup3Lq1QSnpl16wY4ST8EHiwyaTWQ/EG5o4D27q4/Cz2GDdh3QGrrdi6vMndKK+mIxOCngGeLzLYE\nOF7SMZIagIuAudWIr6/cvORmbl5yc9phOJcrmUt4wAxJz0h6GjgT+FsASSMkzQMws+3AZcAjwPPA\nbDNbllbAvTF72WxmL5uddhjO5UrmrtKa2cQS49uBcYnheYA3hDnnui2LNTznnKuIzNXwqiELPYa9\n8QZs2AAtLdDUlHY0zuVDLhNe2j2GtbTAU0/Bzp0wdiwsXOhJz7lqyGXC6+wtLK0OdBYtAv1kEeyA\nrf3CsCc85yovlwkv7R7DmpuhoQG2bg1/m5tTCcO53MllwktbU1M4jV20KCQ7r905Vx2e8FLS1OSJ\nzrlq89tSnHO54QnPOZcbuTyl9R7DnMunXCY87zHMuXzK5Smtd6DjXD55wnPO5UYuE55zLp884Tnn\ncsMTnnMuNzzhOedyI5e3pXiPYc7lU6YSnqT7gM4OY4cAG8xsTJH5XgLeBHYA282ssSfr8R7DnMun\nTCU8M9v1e02SZgIby8x+ppl19GY9nb2FTf3Q1N583DlXozLZhidJwGeAeyqxfO8xzLl8ymTCAz4G\nvGpmvy0x3YBfSXpC0uQqxuWcq2FVP6WVtAAYXmTSFWb2y/j+s5Sv3X3UzNolHQbMl7TczB4rsb7J\nwGSAUaNG7UXkzrlaV/WEZ2ZnlZsuqT/wP4FTyyyjPf5dK2kOcBpQNOGZ2SxgFkBjY6OB9xjmXF5l\n8ZT2LGC5ma0uNlHSQEkHdb4HPg48292Fd/YY9uKLocewlpY+idk5VwMydZU2uoiC01lJI4Afmdk4\n4HBgTriuQX/gbjN7uLsL9x7DnMuvzCU8M5tUZFw7MC6+XwWM7u3yvccw5/Ircwmv0rzHMOfyK3cJ\nD7zHMOcdysU0AAAElUlEQVTyKosXLZxzriI84TnncsMTnnMuNzzhOedywxOecy43POE553LDE55z\nLjc84TnncsMTnnMuNzzhOedywxOecy43POE553LDE55zLjc84TnncsMTnnMuNzzhOedywxOecy43\nUkl4kj4taZmknZIaC6Z9Q9JKSSsknVPi88dI+k9Jv5V0n6SG6kTunKtladXwniX0PbtHX7KSTiL0\nWnYycC5ws6R+RT7/D8CNZnY8sB74YmXDdc7Vg1QSnpk9b2Yriky6ELjXzN4xsxeBlYROtndR6J/x\nT4GfxVE/AT5ZyXidc/Uha534HAk8nhheHcclHQJsMLPtZebZRdJkYHIcfEdStzvtriHDgI60g6iA\net0uqN9tOzHtAMqpWMKTtAAYXmTSFWb2y1IfKzLOejHP7glms4BZMaalZtZYat5a5dtVe+p12yQt\nTTuGciqW8MzsrF58bDUwMjF8FNBeME8HMERS/1jLKzaPc869S9ZuS5kLXCRpP0nHAMcDv0nOYGYG\n/BoYH0d9HihVY3TOuV3Sui3lU5JWA03Av0h6BMDMlgGzgeeAh4G/NrMd8TPzJI2Ii/ga8H8lrSS0\n6f1zN1c9qw83I0t8u2pPvW5bprdLocLknHP1L2untM45VzGe8JxzuVH3CW9vH2OrFZKukvQHSa3x\nNS7tmPaGpHPj97JS0tfTjqevSHpJ0jPxO8r0LRxdkXSbpLXJe1slHSxpfnzsc76koWnGWKjuEx57\n/xhbLbnRzMbE17y0g+mt+D18HzgPOAn4bPy+6sWZ8Tuq9fvwbiccO0lfBxbGxz4XxuHMqPuEtzeP\nsbnUnAasNLNVZrYVuJfwfbkMMbPHgNcLRl9IeNwTMvjYZ90nvDKOBNoSw2UfUasRl0l6Op5qZOpU\noofq8bvpZMCvJD0RH3usN4eb2SsA8e9hKcezh6w9S9srFXyMLVPKbSdwC3ANYRuuAWYCF1cvuj5V\nc99ND3zUzNolHQbMl7Q81pRcFdRFwqvgY2yZ0t3tlPRD4MEKh1NJNffddJeZtce/ayXNIZy+11PC\ne1XSEWb2iqQjgLVpB5SU51PaLh9jqyWxcHX6FOFiTa1aAhwff+i1gXBxaW7KMe01SQMlHdT5Hvg4\ntf09FTOX8LgnZPCxz7qo4ZUj6VPAd4FDCY+xtZrZOWa2TFLnY2zbSTzGVqNmSBpDOPV7CZiSbji9\nZ2bbJV0GPAL0A26Ljx3WusOBOeEnHekP3G1mD6cbUu9JugdoBobFR0W/BdwAzJb0ReD3wKfTi/Dd\n/NEy51xu5PmU1jmXM57wnHO54QnPOZcbnvCcc7nhCc85lxue8JxzueEJzzmXG57wXOZI+lD8EYT9\n49MJyySdknZcrvb5jccukyRdC+wPHACsNrPrUw7J1QFPeC6T4jO0S4AtwOk1/tifywg/pXVZdTBw\nIHAQoabn3F7zGp7LJElzCb90fAxwhJldlnJIrg7U/a+luNoj6S+B7WZ2d+zfYrGkPzWzR9OOzdU2\nr+E553LD2/Ccc7nhCc85lxue8JxzueEJzzmXG57wnHO54QnPOZcbnvCcc7nx3yKrUGJ5Qn8xAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x178e9658f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = ((np.random.rand(10)*2)-1)*10\n",
    "y = ((np.random.rand(10)*2)-1)*10\n",
    "pair = zip(x, y)\n",
    "\n",
    "# The projection matrix\n",
    "P = np.array([[1, 0],\n",
    "              [0, 0]])\n",
    "\n",
    "for item in pair:\n",
    "    x, y = item\n",
    "\n",
    "    # Form the vector z in R2 from the point (x, y)\n",
    "    z = np.array([[x],\n",
    "                  [y]])\n",
    "\n",
    "    # Project z onto the line y = 0\n",
    "    Pz = P @ z\n",
    "\n",
    "    # Plot the result\n",
    "    plt.plot(x, y, 'b.')\n",
    "    plt.plot(Pz[0], Pz[1], 'rx')\n",
    "    plt.plot([x,Pz[0]],[y,Pz[1]],color='green', linestyle='dashed')\n",
    "\n",
    "plt.plot([-10,10],[0,0])\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.axis('square')\n",
    "plt.ylim(-10, 10)\n",
    "plt.xlim(-10, 10)\n",
    "plt.title(\"Projection of 10 Random Coordinates onto y=0\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2\n",
    "\n",
    "This example is similar to the first in that it maps all points in $\\mathbb{R}^2$ to a line. However, now the line is y = x. \n",
    "The projection matrix in this case is: \n",
    "\n",
    "$$\n",
    "   P=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   0.5 & 0.5 \\\\\n",
    "   0.5 & 0.5 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$    \n",
    "You can change the x and y points in the code and it will show how they project onto the line y = x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAEWCAYAAAD7MitWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXJwkBWYMa2SVWxH2pshhrayIgmwtuCC7V\nakv9KWpt+9VqXfhq60K1tlWUoqLiBgguVBYFAflaoxAsyKKFqFF2AsguBMjn98e9AzfDTDJZZu6d\nmc/z8ZhHZubeuXPuzJ137j3n3HNFVTHGmHSQ4XcBjDEmUSzwjDFpwwLPGJM2LPCMMWnDAs8YkzYs\n8IwxaSNpA09E7haR55JluTG878UiskJEtovIjxP9/vEiIteJyEd+lyMRRERFpJN7f6SI3Ot3mUxl\nCQ08ESkVkR/cH/U6EXlBRJrWZlmq+pCq/rKO5SkQkZX1vdxaegwYqqpNVfU/4RNF5EERWSQie0Vk\nWITpV4rItyKyQ0TeFpFDo72R+8Pc4X4Pq0TkryKSWb+rk3giki0iw0Rkubt+pSIyWkTyEl0WVb1R\nVR+s63JEJM/9vrLqo1z1xf1se/pdjpryYw/vAlVtCpwOdAXuCZ9BHEm791lLHYElVUwvAe4AJodP\nEJETgX8C1wCtgJ3A09W836nu93AOcAVwfS3KHDQTgAuBK4EWwKnAfKBHfb5J0MLH1ICqJuwGlAI9\nPY//Arzr3p8N/Bn4N/AD0AloC0wCNuH84H/lee0w4BXP4zOBj4HNwEKgwDPtUOAFYDXwPfA20MR9\nnwpgu3trG2G5F+IE0Wa3jMeHrc/vgc+BLcA4oFGUdc/ACfdvgfXAGJwfZUP3vRXYAXxVzWf4CjAs\n7LmHgNc8j48GyoFmUZahQCfP4/HACM/jXwBfANuAr4Ffe6YVACuB37nrsQb4hWf6Ye53thWYCzwI\nfOSZfhYwz/285gFneabNBv7kfo/bgX+5y3vVXd48IC/KOvV0v88OVXx2VW1PDYG/udvIavd+w7B1\nvhNYC7zsPv8/7vqvxvmHsf9zBV4E/hTjZ9Yf+I+7jiu83y/wnbvc0Daa7z5/vfsdfQ+8B3R0nxfg\nCfd9tuBsmyfV4vMY5m4XY9ztYAnQxZ32Ms7v5ge3THdU91sJe98RwONhz/0L+E0NsuRot9yne9Zl\nA57ffcTXxSvcohSyFDfwgA7uh/OgZ2P/DjgRyAIaAB/i7Kk0Ak4DyoAeni/kFfd+O2Aj0A8nWHq5\nj3Pd6ZNxwqilu9xzvBtiWBm9y+2ME0K93Nfd4W4Y2Z71met+2Ie6G+CNUdb9eve1PwKaAm/i/nAi\nhVAVn2GkwHsHuDPsue3AGVGW4f1hHofzA7w97Ad4NM6P5xycPcbQhlUA7AUecD+Tfu70lu70sTg/\nlCbAScAq3MBzP6PvcfZEs4DB7uPDPNtAifveLYClwDKcMMvC+fG9EGWdHgE+rOazq2p7egD4BDgC\nyMUJ3QfD1vlRnGA8BOgDrHPXsQnwGlUHXlWfWQFwMs62e4q73AHutDx3uVme9Rjgfk7Hu5/LPcDH\n7rTeOHu1Oe73dzzQphafxzBgl1vWTOBh4JNIv+VYfith79sN559Ehvv4cPfzaOU+fhcnNCPd3vUs\n51c4v7nGOKH/WLW/Hx8Cb7tb8G/dD/sQz8b+gGfeDsA+PHsp7of+YoRguhNPeLjPvQdcC7TB+W/U\nMkJ5Cqg68O4FxnumZeD8gAs863O1Z/pwYGSUdf8AuMnz+FhgT2hDpm6B9wFhQestZ4RlKM7exA73\n/uu4ezNR5n8buM3zmf1A5R/gepw97Ex3nY7zTHuIA4F3DTA3bNlFwHWebeCPnmmPA1M9jy8AFkQp\n47PA2CrWobrt6Sugn2dab6DUs87lePbegdHAI57Hnak68CJ+ZlHK+jfgCfd+HgcH3lTghrDtcidO\ntci5OP8kzsQNlFp+HsOAGZ5pJwA/hP2WvYFX5W8lwvt/AfRy7w8FplS37UdZziRgEc6ebNRtOHTz\no55sgKrmqGpHVb1JVX/wTFvhud8W2KSq2zzPfYuzNxeuI3C5iGwO3YCzccKug7uc72tR1rbuewKg\nqhVuGb1lWOu5vxNn763aZbn3s3Dq3OpqO9A87LnmOIci0ZyOU9YrgO44eykAiEhfEflERDa5n2U/\nnP/CIRtVda/ncWi9c3HWyfs9etc5/DMITfd+nus893+I8Dja57sR5/uOprrtKdL309bzuExVd4Ut\nL9p6RixflM8MEekuIrNEpExEtgA3UvnzDtcR+LtnW9+EszfXTlVnAk/hHDauE5FRIhK+bYTKX93v\nK3zbblRF/WUsvxWvl4Cr3ftX4xwm18azOHvZT6rq7upmDlrDgHrurwYOFZFmnueOxPmvEW4Fzh5e\njufWRFUfcacdKiI51bxfJKtxNi7AaUzBCdBIZahOpWXhrMteKv+ga2sJTgU9ACLyI5xDr2VVvUgd\n43H2su5zX9sQmIjTatxKVXOAKTg/qOqU4axTB89zR3ruh38Goem1+TzDzQC6iUj7KNOr254ifT+r\nPY/Dt5U1RF/PmnoNZ0+lg6q2AEZy4POOtI2uwKlX9W7vh6jqxwCq+g9VPQOneqgzTl1juJr8viIJ\nL1dNfyuvABeJyKk4h91ve1471e1BEOk21TNfU5y94eeBYVX1TAgJWuDtp6orcOpRHhaRRiJyCnAD\nTgV2uFeAC0Skt4hkuvMXiEh7VV2DcwjwtIi0FJEGIvIz93XrgMNEpEWUYowH+otIDxFpgFPpvNst\nV029DtwuIke5X9RDwLiw//pRueVuhPOdZbnrGOpK8qq7/j8VkSY4dUVvhv33rsojwBARaQ1k44Rl\nGbBXRPoC58WyEFXdh1M3OUxEGovICTjVCiFTgM5uF5osEbkC51Dp3RjLWdV7zwCmA2+JyBnu8puJ\nyI0icn0M29PrwD0ikisih+P8A3ilirccD1wnIieISGPg/joUvxnO3tYuEemG08ocUoZTJfMjz3Mj\ngbvc1nlEpIWIXO7e7+ruMTbAqbLYhXPoWkkNf1+RrAsrU41+K6q6EqcR6mVgovdIT1X7qtM9K9Kt\nr2cxfwfmq9ONbLL7uVQpsIHnGoxTh7EaeAu4X1Wnh8/kfnkXAXfjbCArcP6rhdbvGpy6pS9x6k5+\n477uS5wN/Wv38KBt2HL/i7O7/SROC9AFON1qymuxLqNxvtw5wDc4G+ItNXj9sziHdIOBP7r3r3HL\nuQTnMOhVd/2aATfFumBVXYRTgf0/bkjeirMBf4/z45tUg3IOxTlUW4tTj/WC5302Aufj/Bg24lRs\nn6+qG2qw/KpchhOq43BaKBcDXXD2/qDq7elPQDFOXdAi4DP3uYhUdSrO3sVMnMr5mXUo903AAyKy\nDSdox3veZydu7wV3Gz1TVd/CaUAZKyJb3fUMBUFznG3le5xDzI04e+uRxPT7iuJhnH8Qm0Xk97X8\nrbyE01hT48NZEbkIp+HoRvep3wKni8hVVb7OrfhLOiLyANBeVVOh/5gxacc90noFp6tRRSLeM+h7\neBG59QMn4OwpGWOSjHvYexvwXKLCDnwOPPe0n/Uistjz3KEiMt09PWi6iLSM8NLPcPoN/cKd79oI\n8xhjAkhEjsfpmtYGp1ogce/t5yGtu0u7HRijqie5zw3HqcB9RET+gNN/7s6w1x2KU9/SBae1aD5O\nJ9vadD0xxqQJX/fwVHUOTh8ir4twKjNx/w6I8NLewHRVDfWvm45TgWmMMVEF8SToVm5XElR1jYgc\nEWGedlTu9LmSKB0cRWQIMASgSZMmZxx33HH1XFxjzO69FXxdtp2dq5dvUNVcv8sTTRADLxaROsFG\nPDZX1VHAKIAuXbpocXFxPMtlTNopWb+NQaM+5Xhg/r29qjvjxFdBbKVdJyJtANy/6yPMs5LKvdzb\nU7lXvDEmAUJhBzB2SHefS1O9IAbeJA70zr8WZySQcO8B57lnTrTEORPgvQSVzxjDwWHX6Yhm1bzC\nf353S3kd5zzOY0VkpYjcgHOaUy8RWY4z1Mwj7rxdxB16XVU34YyzNs+9PeA+Z4xJgGQMO0jiMy1q\nw+rwjKm7qsJOROarahe/yladIB7SGmMCKln37EIs8IwxMUn2sAMLPGNSRlERPPyw87e+pULYQfL2\nwzPGeBQVQY8esFt3kP2nRsyckUl+fv0sO1XCDmwPz5iUMHs27N5bTsWVvdnd+3pmzjpozM9aSaWw\nAws8Y1JCQQE0zMpGvu6NnjqGfx9+Pfsq6hZ6qRZ2YIe0xqSE/Hz44AOYPftevukAz351H9dPgtEX\njiYzI7P6BYRJxbADCzxjkkpRkXP4WlDAQXV0+fmh5+6lw4dw3+z7yG2cy2PnRRvhPbJUDTuwwDMm\naexvmGiwloYPtuaDDw4OvZB7z7mX5g2b0++YfjV6j1QOO7A6PGOSxuzZsOvEUVTcdCy7D/+U2bOr\nnv+2M2/jmMOOQVV59fNXq63TS/WwAws8Y5JGQQFkf9cXduZScfV55J72aUyve/+r97n6rau5flL0\nhox0CDuwwDMmaeTnw6y3O3BHq1m0y8nldwvP49OV1Yde7069eaDgAcYsHBMx9NIl7MACz5ikkp8P\nj97dgaJfzyK3cS69X+nNhp3VX9b33nPuPRB6/+y3P/RK1m9j0FNzYMf2lA87sEYLY5JShxYdmHXt\nLD789kMOb3x4TK+595x74ZtvePDrF7jl3WfJOfMqJ+w2b2Fsr1YpH3Zgw0MZkxJmfTOLJtlN6Nau\nW7XzLps8hozbHmPQZf8LO39wwu6CHvVSjqAPD2V7eMYkuX0V+7ht2m18u+Vbpl8zvdrQy+h6MYMu\nawbbdzC25Xd0uuDKBJXUf1aHZ0ySy8zIZPKVk8ltnEuvl3sxd9XcqPPur7Pb+YMTdiOfgFmzElha\nf1ngGZMCQnV6VYXe/rAL1dk9eDeMHw8DB6ZN6AUy8ETkWBFZ4LltFZHfhM1TICJbPPPc51d5jQkC\nb+iNeeE3lULsQGvszsp1doWFTujNm+dTqRMr8I0WIpIJrAK6q+q3nucLgN+r6vmxLssaLUw6WL9j\nPYd9+jmZVwxGx43jq5O6Vt6zq6cGikis0aLuegBfecPOGBPdEU2OgHN7smLMU5z/1vXw/uM02Lk7\n7mGXDAJ5SBtmEPB6lGn5IrJQRKaKyImJLJQxgTN8eKXD2NLOp7O98UNk7PyBsS1K0z7sIOCBJyLZ\nwIXAGxEmfwZ0VNVTgSeBt6MsY4iIFItIcVlZWfwKa4zfunbd3wBRsn4btz+/lJY7Knhje1HatcZG\nE+jAA/oCn6nquvAJqrpVVbe796cADUTkoC7nqjpKVbuoapfc3Nz4l9gYv7gNECW/uo2Bf5kG27Yx\n9sjNdBo9Iu1aY6MJeuANJsrhrIi0FhFx73fDWZeNCSybMYFTcmIXBl32v+zat5t2Ox8ka4jbppdm\nrbHRBDbwRKQx0At40/PcjSJyo/vwMmCxiCwE/gEM0qA3ORsTR95OxQ83mc+cI5ZTODKf0s2lzgyF\nhXDHHb6W0W+B75ZSn6xbiklVB3UqvqAH8yf9k55F/48mjVtxZWYRFxfm1dulG6MJereUwO7hGWNi\ns388uz17KnU9OePCX/PPVs+wcccWHps0mR494nOR7mSSDP3wjDFRVBq88/YeBw3x9NUPv6b8mYvQ\nra0pz3SGiY/3Xl6Q2R6eMUkqlpGKCwqg4Z7WZGZCdrbzOJ3ZHl6aqupyfyb4Yh2W/cD1au27Bgu8\ntFRUBOf2UMpbfUTDB39a5eX+TPDU9BoUB65Xa+yQNg3Nng27j3+Riut+xq5TRlR7uT8THOl0wZ14\nsMBLQwUF0PC/V8GXA9C+Q9n4oxF+F8nEwMKu7izw0lB+Psycns0DJ4/jZ0cM4PEvhzJiroVekFnY\n1Q/reJzmyveVc8WEK5i8bDLLbllGXk6e30UyYZIp7ILe8dgaLdJcdmY24y4bR/HqYgu7AEqmsEsG\ndkhryM7M5qwOZwEwfsl4O7wNCAu7+md7eGY/VWX8kvFM/GIiADd3u9nnEqUvC7v4sD08s5+I8Nql\nr3HRsRcxdKo1ZPjFwi5+LPBMJdmZ2Yy/fLyFnk8s7OLLAs8cxBt6K7au8Ls4acPCLv6sDs9ElJ2Z\nzYSBE8iUTAA279pMTqMcn0uVuizsEsP28ExUWRlZiAhff/81xz11HE/NfcrvIqUkC7vEscAz1Wrf\nvD1ntj+TW6beYqFXzyzsEssCz1TLW6dnoVd/LOwSL7CBJyKlIrJIRBaIyEHng4njHyJSIiKfi8jp\nfpQzXYSH3ptfvFn9i1JAURE8/HD9D41uYeePoDdaFKrqhijT+gLHuLfuwDPuXxMnodB79KNH6dOp\nj9/FibuiIujRA3Y3XEXDB9vV27iBFnb+CeweXgwuAsao4xMgR0Ta+F2oVJedmc2959xL4waN2bJr\nC+OXjPe7SHEzezbsOuF5Km46lt2t59TLuIEWdv4KcuAp8L6IzBeRIRGmtwO8ncRWus9VIiJDRKRY\nRIrLysriVNT09JeP/8IVE65I2Tq9ggJo+G1/2HIkFYP70uKUOXVanoWd/4IceD9R1dNxDl1vFpGf\nhU2XCK85aKwrVR2lql1UtUtubm48ypm27jvnvpRuyMjPh5mTWnNX25nktezI/yzsy5xvaxd6FnbB\nENjAU9XV7t/1wFtAt7BZVgIdPI/bA6sTUzoD6dF6m58PD93dmqIbZ9KxRUf6v9afddvX1WgZFnbB\nEcjAE5EmItIsdB84D1gcNtsk4Odua+2ZwBZVXZPgoqY9b+gN//dwtu3e5neR4qJ109bMvHYmo84f\nRaumrWJ+nYVdsAS1lbYV8JaIgFPG11R1mojcCKCqI4EpQD+gBNgJ/MKnsqa9UOit37GeZg2boaq4\n311Kad20NYNPHgzArG9mkZmRyc86hte0HGBhFzw2xLupVxVawdApQzkh9wSGdhvqd3HiokIr6P5c\nd5aWLWXqVVMjhl66hl3Qh3gP5CGtSV77KvaxetvqlK3TA8iQDP41+F90bNGRvq8e3JCRrmGXDCzw\nTL1qkNkg5Rsy4ECdXscWHen7Yk/mvPU3wBN25bsZm7nYwi5gLPBMvQtvvU3VQUT3h17jtrz+6l2U\n/OuDA2H32l10OvNUv4towgS10cIkuVDoXfPWNSl9NbTWTVvz0dDPKDu2mEHT10HjQxg74X46Pft3\nKCz0u3gmjAWeiZvQJSBDlm1cRufDOvtYovjYtK0BV32e4YTdc7fS6aZfWNgFlB3SmoSY+c1Mjh9x\nfMrV6VWqs5twvxN2zzwDs2b5XTQTgQWeSYizjzybCzpfkFINGZXC7rW7nMPYBx6A8eNh4EALvQCy\nwDMJkWqnoVXqenLI8sp1doWFTujNm+djCU0k1vHYJFT5vnIGvjGQd/77Dp/c8And2yffEIbWzy66\noHc8tkYLk1ChPb2JSyfSrV34eBDBZ2GX3OyQ1iRcdmY2g08ejIiwcO1Cnp3/rN9FiomFXfKzPTzj\nq79/+ndeWPAC5fvKubnbzX4XJyoLu9RggWd8NfL8kWz6YRNDpzoDDQQx9CzsUocd0hpfeVtvh04d\nGozT0IYP39+lxM6NTS0WeMZ33tB7d/m7VGiFvwXq2hUGDrRzY1OQHdKaQAiFXoVWkCEZ7Nm3hwaZ\nDfwpTGEhJaPHcsWMtXBIY8bZubEpw/bwTGBkZ2bTKKsRW3Zt4azRZ/l2eFuyfhuDFipySGPGPXcr\nna4cwK6f1sMFaY3vLPBM4BzS4BDaN2/vS51epHNjH/n0cc7++6l8/8P3CS2LqX+BCzwR6SAis0Tk\nCxFZIiK3RZinQES2iMgC93afH2U18REaZWXAcQMSGnrRzo09+Vf3sGjLMno9faaFXpILXOABe4Hf\nqerxwJk416Q9IcJ8/6eqp7m3BxJbRBNv4aH3wn9eiOv7VXVubP/L7uLNUx5i0fav6PVyLwu9JBa4\nwFPVNar6mXt/G/AF0M7fUhk/hELv5q43c+5R59bvwqvqenLXbw5qoOh/2V28OfgdFq1fRO9XerO3\nYm/9lsckRKBbaUUkD/gx8GmEyfkishDn4tu/V9UlUZYxBBgCcOSRR8anoCZusjOzeaqfM7JKhVYw\n65tZ9PhRj7ovONT1ZPRYBi3UyoexUfTv3J83B75J2c4ysjIC/dMx0ahqIG9AU2A+cEmEac2Bpu79\nfsDyWJZ5xhlnqEleT899WhmGPvnpk/WyvOWTZuhpt76ip945QZcffbLqzJk1ev2c0jm6aeemeilL\nqgCKNQD5Ee0WuENaABFpAEwEXlXVN8Onq+pWVd3u3p8CNBCRwxNcTJNgN5x+AwOOG1Av4+mFup5s\nOSSDlfyGrVd2r1E/u827NnPB6xdYnV6SCVzgiXPJ+ueBL1T1r1Hmae3Oh4h0w1mPjYkrpfGDtyGj\nLqHnrbN7YeqfadpoGz33PM9nk/4Z8zJyGuXw6iWvsmj9Igu9JBK4wAN+AlwDnOvpdtJPRG4UkRvd\neS4DFrt1eP8ABrm70ybFeUPvjul3sGrrqhq9PrzryTl/e5LZv1lA20M7sv3eO2s0LHuoTs9CL3nY\niMcmKZXvK2fRukWc0faMmF9TqetJ5mLn3Fj3MHZfxT4yP5wD8+axcegNHNb4sJiXO3nZZC4ZfwnD\nzhnGXT+9q2YrkmKCPuKxBZ5JeqP/M5qde3YytNvQqPPEOsTTiwte5Lfv/ZYZP5/B6W1Oj7kMC9Yu\n4JRWp5AhQTxoSpygB156fzsm6akqk5dP5papt0Q9I6Mm49kV5BXQvGFzeo7pyWdrPou5HKe1Po0M\nyaB0cymXjb/MDm8DygLPJDUR4fVLX496GlpNB+/My8lj9nWzaxV6AMs3Ludfy/5ldXoBZYFnkt5B\n596+6px+XdvBO8NDr2xHWcxl6XV0L2vICDALPJMS9ofeYT9hy8vP1XnwzlDoDe81nNwmuTV6rbXe\nBpcFnk+KiuDhh52/pn5kZ2Yz4aYPGXjzJAZNXwfbtx04XawWg3fm5eTxy9N/CcDHKz6u0eFtKPQa\nZTWq8fua+LHA80FREfToAfeMmM+5/TelTeglIuS/2bDTOTe28SGMdQfvrOtIxRVawY3v3ljjOr3+\nnfvzf7/4P1oe0pJde3exedfmOpXD1J0Fng9mz4bduoOKwf3YNbAnU2Zt8rtIcVdUBOeet4t7XppM\njx7xCb1Ig3fyzDM16kwcSYZkMGnwpFo1ZIgIqsoVE66g55iednjrMws8HxQUQENpQsaklyB3KeMb\n9WTTD6kderNnw+4uw6kYfD67TnyW2bPrd/nRBu9k/HgYOLDOoVeX1lsRYcjpQ6xOLwAs8HyQnw8f\nfAB/uq4Pf+32Nt/uXErPMakdegUF0LD4DljeDz1/CFs6PVtvy65q8E4KC53Qmzevzu/jDb2RxSNr\n9FpryAgGO9MiAKaVTGPA2AHcfubtPNzzYb+LEzdFRTB91i6mtbiUog1TGHX+KH51xq/qtEw/LpK9\nZtsacpvkkpWRharijmMRk9BpaIV5hUy7elocS+mPoJ9pYYEXEMWrizm11an+XZowgXbt3cWl4y9l\n7qq5LL9lOTmNcmq1HD/Czmvt9rVcOv5Snuz7ZI1OQ5u6fCrtmrfjlFanxLF0/gh64NkhbUB0aduF\nBpkNWL9jPVe9eVVKH942ymrExIET+ff1/07asAMnuFdtXVXjOr2+x/TllFanoKqMLB5ph7cJZIEX\nMEvWL2Hi0okpX6fXKKsRnQ/rjKpy/6z7eXZ+7HV6QQg7qIfT0DYt57Zpt1mdXgJZ4AVM4VGFvD3o\nbZaWpX5DBsA+3UfxmmKGvDskptALStiFhIfegrULYn5t58M6W0NGglngBVCfTn3SJvSyMrKYOHAi\n/Y7pV23oBS3sQkKhl98hnzZN29TotdZ6m1gWeAEVCj1F2b13t9/FiatQnV5VoRfUsAvJy8lj8pWT\nadW0FXv27WHZxmUxvzYUekvLlvLvFf+OYymNtdIG3L6KfWRmZLK3Yi/by7fXupI/Gezau4uBbwxk\n4DeNuTr/1/v70pWs38agp+bAnj2Mvb1H4MIu3K1Tb+WVz1+p8SCia7evpXXT1gA17u4SFNZKW0si\n0kdE/isiJSLyhwjTG4rIOHf6p+41bFNOZkYmADdMuoFzXzo3pQ9vG2U14p1B7zhhN3Agq6aNPxB2\nm7cw9qxmgQ87gN/m/7ZWDRmhsJu8bDJnv3C2Hd7GQSADT0QygRFAX+AEYLCInBA22w3A96raCXgC\neDSxpUyswScNTos6PRGBwkKKnrufzh/9jgF/neaEXa9WdLqgHi7AnQB1bb0Fp1+m1enVv2oDT0SG\nikjLRBTGoxtQoqpfq2o5MBa4KGyei4CX3PsTgB6SjMcAMUqnhgyA5l0H027vX2iwczdjW5QmTdiF\neEPvknGXUL6vPObXhjdkvDfnextKrJ7EsofXGpgnIuPdw8xEhEo7YIXn8Ur3uYjzqOpeYAtw0KWm\nRGSIiBSLSHFZWewj1waRN/QuHncxqVr/WrJ+G9eO/ITmOyt4o8V3dBr5RJ1P/vdDKPTGXTaO7Mzs\nGr02FHqfr11E39d6cc8D2+M2ykw6qTbwVPUe4Bici2NfBywXkYdE5Og4litSqIb/umOZB1Udpapd\nVLVLbm7NRq4Noj6d+vDOoHf487l/TspK7epUqrPr1YpOD95dbyOe+CEvJ4/u7bsD8NTcp2o8nt7g\nzDfh259SsasJ5eXU+ygz6SamOjz3Itdr3dteoCUwQUSGx6lcK4EOnsftgdXR5hGRLKAFkNrHea7e\nnXpz9pFnA/Dq56+ycedGn0tUP/Z3Pdmzp3KdXT2OeOKXbbu38djHj9W4Tu/Gc/vT6MMnyMwUsrOd\nUWdM7VXbLUVEbgWuBTYAzwFvq+oeEckAlqtqve/puQG2DOgBrALmAVeq6hLPPDcDJ6vqjSIyCLhE\nVQdWtdxk7JZSlRVbVnDMk8dwQu4JzPj5DA495FC/i1RrQe9nVx9KN5dS8GIBW3dvrVGXlaIiZ8+u\noMAZWizIUqFbyuE4YdJbVd9Q1T0AqloBnB+PQrl1ckOB94AvgPGqukREHhCRC93ZngcOE5ES4LfA\nQV1XUl10XTWcAAAOj0lEQVSHFh1SoiEjHcIOat96m58Pd90V/LBLBtbxOAWExtNLxj29dAk7r9LN\npRS+VMjdZ99d5/EAgyYV9vBMwIVab7/Y8AXTSpJnUMl0DDtw9vQW/7/F+8Mu1U8dDBILvBTRp1Mf\nlg1dxpUnXwlQqy4ribx0ZLqGXUiT7CYAfPTdR3R6slOtOiebmrPASyEdWjgN2x999xHnvHhOjer0\nioqg8Kpi7nnlrbj390r3sPNq37w9mZJZ6zMyTM1Y4KWg7eXbmbtqbo0aMmbPht1n3U/FpQPZ/aOJ\ncevvZWFXWX2chmZiZ4GXgmpzGlpBATR693VY3ZWKSwahx0+s93JZ2EUWHnolm0r8LlLKssBLUeGh\nV91V7/PzYebU5tz3o2mcdGhX7l80iIlL6y/0LOyqFgq96398PUflHOV3cVKWdUtJcdNKpjF28Vie\nu/A5sjKyYnrN1t1b6fNKH9o0a8OEyyfU+RQ2C7uaW7V1Fet3rOfHbX7sd1FqJOjdUizw0siabWto\nmNUwpn5623Zvo2FWQ7Izs/cPQlobFna1c97L51G8urjGg4j6LeiBZ4e0aWJfxT76vto35jq9Zg2b\nkZ2ZzYadG+j+XHfe/OLNGr+nhV3tjbpglDVkxIEFXprIzMjkkZ6P1Pg0tOzMbBpmNeSKCVfUKPQs\n7OrGWm/jwwIvjdSm9bZ5w+ZMvWoq3dp1izn0LOzqhzf0fvf+71J2/MNEssBLM97Qu3XqrTG9Jjz0\nJi+bHHVeC7v6lZeTx4fXfcgbl7+RkuMfJpoFXhrq06kP066exl97/zXm14RCb9BJgzit9WkR57Gw\ni4+OOR05vPHhlO8r59q3r7XD2zqwwEtTBXkFHNHkCMr3lfOHGX+I+fD25Ytfpl3zdux79BE+fefp\n/dP2j1S8Y7uFXZyU7Sjjw9IPrU6vDizw0tzCtQv52yd/q/F4eo+0Xs7Z82/mzTf+NykvpZiM2jVv\nZw0ZdWSBl+a6tutaq0FEbxn0BN1anshVnz/HxX99L+kupZisrPW2bizwTK1bb0dcPoP2e/5C1s5d\nSXkpxWQVCr2OOR2p0Aq/i5NULPAMcCD0Vm1bRenm0mrnL1m/jV+Omk+LJL+UYrLKy8lj/pD5dGnr\nnNSwbvs6n0uUHAIVeCLyFxH5UkQ+F5G3RCQnynylIrJIRBaISPqeK1bP+nTqw9e3fr3/VKZoI/Gm\n2qUUk1WGOD/fEXNHcPyI4+3wNgaBCjxgOnCSqp6Cc9Wyu6qYt1BVTwvyeXvJKDQS7+MfP07+8/kH\nHd6m8qUUk1X/zv2tTi9GgQo8VX3fvWIZwCc416M1PjjxiBNZunYRPZ/uvj/0KnU9ub3HwXV2hYVw\nxx0+lDa9WUNG7AIVeGGuB6ZGmabA+yIyX0SGVLUQERkiIsUiUlxWVlbvhUxVfTr14e2T/sTSzSX0\nfLo7xd+ttK4nAeYNvV4v9+L7H773u0iBlPDhoURkBtA6wqQ/quo77jx/BLrgXA/3oAKKSFtVXS0i\nR+AcBt+iqnOqe+90Hx6qNqZNfJRL//MUbfYOJ2enWteTgCvdXMrHKz7efzGnRAv68FCxjQhZj1S1\nZ1XTReRanAt894gUdu4yVrt/14vIW0A3oNrAMzXX6ac3cfS8o8ncuYuxLb6j0wX+/JBMbPJy8sjL\nyQNgxtfONYqTaTy9eAvUIa2I9AHuBC5U1Z1R5mkiIs1C94HzgMWJK2X6CNXZZe8st64nSWZvxV5u\nmXqL1emFCVTgAU8BzYDpbpeTkeAcworIFHeeVsBHIrIQmAtMVtXkufp0krCuJ8ktKyOLqVdNtYaM\nMDbEuznI/q4nO7Y7DRTeOrtZs5yuJ9YamxRKN5dS8GIBW3dvTchw8UGvwwvaHp7xWaUhnqzrSdLz\ntt6OWTjG7+L4LuGNFia4bDy71JSXk8env/yU3Ca5AKhq2g4mant4BrCwS3WtmrYiQzL4dvO3dH22\na9rW6VngGQu7NKIoG3ZuSNuGDAu8NGdhl17CT0N78b3PePhhKCryu2SJYYGXxizs0lMo9BrSnF/M\n6sk9T35Ojx7pEXoWeGnKwi695eXkcc2+2VBaSMX3HSgvh9mz/S5V/FngpSELOwNwcWEeh0yaSOae\nlmRnQ0GB3yWKP+uWkmYs7ExIfj588IGzZ1dQ4DxOdRZ4acTCzoTLz0+PoAuxQ9o0YWFnjAVeWrCw\nM8ZhgZfiLOyMOcACL4VZ2BlTmQVeirKwM+ZgFngpyMLOmMgs8FKMhZ0x0VngpRALO2OqFrjAE5Fh\nIrLKvabFAhHpF2W+PiLyXxEpEZE/JLqcQWNhZ0z1gnqmxROq+li0iSKSCYwAegErgXkiMklVlyaq\ngEFiYWdMbAK3hxejbkCJqn6tquXAWOAin8vkCws7Y2IX1MAbKiKfi8hoEWkZYXo7YIXn8Ur3uYOI\nyBARKRaR4rKysniU1TcWdsbUjC+BJyIzRGRxhNtFwDPA0cBpwBrg8UiLiPBcxOtNquooVe2iql1y\nc3PrbR38ZmFnTM35Uoenqj1jmU9EngXejTBpJdDB87g9sLoeipYULOyMqZ3AHdKKSBvPw4uBxRFm\nmwccIyJHiUg2MAiYlIjy+c3CzpjaC2Ir7XAROQ3nELUU+DWAiLQFnlPVfqq6V0SGAu8BmcBoVV3i\nV4ETxcLOmLoJXOCp6jVRnl8N9PM8ngJMSVS5/GZhZ0zdBe6Q1hzMws6Y+mGBF3AWdsbUHwu8ALOw\nM6Z+WeAFlIWdMfXPAi+ALOyMiQ8LvICxsDMmfizwAsTCzpj4ssALCAs7Y+LPAi8ALOyMSQwLPJ9Z\n2BmTOBZ4PrKwMyaxLPB8YmFnTOJZ4PnAws4Yf1jgJZiFnTH+scBLIAs7Y/xlgZcgFnbG+M8CLwEs\n7IwJBgu8OLOwMyY4AjXEu4iMA451H+YAm1X1tAjzlQLbgH3AXlXtkrBC1oCFnTHBEqjAU9UrQvdF\n5HFgSxWzF6rqhviXqnYs7IwJnkAFXoiICDAQONfvstSGhZ0xwRTUOryfAutUdXmU6Qq8LyLzRWRI\nAstVLQs7Y4Ir4Xt4IjIDaB1h0h9V9R33/mDg9SoW8xNVXS0iRwDTReRLVZ0T5f2GAEMAjjzyyDqU\nvHoWdsYEm6iq32WoRESygFXAGaq6Mob5hwHbVfWx6ubt0qWLFhcX172QEVjYGQMiMj+ojYgQzEPa\nnsCX0cJORJqISLPQfeA8YHECy3cQCztjkkMQA28QYYezItJWRKa4D1sBH4nIQmAuMFlVpyW4jPtZ\n2BmTPALXSquq10V4bjXQz73/NXBqgosVkYWdMckliHt4ScHCzpjkY4FXCxZ2xiQnC7wasrAzJnlZ\n4NWAhZ0xyc0CL0YWdsYkPwu8GFjYGZMaLPCqYWFnTOqwwKuChZ0xqcUCLwoLO2NSjwVeBBZ2xqQm\nC7wwFnbGpC4LPA8LO2NSmwWey8LOmNRngYeFnTHpIu0Dz8LOmPSR1oFnYWdMeknbwLOwMyb9pGXg\nWdgZk57SLvAs7IxJX74EnohcLiJLRKRCRLqETbtLREpE5L8i0jvK648SkU9FZLmIjBOR7Fjed/fe\nCgs7Y9KYX3t4i4FLgEoXzxaRE3CuWnYi0Ad4WkQyI7z+UeAJVT0G+B64IZY3/bpsO2BhZ0y68iXw\nVPULVf1vhEkXAWNVdbeqfgOUAN28M4iIAOcCE9ynXgIGxPreFnbGpK+gXaaxHfCJ5/FK9zmvw4DN\nqrq3inn2E5EhwBD34e5jWjX39aLdcXI4sMHvQsRBqq4XpO66Het3AaoSt8ATkRlA6wiT/qiq70R7\nWYTntBbzHJigOgoY5ZapWFW7RJs3Wdl6JZ9UXTcRKfa7DFWJW+Cpas9avGwl0MHzuD2wOmyeDUCO\niGS5e3mR5jHGmIMErVvKJGCQiDQUkaOAY4C53hlUVYFZwGXuU9cC0fYYjTFmP7+6pVwsIiuBfGCy\niLwHoKpLgPHAUmAacLOq7nNfM0VE2rqLuBP4rYiU4NTpPR/jW4+qx9UIEluv5JOq6xbo9RJnh8kY\nY1Jf0A5pjTEmbizwjDFpI+UDr66nsSULERkmIqtEZIF76+d3mepCRPq430uJiPzB7/LUFxEpFZFF\n7ncU6C4c1RGR0SKyXkQWe547VESmu6d9TheRln6WMVzKBx51P40tmTyhqqe5tyl+F6a23O9hBNAX\nOAEY7H5fqaLQ/Y6SvR/eizi/Ha8/AB+4p31+4D4OjJQPvLqcxmZ80w0oUdWvVbUcGIvzfZkAUdU5\nwKawpy/COd0TanjaZyKkfOBVoR2wwvO4ylPUksRQEfncPdQI1KFEDaXidxOiwPsiMt897THVtFLV\nNQDu3yN8Lk8lQTuXtlbieBpboFS1nsAzwIM46/Ag8DhwfeJKV6+S7rupgZ+o6moROQKYLiJfuntK\nJgFSIvDieBpboMS6niLyLPBunIsTT0n33cRKVVe7f9eLyFs4h++pFHjrRKSNqq4RkTbAer8L5JXO\nh7TVnsaWTNyNK+RinMaaZDUPOMYd6DUbp3Fpks9lqjMRaSIizUL3gfNI7u8pkkk4p3tCAE/7TIk9\nvKqIyMXAk0AuzmlsC1S1t6ouEZHQaWx78ZzGlqSGi8hpOId+pcCv/S1O7anqXhEZCrwHZAKj3dMO\nk10r4C1nSEeygNdUdZq/Rao9EXkdKAAOd08VvR94BBgvIjcA3wGX+1fCg9mpZcaYtJHOh7TGmDRj\ngWeMSRsWeMaYtGGBZ4xJGxZ4xpi0YYFnjEkbFnjGmLRhgWcCR0S6uoMgNHLPTlgiIif5XS6T/Kzj\nsQkkEfkT0Ag4BFipqg/7XCSTAizwTCC559DOA3YBZyX5aX8mIOyQ1gTVoUBToBnOnp4xdWZ7eCaQ\nRGQSzkjHRwFtVHWoz0UyKSDlR0sxyUdEfg7sVdXX3OtbfCwi56rqTL/LZpKb7eEZY9KG1eEZY9KG\nBZ4xJm1Y4Blj0oYFnjEmbVjgGWPShgWeMSZtWOAZY9LG/wevv5SVVEUupAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x178e8aada20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = ((np.random.rand(10)*2)-1)*10\n",
    "y = ((np.random.rand(10)*2)-1)*10\n",
    "pair = zip(x, y)\n",
    "\n",
    "# The projection matrix\n",
    "P = np.array([[0.5, 0.5],\n",
    "              [0.5, 0.5]])\n",
    "\n",
    "for item in pair:\n",
    "    x, y = item\n",
    "\n",
    "    # Form the vector z in R2 from the point (x, y)\n",
    "    z = np.array([[x],\n",
    "                  [y]])\n",
    "\n",
    "    # Project z onto the line y = 0\n",
    "    Pz = P @ z\n",
    "\n",
    "    # Plot the result\n",
    "    plt.plot(x, y, 'b.')\n",
    "    plt.plot(Pz[0], Pz[1], 'rx')\n",
    "    plt.plot([x,Pz[0]],[y,Pz[1]],color='green', linestyle='dashed')\n",
    "\n",
    "plt.plot([-10,10],[-10,10])\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.axis('square')\n",
    "plt.ylim(-10, 10)\n",
    "plt.xlim(-10, 10)\n",
    "plt.title(\"Projection of 10 Random Coordinates onto y=x\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3\n",
    "\n",
    "In this third example, we again project onto the line $y=0$ but in some direction defined by a variable $\\alpha$. This creates an oblique projection. The projection matrix is\n",
    "\n",
    "$$\n",
    "P=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & \\alpha \\\\\n",
    "   0 & 0 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$\n",
    "\n",
    "We can prove this is an operator by showing that $P^2 = P$:\n",
    "\n",
    "$$\n",
    "P^2 = P \\cdot P =\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   0 & 0 \\\\\n",
    "   \\alpha & 1 \\\\\n",
    "  \\end{array} } \\right]\n",
    "  \\cdot \\left[ {\\begin{array}{cc}\n",
    "   0 & 0 \\\\\n",
    "   \\alpha & 1 \\\\\n",
    "  \\end{array} } \\right]\n",
    "  = \\left[ {\\begin{array}{cc}\n",
    "   0 \\cdot 0 + 0 \\cdot \\alpha & 0 \\cdot 0 + 0 \\cdot 1 \\\\\n",
    "   \\alpha \\cdot 0 + 1 \\cdot \\alpha & \\alpha \\cdot 0 + 1 \\cdot 1 \\\\\n",
    "  \\end{array} } \\right]\n",
    "  = \\left[ {\\begin{array}{cc}\n",
    "   0 & 0 \\\\\n",
    "   \\alpha & 1 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$\n",
    "\n",
    "The projection is an orthgonal projection only if $\\alpha = 0$.\n",
    "\n",
    "In the code sample below, you can play with the parameter ``a`` to see how it affects the direction of projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAEWCAYAAAD7MitWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFPWZ+PHPMyeHDOeIyMg9w6WCwiqjIiMoV6ZFjUYh\nahLdoElMsrua9Uiy8ENXEpPobmIiUYNJvPFAAUFFcDCugwqKCnILCIIccl8Dwzy/P6qGNGPP1VdV\ndz3v16tf09VV/f0+1VX1zLeOb5WoKsYYEwQZXgdgjDHJYgnPGBMYlvCMMYFhCc8YExiW8IwxgWEJ\nzxgTGCmb8ETkLhF5NFXKbUC9l4vIRhHZLyJnJbv+RBGR74rI217HkQwioiLSw30/RUR+6XVM5kRJ\nTXgisl5EDrkb9VYReUxEToqmLFW9V1X/NcZ4SkRkU7zLjdJvgVtU9SRV/bDmSBG5W0Q+EZFKEZkY\nYfw4EdkgIgdE5CURaVNbRe6GecBdDl+IyP0ikhnf2Uk+EckRkYkistqdv/UiMlVEuiQ7FlW9WVXv\njrUcEeniLq+seMQVL+5ve7HHMbQRkenust4gIuPq+44XLbyQqp4EnA38C/CLmhOII2Vbn1HqDCyr\nY/wa4D+BV2qOEJG+wJ+B64D2wEHgT/XU189dDkOAq4EboojZb54HLgXGAS2BfsBiYFg8K/Fb8gmw\nPwJHcNb5bwMPudtC7VQ1aS9gPXBx2PBvgFnu+zLgv4H/Aw4BPYBTgRnATpwN/vth350IPBE2PAh4\nB9gNfASUhI1rAzwGbAZ2AS8Bzd16qoD97uvUCOVeipOIdrsx9q4xP7cBHwN7gGeBJrXMewZOct8A\nbAP+jrNR5rp1K3AAWFvPb/gEMLHGZ/cCT4UNd3dXhBa1lKFAj7DhacAfw4a/BywH9gGfATeFjSsB\nNgG3uvOxBfhe2Pi27jLbC7wH3A28HTb+POB99/d6HzgvbFwZcI+7HPcDM93ynnTLex/oUss8Xewu\nz9Pq+O3qWp9ygf9x15HN7vvcGvN8O/Al8Lj7+c/c+d+M8w/j+O8K/BW4p4G/2TeAD9153Bi+fIHP\n3XKr19Fi9/Mb3GW0C3gN6Ox+LsADbj17cNbN06P4PSa668Xf3fVgGTDQHfc4znZzyI3pP+vbVmrU\n+0fgdzU+mwn8WyNySXOcdbwo7LPHgV/V+b1EJLY6glyPm/CA09wf5+6wlf1zoC+QBWQDC3BaKk2A\n/sB2YFjYAnnCfd8R+AoYjZNYLnGH893xr+Ako9ZuuUPCV8QaMYaXW4SThC5xv/ef7oqREzY/77kr\nTht3Bby5lnm/wf1uN+Ak4EXcDSdSEqrjN4yU8F4Gbq/x2X5gQC1lhG+YvXA2wH+vsQF2x9l4huC0\nGM8O+80qgUnubzLaHd/aHf8MzobSHDgd+AI34bm/0S6clmgWMNYdbhu2Dqxx624JfAqswklmWTgb\n32O1zNOvgAX1/HZ1rU+TgIXAyUA+TtK9u8Y8/xonMTYFRgJb3XlsDjxF3Qmvrt+sBDgDZ9090y33\nMndcF7fcrLD5uMz9nXq7v8svgHfccSNwWrWt3OXXG+gQxe8xETjsxpoJTAYWRtqWG7Kt1Kj3HJx/\nEhnucDv392jvDs/CSZqRXtUNpLOAQzXKvQ2Y6beEt98NfIP7YzcNW9knhU17GnCMsFaK+6P/NUJi\nup2w5OF+9hrwHaADzn+j1hHiKaHuhPdLYFrYuAycDbgkbH6uDRt/HzCllnmfB/wwbLgncLR6RSa2\nhDePGok2PM4IZShOa+KA+/5p3NZMLdO/BPw07Dc7xIkb4DacFnamO0+9wsbdyz8T3nXAezXKLge+\nG7YO/Dxs3O+AOWHDIWBJLTE+AjxTxzzUtz6tBUaHjRsBrA+b5yOEtd6BqYS1JnA2+LoSXsTfrJZY\n/wd4wH3fha8nvDnAjTXWy4M4h0WG4vyTGISbUKL8PSYCb4SN60NYguHrCa/ObSVC/cuBS9z3twCz\n61v3a3x/MPBljc++D5TV9T0vjpNdpqqtVLWzqv5QVQ+FjdsY9v5UYKeq7gv7bANOa66mzsBVIrK7\n+gVcgJPsTnPL2RVFrKe6dQKgqlVujOExfBn2/iBO663estz3WTjHH2K1H8ir8Vkezq5Ibc7GifVq\n4FycVgoAIjJKRBaKyE73txyN81+42leqWhk2XD3f+TjzFL4cw+e55m9QPT7899wa9v5QhOHaft+v\ncJZ3bepbnyItn1PDhrer6uEa5dU2nxHjq+U3Q0TOFZE3RWS7iOwBbubE37umzsD/hq3rO3Facx1V\ndT7wIM5u41YReVhEaq4b1fHXt33VXLeb1HH8siHbSri/Ade676/F2R1tjGjWed9dlqJh7zcDbUSk\nRdhnnXD+a9S0EaeF1yrs1VxVf+WOayMireqpL5LNOCsX4JxMwUmgkWKozwll4cxLJSdu0NFahnOA\nHgAR6Yaz67Wqri+pYxpOK+u/3O/mAi/gnDVur6qtgNk4G1R9tuPM02lhn3UKe1/zN6geH83vWdMb\nwDkiUlDL+PrWp0jLZ3PYcM11ZQu1z2djPYVzLO00VW0JTOGfv3ekdXQjznHV8PW9qaq+A6Cqv1fV\nATiHh4pwjjXW1JjtK5KacTV2W3kCGCMi/XB2u18K++4c9wqCSK857mSrgCwRKQwrsx91n/jzXcI7\nTlU34hxHmSwiTUTkTOBGnAPYNT0BhERkhIhkutOXiEiBqm7B2QX4k4i0FpFsEbnQ/d5WoK2ItKwl\njGnAN0RkmIhk4xx0rnDjaqyngX8Xka7upTj3As/W+K9fKzfuJjjLLMudx+pLSZ5053+wiDTHOVb0\nYo3/3nX5FTBeRE4BcnCS5XagUkRGAcMbUoiqHsM5NjlRRJqJSB+cwwrVZgNF7iU0WSJyNc6u0qwG\nxllX3W8Ac4HpIjLALb+FiNwsIjc0YH16GviFiOSLSDucfwBP1FHlNOC7ItJHRJoBE2IIvwVOa+uw\niJyDc5a52nacQzLdwj6bAtxZfUZSRFqKyFXu+39xW4zZOIcsDuPsup6gkdtXJFtrxNSobUVVN+Gc\nhHoceCF8T09VR6lzeVak1yh3mgM469okEWkuIucDY6inpejbhOcai3MMYzMwHZigqnNrTuQuvDHA\nXTgryEac/2rV83cdzrGlFTjHTv7N/d4KnBX9M3f34NQa5a7EaW7/AdiBcwwppKpHopiXqTgL4y1g\nHc6K+ONGfP8RnF26scDP3ffXuXEuw9kNetKdvxbADxtasKp+gnMA+2dukvwJzgq8C2fjm9GIOG/B\n2VX7Euc41mNh9XwFlOJsDF/hHNguVdUdjSi/LlfiJNVncc5QLgUG4rT+oO716R5gEc5ZzU+AD9zP\nIlLVOTjH2ubjHJyfH0PcP8TZcPfhJNppYfUcxL16wV1HB6nqdJwTKM+IyF53Pke5X8nDWVd24exi\nfoXTWo+kQdtXLSbj/IPYLSK3Rbmt/A3nZE1jd2er/RDnBNI2nO34B+62UCtxD/alHBGZBBSoajpc\nP2ZM4Lh7Wk/gXGpUlYw6/d7Ci8g9PtAHp6VkjEkx7m7vT4FHk5XswOOE53b72SYiS8M+ayMic8Xp\nHjRXRFpH+OoHONcNfc+d7jsRpjHG+JCI9Ma5NK0DzmGB5NXt5S6t26TdD/xdVU93P7sP5wDur0Tk\nDpzr526v8b02OMdbBuKcLVqMc5FtNJeeGGMCwtMWnqq+hXMNUbgxOAczcf9eFuGrI4C5qlp9fd1c\nnCvfjTGmVn7sBN3evZQEVd0iIidHmKYjJ170uYlaLnAUkfHAeIDmzZsP6NWrV5zDNcZUW7x48Q5V\nzfc6jtr4MeE1RKSLYCPum6vqw8DDAAMHDtRFixYlMi5jAk1E6utx4ik/nqXdKiIdANy/2yJMs4kT\nr3Iv4MSr4o0x5mv8mPBm8M+r87+DcyeQml4Dhrs9J1rj9AR4LUnxGWNSlNeXpTyN04+zp4hsEpEb\ncbo5XSIiq3FuNfMrd9qB4t56XVV34txn7X33Ncn9zBhjapWyPS2iYcfwjEksEVmsqgO9jqM2ftyl\nNY1UXg6TJzt/jTG1S9WztMZVXg7DhkFF1SFy727KvHlQXOx1VMb4k7XwUlxZGVR0fYmq29pQ0Xwt\nZWVeR2SMf1nCS3ElJZC980zIPkxmr1mUlHgdkTH+ZQkvxRUXw5svduNk6UO/b8203Vlj6mAJLw0U\nF8N3i0tZsnsBew7v8TocY3zLEl6aCPUMUVlVyWtr7fprY2pjCS9NFBcUc+/QexnQYYDXoRjjW3ZZ\nSprIzMjkzsF3eh2GMb5mLbw0UlFZwYyVM1i5Y6XXoRjjS5bw0sjhysN8c9o3mfrhVK9DMcaXLOGl\nkZZNWjKk8xBmrprpdSjG+JIlvDQTKgqxfMdy1u5c63UoxviOJbw0U1pUCsCsVbM8jsQY/7GEl2a6\nt+lO73a9eWfTO16HYozv2GUpaWj+d+bTvnl7r8Mwxncs4aWhU046xesQjPEl26VNU7fPvZ1bX7vV\n6zCM8RVfJjwR6SkiS8Jee0Xk32pMUyIie8Km+S+v4vWjLfu38NeP/kplVaXXoRjjG75MeKq6UlX7\nq2p/YABwEJgeYdJ/VE+nqpOSG6W/hYpC7Dy0k/KNdt93Y6r5MuHVMAxYq6q+fsCv34zoMYKsjCy7\nPMWYMKmQ8K4Bnq5lXLGIfCQic0SkbzKD8ru83DzrdWFMDb5OeCKSA1wKPBdh9AdAZ1XtB/wBeKmW\nMsaLyCIRWbR9+/bEBetD1/e7niGdh3Dk2BGvQzHGF3z9XFoRGQP8SFWHN2Da9cBAVd1R2zT2XFpj\nEsueSxubsdSyOysip4iIuO/PwZmXr5IYW0qo0iq7XZQxLt8mPBFpBlwCvBj22c0icrM7eCWwVEQ+\nAn4PXKN+bq565O4Fd3P6Q6fbsy6MwccJT1UPqmpbVd0T9tkUVZ3ivn9QVfuqaj9VHaSq1nk0gmHd\nhsXlWRfl5TB5svPXmFRlXcvSXHFBMW2btmXmqpl8q++3oiqjvByGDYMjRyAnB+bNwx4HaVKSb1t4\nJj4yMzIZXTia2atnR93roqwMKlqs5Nhl46hovpaysriGaEzSWMILgFh7XZSUQHZGNpzxNJm9Z1JS\nEtfwjEkaS3gBMKLHCGaNncXAU6O7WqC4GN58sRsnSx/OvGqm7c4GVDocx7WEFwB5uXl8o+gbNM1u\nGnUZxcXwvfNCfLT7LTvjG0Dl5XDRlav4xRPTGTYsdZOeJbyA2LxvMxPenMD63eujLiNUFKKyqpJX\n17wav8BMSigrgyMD7qfqsuuoOHY4ZY/jWsILiENHDzHprUm8vOLlqMsYVDCI4oJiqrQqjpGZVFBS\nAtmfhSDnAFk9FqTscVxLeAFR/ayLWG4mkJmRyTs3vsPYM8bGMTKTCoqL4bUpQ8mmKaFbU/c4riW8\nAAkVhViwYUHMx+COHjvK3oq9cYrKpIqSC5oyquclvL93JqnaqckSXoCEeoZi7nVx5NgRCh4o4J63\n7oljZCZVlBaWsnnfZtbtXud1KFGxhBcgxQXFdGrZiS/2fhF1GTmZOZzZ/ky7z15AjT1jLDt+toNu\nrbt5HUpUrGtZgGRmZPLZTz4jMyMzpnJCRSF++upPWbNzDT3a9IhTdCYVnJRzktchxMRaeAFTnexi\nOdNaWlQKYLePD6h3Nr7DhY9dyLYD27wOpdEs4QXMsapjnPvoudw1766oy+jWuht98vvYbm1ANc1q\nyj8+/wezV8/2OpRGs4QXMJkZmbTIacGMlTNiKuc3l/yGiUMmxicok1L6n9KfgryClPyHZwkvgEJF\nIZbvWM7anWujLmN04WgGdx4cx6hMqhARSgtLeX3t61RUVngdTqNYwgugUM8QEPsxuPKN5Tz58ZPx\nCMmkmNKiUvYf2U/Z+jKvQ2kUS3gBFK9jcFMWT+HHc34c9X32TOoa2nUoI3uMJCczx+tQGsUuSwmo\nO86/I+Y+saGiEH//6O/OWbvOF8YpMpMKmmY3Zc6353gdRqP5NuG5j13cBxwDKms++s19Ytn/AqOB\ng8B3VfWDZMeZqq7rd13MZQzvPpzsjGxmrZplCS+gvjroPCiwbbO2HkfSMH7fpb1IVfvX8pzLUUCh\n+xoPPJTUyNLApr2bmL9uftTfz8vNY0iXISl5ts7EbuehnbT/bXv+vPjPXofSYH5PeHUZA/xdHQuB\nViLSweugUsld8+7iqueuiukYXKgoxJf7v2T7ge1xjMykgjZN23BWh7NS6h+enxOeAq+LyGIRGR9h\nfEdgY9jwJvezE4jIeBFZJCKLtm+3jTJcaVEpOw/tZOGmhVGX8f2zv8+227aR3zw/jpGZVBEqCvHu\npndTpteFnxPe+ap6Ns6u649EpOZBIonwna/ds0ZVH1bVgao6MD/fNspwI7qPICsji5kro/8P3TS7\nKdmZ2XGMyqSSUFEIRVOm14VvE56qbnb/bgOmA+fUmGQTcFrYcAGwOTnRpYeWTVoypHPsx+Bmr55N\n3z/1Zffh3XGKzKSKVOt14cuEJyLNRaRF9XtgOLC0xmQzgOvFMQjYo6pbkhxqyqvudbFp76aoy2jV\npBWfbv+U19ZEf589k5pEhL9d9jfuu/g+r0NpEF8mPKA98LaIfAS8B7yiqq+KyM0icrM7zWzgM2AN\n8AjwQ29CTW3X97uezf+xmYK8gqjLOLfjubRr1i5l/sub+BradSjd23T3OowG8eV1eKr6GdAvwudT\nwt4r8KNkxpWOWjdtHXMZmRmZjC4czcyVM6msqiQrw5erlUmg55Y9x/4j+70Oo15+beGZJFq4aSGh\np0MxPesiVBRi1+FdlG9M0QeWmpg8/vHjTFww0esw6mUJz1BZVcmsVbNietbF8O7DGXfGuJS/I66J\nTqgoxOd7Pvc6jHpZwjMUFxTTtmnbmI7B5eXm8eQVT3JWh7PiGJlJFdV3wfY7S3jm+DG42atnx3zn\nk1VfrWLHwR1xisykig4tOjDw1Eg9QP3FEp4BnF2SWHtdbNi9gZ4P9uSJj5+IY2QmVZQW+r+VZwnP\nAM4xuEEFgzhceTjqMjq36mzPugiwn1/4c69DqJclPAM4vS7Kbyzn4m4Xx1ROqCjEWxveiumMr0lN\nqXA5kiU8c4KDRw+yr2Jf1N8PFYWorKqM6YyvMYliCc8ct+3ANtre15a/fPiXqMsYVDDIel0Y37KE\nZ447ufnJdGvdLaZklZmRycyxM/n9yN/HMTJj4sMSnjlBPI7BDSoYFJcua8bEmyU8c4LSotK4HIP7\nw7t/4JHFj8QpKmPiwxKeOUE8el0AzFg1g/sX3h+nqIyJD0t45gSZGZk8eumj/Oy8n8VUTqgoxIod\nK1izc02cIjMmdpbwzNdc1usyzmx/ZkxlVPetnLVqVjxCMiYuLOGZiGavns2Ly1+M+vvdWnezXhfG\ndyzhmYgeWPgAv5j/i5jKuKLXFWRKJlVaFaeojImNJTwTUWlhKct3LGftzrVRlzHpokm8ft3rZIit\nZsYffLcmishpIvKmiCwXkWUi8tMI05SIyB4RWeK+/suLWNNZqGcIIKZdUhHnSZoVlRVxicmYWPku\n4QGVwK2q2hsYhPNM2j4RpvuHqvZ3X5OSG2L6i9cxuMn/mEyn/+kU8332jIkH3yU8Vd2iqh+47/cB\ny4GO3kYVTKGiEOt2rePosaNRl1HYtpBtB7bZsy6ML/gu4YUTkS7AWcC7EUYXi8hHIjJHRPrWUcZ4\nEVkkIou2b9+eoEjT04QhE1j7k7VkZ2ZHXcbw7sPJzsi2s7XGF3yb8ETkJOAF4N9UdW+N0R8AnVW1\nH/AH4KXaylHVh1V1oKoOzM/PT1zAaahpdtPjx+GilZebx5AuQwKR8MrLYfJk56/xJ18mPBHJxkl2\nT6rq1y4GU9W9qrrffT8byBaRdkkOMxCmfjiVMx46g2NVx6Iuw8teF8lKQuXlMHT0Ln5x/xqGDbOk\n51e+u0WpOE2KvwDLVTViZ0wROQXYqqoqIufgJO6vkhhmYJyUcxJLty2lfFM5F3S6IKoyLut1GUeP\nHaVlbss4R1e38nIYNgwqcraQe3cH5s2D4uLE1FVWBoevLYad3Tny7CuUlSWuLhM9P7bwzgeuA4aG\nXXYyWkRuFpGb3WmuBJaKyEfA74FrVFW9Cjidjeg+gqyMLGaujH6XtFPLTtx63q3kN0/uIYWyMjh8\nzj1U/bgrFXqAsrLE1VVSApnrR0K3eWQ3P0BJSeLqMtHzXcJT1bdVVVT1zLDLTmar6hRVneJO86Cq\n9lXVfqo6SFXf8TrudNWySUuGdI79GNzeir089clTSX3WRUkJZH95HmRVkFX4RkKTUHExPHBTCLIq\nmPj3N6x151O+S3jGf0JFoZh7XXyy9RO+/eK3eXXNq3GMrG7FxfDGXwaTSx4jfjIz4UnoplGDycvN\nY43YDRP8yhKeqdelPS/lpgE3xXTG1qtnXQw+P5sxfUfy3u5ZCe/Tm5OZw8geI5m1OvF1mehYwjP1\n6tq6K1NKp9Ctdbeoy8jMyGR04Whmr56d9F4XoaIQWw9sZdHmRQmva0zLCYyrWMC7C23T8iNbKqZB\nqrSK9794P+ZHOO46vCvpvS6+UfgNnrziSXq165XQesrL4Xu3ruGBZ5bYpSk+ZQnPNMjCTQs559Fz\nmLNmTtRlVPe6ePvzt+MYWf1aN23NuDPGkZebl9B6ysrgyFl/QIfdxuEL7kjoWWETHUt4pkHO7Xgu\nbZq2iekYXF5uHut+uo47B98Zx8gaZvuB7dxffj+b921OWB0lJZD9WQhabUTP/zV9Bm1KWF0mOpbw\nTIPE6xhcxzxv7gOx4+AObn39Vl5e8XLC6iguhqcnho4Pb2lhZ2v9xhKeabBQUYidh3aycNPCqMuo\nqKxg3AvjePSDR+MYWf16tetF99bdE36W+JvDutKnXR+aZDUJRP/hVGMJzzRYda+LWB7Mk5uVy8db\nP+bppU/HMbL6iQilRaXMXzefA0cOJLSuS3teSm5mLvM+m5fwukzjWMIzDdaySUve+u5bTBgyIaZy\nQkUh3trwVlJ7XVTXW3GsgrmfzU1oPRNKJvD8Vc/TuVVnNuzZkNC6TONYwjONUnxaMU2zm8ZURmlR\nKZVVlUntdQEwuPNg2jRtw8odKxNaT5OsJgzrNowVP1pBn/xIN+s2XrGEZxrlWNUxJrw5gWeXPht1\nGV71usjJzGHTv2/i9gtuT3hdUz+cSp8/9eHQ0UPW68JHLOGZRsnMyOT55c/zyAePxFTGTQNuomfb\nnnGMrGFibZ02VMsmLVmxYwUn//ZkFm9enJQ6Tf0s4ZlGKy0sZcGGBTEdg7tn6D38csgv4xhVwxyr\nOsYlj1/CpAWJfe7T8O7DyZIsDhw5YGdrfcQSnmm0UM8QlVWVvLb2tZjKOXrsKOt3r49PUA2UmZHJ\n4crDTF8xPaH15OXmUdK1xC5P8RlLeKbRiguKadu0bcwb8lXPXcWoJ0fFKaqGCxWFWPLlEjbu2Zjw\neg5VHkpKXaZhLOGZRsvMyOSqPleRlRHbEwKGdh3qybMuQkVOb4hXVr+S0Hou7Xkp404fl5S6TMNY\nwjNReaj0IR4b81hMZVQnnlguZI5GsnpddGnVhSeueIJJJZMoLrBbIPuBbxOeiIwUkZUiskZE7ogw\nPldEnnXHv+s+wzY47rsPbroJ3nzzn5+9+abz2X33Jaf+N9/k4NGDJ9bfiLq7PvIcfZt1OTHxNLKM\naMhvfsOtbUsZ2X1kYuu97z70pvFc8WXrf95LMJZl5PUyTweq6rsXkAmsBboBOcBHQJ8a0/wQmOK+\nvwZ4tr5yBwwYoGlj/nzVvDzVli2d9zWHk1D/D77ZRHvf1/mf8bRr17i658/XO0JNNev/ZeruQ7uj\nKyMaNetJVL3z5+snXZrrlubo/X+8Tj+e+Whsy8jrZd4AwCL1QQ6p7SVaz8O+ROQWnOfD7kpU0o1Q\nZzEwUVVHuMN3Aqjq5LBpXnOnKReRLOBLIF/rmKE2nXvrJXdNTWzwybR7NyxdCqpUCWSIQN/ToVWr\npFT/xdbVrDm8mUHakdzN26BPn0bXfXjnVg6vW03Lth2RzVuiKiMqu3dTuWIZFafk03zzjsTVu3sX\nS7d8zO4m0P4AFO7JjG0ZhS1zADKSu8zrM+3m8xar6kCv46hNQ3ZpTwHeF5Fp7m5mbI+ib5iOQPhp\nrU3uZxGnUdVKYA/QtmZBIjJeRBaJyKKjR48mKFyPtGqFFnTk3VOr+LxFFXQsSOqK3z6/K8VaQO6G\nL+DUU6Oqu0mb9rRqW4Bs+DzqMqLSqhVHTzmZ5hu2JLbeVq3p3vRUzt4ChV8R+zJq1QoKCqCqynkl\neZmnvIY0AwEBRgDPAGuAe4HuiWp2AlcBj4YNXwf8ocY0y4CCsOG1QNu6yk2rXVrV47s0F9woeum4\njOTv2lTvCv7yl9HvEsajjGgkq97q3c5mzVSbNo19GbnlVTVrqou65GpF6zzf7M6q+n+XtjFJqB/w\nP8AK4CHgQ+C+hAQFxcBrYcN3AnfWmOY1oNh9nwXsAGcXvbZXWiW86hW/ZUsd/q9N9PtTvpH0Y3gx\nHwdL1rE0r+qN9zG3sO+P+2VfZSI69/RmdgyvEa96d2lF5Ccishi4D/g/4AxV/QEwAPhmdO3Ker0P\nFIpIVxHJwTkpMaPGNDOA77jvrwTmuz94MLz/PlxzDTJ9OnnDS5m1/wOqpr8IV1/tjEtG/dOmwUUX\nOcMXXeQMN6bueJQRjWTV6y4jpk936rjoInjppeiXUVh5+cUXA/DCj4clb5mng/oyIjAJ6FzLuN6J\nysTAaGAVzq7qz8NiudR93wR4DmcX+z2gW31lplULL8zflvxNmYgOf3y416GYJFmwfoEyET35Nydr\nVVWV1+Ech89bePWepU0nAwcO1EWLEv9s0mTbcXAH+b/JJ1My2XPHHprnNPc6JJNglVWVtJzckoOV\nB1n6g6X0Pbmv1yEBICIpf5bW+Fy7Zu3om9+XY3qMNz57w+twTBJkZWQxvMdwAF5embgHE6UbS3hp\n4vHLH6dFTgu7M0eA3Fp8K7cV38bNA2/2OpSUEVvvb+MbZ3U4i1GFo5i1ahZVWkWG2P+ydHdBpwu4\noNMFXoeXAmA0AAAQLklEQVSRUmyrSCOd8jrRvU13Dh095HUoJknW7FzDldOuZM7qOV6HkhKshZdG\njhw7wgdbPiA5nWGMH7y84mVeWP4CB48eZFRh8u8tmGqshZdGQj1DHK48zIPvPeh1KCZJxvQaA8D8\ndfM5cuyIx9H4nyW8NHJh5wvJzczl9jdutzvsBkSPNj3o2KIjFccqWLB+gdfh+J4lvDSSk5nDkC5D\nAJixqmbHFJOurux9JQAvLH/B40j8zxJemrnujOtoktWE6csT+5Aa4x9X9LkCwFr1DWAnLdLMtf2u\nZfGWxTy06CEOHDlgvS4C4LzTzmPLf2zhlBaneB2K71kLLw2FeoaoOFZhvS4CIisj63iysxMXdbOE\nl4Y27d1EVkYWZ51yltehmCRZs3MNHe/vyDmPnON1KL5mCS8Nnd3hbCqrKnl17ateh2KSpF2zdmzZ\nt4WPt37MjoM7vA7HtyzhpaG++X0pyCtg8tuTWfLlEq/DMUnQqkkrzu5wNooye/Vsr8PxLUt4aUhE\nGNl9JOt3r+eZpc94HY5JkrGnjwWwZV4HS3hp6lt9vwXAtGXTPI7EJMulPS8FrNdFXSzhpakhXYYw\nqsco1u1ex+d7Pvc6HJMEhW0LGdZ1GOPOGEdlVaXX4fiSJbw0lZOZw++G/w6AV1a94nE0JlneuP4N\npo6ZSrPsZl6H4ku+Sngi8hsRWSEiH4vIdBGJ+MBNEVkvIp+IyBIRSb97tsdJ55adKWxTyMa9dgV+\nkKzbvY4/vvdHgvT4hobyVcID5gKnq+qZOA/wubOOaS9S1f5+vn++1w4cPcCanWvIyczxOhSTJFVa\nRf8p/bllzi18uv1Tr8PxHV8lPFV9XVWrDz4sBAq8jCfV5TfPp/i0YmasnGE3BQ2IDMlgaNehgD3r\nIhJfJbwabgBqu42rAq+LyGIRGV9XISIyXkQWicii7du3xz1IvystLOXDLz/khhk3eB2KSZKr+14N\n2OUpkSQ94YnIGyKyNMJrTNg0PwcqgSdrKeZ8VT0bGAX8SEQurK0+VX1YVQeq6sD8/Py4zksqqL5U\nYfrS2fzfO1UeR2OSYWSPkQjCJ9s+sV4XNSQ94anqxap6eoTXywAi8h2gFPi21nLUVVU3u3+3AdMB\n60BYiz1r+sCBdlSwl2HXLaK83OuITKJV97oAePvztz2Oxl98tUsrIiOB24FLVfVgLdM0F5EW1e+B\n4cDS5EWZWhYsEOSv/4CqDI50nUlZmdcRmWR4cPSDvPrtV7ms12Veh+Irvkp4wINAC2Cue8nJFAAR\nOVVEqjsItgfeFpGPgPeAV1TVesnXoqQEmuzvBRvPh6KZlJR4HZFJhkEFgxjRY4TXYfiOr24Aqqo9\navl8MzDaff8Z0C+ZcaWy4mKYNw/+PP8OCvscYNAgBeypZkHwzCfPcM8/7uHB0Q9S0qXE63B8wVcJ\nzyRGcTEUF4/2OgyTZAs3LWTZ9mU8s/QZS3guv+3SmgRa9dUqu5lAgFze+3IApq+Ybr0uXJbwAuRP\n7/+J66dfz4EjB7wOxSTB+Z3Op1lWM7Yd2Ga9LlyW8AIkVGTPugiSrIwshvcYDlivi2qW8AJkcOfB\n5OXmMXPVTK9DMUlydd+raZbdjKwMO1wPdtIiUHIycxjZYySzVs2iSqvIEPt/l+6u6nMVV/e9GhE7\nMw/WwgucUFGInYd2suqrVV6HYpIgMyMTEeHg0YNs2bfF63A8ZwkvYK7ofQU7/nMHvdr18joUkySv\nr32dvMl53DjjRq9D8Zzt0gaM3Qk3eLq26soxPXb8WRdBvj+itfAC6N1N7zL4scFs3GN3Qg6CwraF\ndGzRkYpjFSxYv8DrcDxlCS+A8nLzePvzt5m1apbXoZgkubL3lQC8sPwFjyPxliW8AOrVrhfdW3e3\ny1MCxHpdOCzhBZCIECoKMX/dfOt1ERDndzqfsaeP5XeX/M7rUDxlCS+gSotKrddFgGRlZPHUN5/i\n2n7XBvqaPEt4ATW482BCRSFOyjnJ61BMklRpFQ8vfpi75t3ldSiekSDtzw8cOFAXLbLH2Jpg2n14\nN21+3QaAbT/bRrtm7eJeh4gs9vOjU62FF3A7Du5g+4HgPc0tiFo1acWADgNQlNmrZ9f/hTRkCS/A\n9lXso+P9Hfn9u7/3OhSTJNecfg0Q3Ec4+i7hichEEfnCfabFEhGJeKteERkpIitFZI2I3JHsONNB\ni9wWnNvxXGattuvxgqL6sZ3VvS6CxncJz/WAqvZ3X19re4tIJvBHnOfS9gHGikifZAeZDkqLSlny\n5RLrdREQ1b0uFGXNzjVeh5N0fk149TkHWKOqn6nqEeAZYEw93zERhIpCANbrIkDmXjeXXbfvok9+\n8NoIfk14t4jIxyIyVURaRxjfEQhvkmxyP/saERkvIotEZNH27XZwvibrdRE8vfN70yy7GaoauF4X\nniQ8EXlDRJZGeI0BHgK6A/2BLUCkS8MjXTkZccmp6sOqOlBVB+bn58dtHtKFiDB1zFQeHP2g16GY\nJPrBKz+g9a9bB+5ZF57cHkpVL27IdCLyCBBpX2sTcFrYcAGwOQ6hBdKFnS/0OgSTZF/u/5I9FXuY\nsXIGfU/u63U4SeO7XVoR6RA2eDmwNMJk7wOFItJVRHKAa4AZyYgvXT3/6fP8edGfvQ7DJMnVfa8G\ngnd5iu8SHnCfiHwiIh8DFwH/DiAip4rIbABVrQRuAV4DlgPTVHWZVwGng+c+fY4JZROo0iqvQzFJ\nMLLHSAThk22fsOPgDq/DSRrfJTxVvU5Vz1DVM1X1UlXd4n6+WVVHh003W1WLVLW7qv63dxGnh1BR\niK0HtrJos3W9C4Kg9rrwXcIz3hjVYxQZksHMlXa2NihuOOsGurfuTpdWXbwOJWns5gHmuP7/eyHb\n9uzlhUuWUFzsdTQmFdnNA0xKKC+HZdNL2bK1kqEjD1Be7nVEJhlUlbL1ZazdudbrUJLCEp4BoKwM\nqv7vVvjTUo4eaE5ZmdcRmWS49x/3ctHfLuK37/zW61CSwhKeAaCkBHKzM8nMhJwcZ9ikvws6XQDA\niyteDESvC3surQGguBjmzXNaeiUl2DG8gDi/0/k0y2rGtgPb+HT7p2l/EbK18MxxxcVw552W7IIk\nKyOL4T2GA/Dyypc9jibxLOEZE3DVvS6e//R5jyNJPNulNSbgRvYYyX8P/W+uPfNar0NJOEt4xgRc\nqyatuGtwMJ5kZru0xhh2HtrJ1c9dzU9f/anXoSSUJTxjDF8d/Ippn07jz4v+nNbPurCEZ4w5/qyL\nimMVLFi/wOtwEsYSnjEGgCt7XwnAi8tf9DiSxLGEZ4wB4PLelwPp3evCEp4xBnB6XTTPbk5ebh77\njuzzOpyEsMtSjDGA0+ti621baZ7T3OtQEsZaeMaY46qT3a5DuzyOJDF81cITkWeBnu5gK2C3qvaP\nMN16YB9wDKj08w0HjUklVVpF/yn9WbptKVtv20p+8/R6tKmvWniqerWq9neT3AtAXaeLLnKntWRn\nTJxkSAY5mTkoypw1c7wOJ+58lfCqiYgA3wKe9joWY4LmmtOvAdLzEY6+THjAYGCrqq6uZbwCr4vI\nYhEZn8S4jEl7Y3qOAWDe2jLuvvdIWt3uP+kJT0TeEJGlEV5jwiYbS92tu/NV9WxgFPAjEbmwjvrG\ni8giEVm0ffv2OM2FMemrsG0hrfZcwBE9xITHFjBsGGmT9JJ+0kJVL65rvIhkAVcAA+ooY7P7d5uI\nTAfOAd6qZdqHgYfBeWpZlGEbExjl5bDnpbuh3TJ0y5lUVDh3wk6HG8P66iyt62JghapuijRSRJoD\nGaq6z30/HJiUzACNSWdlZSAbStB1JQBkZqfPM078eAzvGmrszorIqSJS/Xj09sDbIvIR8B7wiqq+\nmuQYjUlbJSWQmwsZGZCVBQ8+mB6tO/BhC09Vvxvhs83AaPf9Z0C/JIdlTGCk8wOdfJfwjDHeKy5O\nr0RXzY+7tMYYkxCW8IwxgWEJzxgTGJbwjDGBYQnPGBMYlvCMMYFhCc8YExiW8IwxgWEJzxgTGJbw\njDGBYQnPGBMYlvCMMYFhCc8YExiW8IwxgWEJzxgTGJbwjDGBYQnPGBMYlvCMMYHhScITkatEZJmI\nVInIwBrj7hSRNSKyUkRG1PL9riLyroisFpFnRSQnOZEbY1KZVy28pTjPnj3hWbIi0gfnqWV9gZHA\nn0QkM8L3fw08oKqFwC7gxsSGa4xJB54kPFVdrqorI4waAzyjqhWqug5Yg/OQ7eNERIChwPPuR38D\nLktkvMaY9OC3p5Z1BBaGDW9yPwvXFtitqpV1THOciIwHxruDFSKyNE6x+kk7YIfXQSRAus4XpO+8\n9fQ6gLokLOGJyBvAKRFG/VxVX67taxE+0yim+ecI1YeBh92YFqnqwNqmTVU2X6knXedNRBZ5HUNd\nEpbwVPXiKL62CTgtbLgA2Fxjmh1AKxHJclt5kaYxxpiv8dtlKTOAa0QkV0S6AoXAe+ETqKoCbwJX\nuh99B6itxWiMMcd5dVnK5SKyCSgGXhGR1wBUdRkwDfgUeBX4kaoec78zW0ROdYu4HfgPEVmDc0zv\nLw2s+uE4zoaf2HylnnSdN1/PlzgNJmOMSX9+26U1xpiEsYRnjAmMtE94sXZjSxUiMlFEvhCRJe5r\ntNcxxUJERrrLZY2I3OF1PPEiIutF5BN3Gfn6Eo76iMhUEdkWfm2riLQRkblut8+5ItLayxhrSvuE\nR+zd2FLJA6ra333N9jqYaLnL4Y/AKKAPMNZdXuniIncZpfp1eH/F2XbC3QHMc7t9znOHfSPtE14s\n3diMZ84B1qjqZ6p6BHgGZ3kZH1HVt4CdNT4eg9PdE3zY7TPtE14dOgIbw4br7KKWIm4RkY/dXQ1f\n7Uo0Ujoum2oKvC4ii91uj+mmvapuAXD/nuxxPCfwW1/aqCSwG5uv1DWfwEPA3TjzcDfwO+CG5EUX\nVym3bBrhfFXdLCInA3NFZIXbUjJJkBYJL4Hd2HylofMpIo8AsxIcTiKl3LJpKFXd7P7dJiLTcXbf\n0ynhbRWRDqq6RUQ6ANu8DihckHdp6+3Glkrclava5Tgna1LV+0Che6PXHJyTSzM8jilmItJcRFpU\nvweGk9rLKZIZON09wYfdPtOihVcXEbkc+AOQj9ONbYmqjlDVZSJS3Y2tkrBubCnqPhHpj7Prtx64\nydtwoqeqlSJyC/AakAlMdbsdprr2wHTnlo5kAU+p6qvehhQ9EXkaKAHauV1FJwC/AqaJyI3A58BV\n3kX4dda1zBgTGEHepTXGBIwlPGNMYFjCM8YEhiU8Y0xgWMIzxgSGJTxjTGBYwjPGBIYlPOM7IvIv\n7k0Qmri9E5aJyOlex2VSn114bHxJRO4BmgBNgU2qOtnjkEwasIRnfMntQ/s+cBg4L8W7/RmfsF1a\n41dtgJOAFjgtPWNiZi0840siMgPnTsddgQ6qeovHIZk0kPZ3SzGpR0SuBypV9Sn3+RbviMhQVZ3v\ndWwmtVkLzxgTGHYMzxgTGJbwjDGBYQnPGBMYlvCMMYFhCc8YExiW8IwxgWEJzxgTGP8fxNnUZFT+\n+3QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x178e9131eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = 0.3\n",
    "\n",
    "x = ((np.random.rand(10)*2)-1)*10\n",
    "y = ((np.random.rand(10)*2)-1)*10\n",
    "pair = zip(x, y)\n",
    "\n",
    "# The projection matrix\n",
    "P = np.array([[1, a],\n",
    "              [0, 0]])\n",
    "\n",
    "for item in pair:\n",
    "    x, y = item\n",
    "\n",
    "    # Form the vector z in R2 from the point (x, y)\n",
    "    z = np.array([[x],\n",
    "                  [y]])\n",
    "\n",
    "    # Project z onto the line y = 0\n",
    "    Pz = P @ z\n",
    "\n",
    "    # Plot the result\n",
    "    plt.plot(x, y, 'b.')\n",
    "    plt.plot(Pz[0], Pz[1], 'rx')\n",
    "    plt.plot([x,Pz[0]],[y,Pz[1]],color='green', linestyle='dashed')\n",
    "\n",
    "plt.plot([-10,10],[0,0])\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.axis('square')\n",
    "plt.ylim(-10, 10)\n",
    "plt.xlim(-10, 10)\n",
    "plt.title(\"Projection of 10 Random Coordinates onto y=0\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An Engineering Application\n",
    "\n",
    "### Example 1: Nearest point projection\n",
    "\n",
    "If we are given a point $x \\in S$, suppose that we would like to approximate $x$ with a point in $V \\subset S$. Assuming that $x \\notin V$, we need to find the point in $V$ that is closest to $x$; that is, $\\hat{x}$. This point is given by the orthogonal projection of $x$ onto $V$. \n",
    "\n",
    "In this example, let $V$ be the plane in $\\mathbb{R}^3$ spanned by the vectors \n",
    "\n",
    "$$\n",
    "y = \n",
    "    \\left[ {\\begin{array}{c}\n",
    "    1  \\\\\n",
    "    0  \\\\\n",
    "    1  \\\\\n",
    "    \\end{array} } \\right]\n",
    "\\textrm{and }\n",
    "z =\n",
    "    \\left[ {\\begin{array}{c}\n",
    "    0 \\\\\n",
    "    1 \\\\\n",
    "    1 \\\\\n",
    "    \\end{array} } \\right]\n",
    "$$.\n",
    "\n",
    "If \n",
    "\n",
    "$$\n",
    "A=\n",
    "  \\left[ {\\begin{array}{cc}\n",
    "   1 & 0 \\\\\n",
    "   0 & 1 \\\\\n",
    "   1 & 1 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$\n",
    "\n",
    "then the projection matrix that projects points in $\\mathbb{R}^3$ orthogonally onto the plane spanned by $y$ and $z$ is described by $P = A(A^HA)^{-1}A^H$. Note that the columns of $A$ are the vectors that span $V$. In the python code below, we project the vector \n",
    "\n",
    "$$\n",
    "  q =\\left[ {\\begin{array}{c}\n",
    "   7 \\\\\n",
    "    8 \\\\\n",
    "    9 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$\n",
    "\n",
    "onto the plane $V = span(y,z)$, and get \n",
    "\n",
    "$$\n",
    "  \\hat{x} = Pq =\\left[ {\\begin{array}{c}\n",
    "   5 \\\\\n",
    "    6 \\\\\n",
    "    11 \\\\\n",
    "  \\end{array} } \\right]\n",
    "$$.\n",
    "\n",
    "This vector is the closest approximation of $q$ on $V$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection matrix:\n",
      "[[ 0.66666667 -0.33333333  0.33333333]\n",
      " [-0.33333333  0.66666667  0.33333333]\n",
      " [ 0.33333333  0.33333333  0.66666667]] \n",
      "\n",
      "Closest approximation in V:\n",
      "[[  5.]\n",
      " [  6.]\n",
      " [ 11.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 0],\n",
    "              [0, 1],\n",
    "              [1, 1]])\n",
    "\n",
    "P = A @ np.linalg.inv(A.conj().T @ A) @ A.conj().T\n",
    "print('Projection matrix:')\n",
    "print(P, \"\\n\")\n",
    "\n",
    "q = np.array([[7],\n",
    "              [8],\n",
    "              [9]])\n",
    "Pq = P @ q\n",
    "print('Closest approximation in V:')\n",
    "print(Pq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: Representing 3D space in 2D\n",
    "\n",
    "Suppose we want to represent a 3D coordinate axis on a screen. This is the most general case of rendering any 3D object on a computer screen; if we can represent the coordinate axis, then we can represent any object which we might place within that Euclidean space. Consider 3D plots generated by a program such as MATLAB or Matplotlib; what we see as a good 3D representation of a line or surface is simply a projection onto a 2D surface: your computer screen. So instead of relying on the plotting library to figure that out for us--what if we tried to calculate where the projection would fall ourselves?\n",
    "\n",
    "Suppose we have $x = y = z = \\{0, 1, 2, 3, 4, 5\\}$ and we want to plot $(x, 0, 0), (0, y, 0), (0, 0, z)$ for all values of $x, y, z$ in their respective sets. This is, of course, simply the first five points all the coordinate axis in each direction. Suppose we project onto an observation plane that is 10 units from the origin. The plane will be defined by a vector $\\vec{r} = (r, \\phi, \\theta)$ from the origin, where $r$ is the magnitude of the vector, $\\phi$ is the angle from $x$ towards $y$ in the $x-y$ plane, and $\\theta$ is the angle measured from $z$, which points directly upwards.\n",
    "\n",
    "<img src=\"files/3d_coordinate_axis.jpg\" width=\"300\" />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Manipulation parameters (angles in degrees)\n",
    "theta = 45 # Between {0, 90deg}\n",
    "phi = 45 # Between {0, 90deg}\n",
    "r = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two ways we could approach this problem.\n",
    "\n",
    "1. Draw a line from the starting point to the plane in the direction of the plane's normal vector and see where it intersects.\n",
    "\n",
    "2. Find the projection to a space *orthogonal* to the plane of interest, $W$, and then find $I-W$ which is the projection onto the plane (since the range and null space of a projection are disjoint).\n",
    "\n",
    "Let's demonstrate both methods with code examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Approach 1\n",
    "\n",
    "First, we'll see where a vector from the point of interest in the direction of the plane's normal vector intersects the plane.\n",
    "\n",
    "Our normal unit vector can be defined as $\\hat{n} \n",
    "= \\left\\langle \\frac{r \\sin \\theta \\cos \\phi}{\\| \\vec{r} \\|}, \\frac{r \\sin \\theta \\sin \\phi}{\\| \\vec{r} \\|}, \\frac{r \\cos \\theta}{\\| \\vec{r} \\|} \\right\\rangle \n",
    "= \\langle \\sin \\theta \\cos \\phi, \\sin \\theta \\sin \\phi, \\cos \\theta \\rangle$.\n",
    "\n",
    "We can write a parametric equation for a line from any point $(x, y, z) \\in \\mathbb{R}^3$ in the direction of $\\hat{n}$ parameterized by $t$. The form of that equation is\n",
    "\n",
    "$$\n",
    "X(t) = \\left( x + t \\sin \\theta \\cos \\phi, y + t \\sin \\theta \\sin \\phi, z + t \\cos \\theta \\right), t \\in \\mathbb{R}\n",
    "$$\n",
    "\n",
    "We want to find the value $t_0$ such that $X(t_0)$ satisfies the plane equation. \n",
    "\n",
    "A plane defined by a point $p = (x_1, y_1, z_1)$ and a normal vector $\\hat{n} = \\langle A, B, C \\rangle$ has the form\n",
    "\n",
    "$$\n",
    "A(x - x_1) + B(y - y_1) + C(z - z_1) = 0\n",
    "$$\n",
    "\n",
    "In this case, $p = X(t_0)$. Substituting appropriately yields, as the projection of a point onto the plane,\n",
    "\n",
    "<!--\n",
    "\\frac{r \\sin \\theta \\cos \\phi}{\\| r \\|} \\left[ x - \\left(x_0 + \\frac{r \\sin \\theta \\cos \\phi}{\\| r \\|}t_0 \\right) \\right] + \\frac{r \\sin \\theta \\sin \\phi}{\\| r \\|} \\left[ y - \\left(y_0 + \\frac{r \\sin \\theta \\sin \\phi}{\\| r \\|}t_0 \\right) \\right] + \\frac{r \\cos \\theta}{\\| r \\|} \\left[ z - \\left(z_0 + \\frac{r \\cos \\theta}{\\| r \\|}t_0 \\right) \\right] = 0\n",
    "-->\n",
    "\n",
    "$$\n",
    "\\sin \\theta \\cos \\phi \\left(x + t_0 \\sin \\theta \\cos \\phi \\right) + \\sin \\theta \\sin \\phi \\left(y + t_0 \\sin \\theta \\sin \\phi \\right) + \\cos \\theta \\left(z + t_0 \\cos \\theta \\right) = 0\n",
    "$$\n",
    "\n",
    "We define the following relations\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "A &= \\sin \\theta \\cos \\phi \\\\\n",
    "B &= \\sin \\theta \\sin \\phi \\\\\n",
    "C &= \\cos \\theta\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Solving for $t_0$ then yields\n",
    "\n",
    "$$\n",
    "t_0 = - \\frac{A x + B y + C z}{A^2 + B^2 + C^2}\n",
    "$$\n",
    "\n",
    "We can then formulate an equation for the projection of any starting point in the direction of the normal vector to some point in the plane as\n",
    "\n",
    "$$\n",
    "P(x, y, z) = (x + A t_0, y + B t_0, z + C t_0)\n",
    "$$\n",
    "\n",
    "The projection matrix could then be the stacked values of the vectors\n",
    "\n",
    "$$\n",
    "P(1, 0, 0) \\\\\n",
    "P(0, 1, 0) \\\\\n",
    "P(0, 0, 1) \\\\\n",
    "$$\n",
    "\n",
    "To actually render the result on a screen, we'd have to find the rotation matrix that brings the plane we projected onto back into our $(x, y)$ coordinate system and then select out just the $x$ and $y$ components of that $3 \\times 3$ matrix.\n",
    "\n",
    "*[Contributing Source](https://math.stackexchange.com/a/1140391)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection Matrix:\n",
      "[[ 0.75       -0.25       -0.35355339]\n",
      " [-0.25        0.75       -0.35355339]\n",
      " [-0.35355339 -0.35355339  0.5       ]]\n",
      "x-axis projection:\n",
      "[[ 0.          0.75        1.5         2.25        3.          3.75      ]\n",
      " [ 0.         -0.25       -0.5        -0.75       -1.         -1.25      ]\n",
      " [ 0.         -0.35355339 -0.70710678 -1.06066017 -1.41421356 -1.76776695]]\n",
      "y-axis projection:\n",
      "[[ 0.         -0.25       -0.5        -0.75       -1.         -1.25      ]\n",
      " [ 0.          0.75        1.5         2.25        3.          3.75      ]\n",
      " [ 0.         -0.35355339 -0.70710678 -1.06066017 -1.41421356 -1.76776695]]\n",
      "z-axis projection:\n",
      "[[ 0.         -0.35355339 -0.70710678 -1.06066017 -1.41421356 -1.76776695]\n",
      " [ 0.         -0.35355339 -0.70710678 -1.06066017 -1.41421356 -1.76776695]\n",
      " [ 0.          0.5         1.          1.5         2.          2.5       ]]\n"
     ]
    }
   ],
   "source": [
    "coords = np.array([1, 2, 3, 4, 5])\n",
    "\n",
    "x_axis = np.zeros((3, 1))\n",
    "y_axis = np.zeros((3, 1))\n",
    "z_axis = np.zeros((3, 1))\n",
    "for val in coords:\n",
    "    x_axis = np.concatenate((x_axis, np.array([[val, 0, 0]]).T), axis=1)    \n",
    "    y_axis = np.concatenate((y_axis, np.array([[0, val, 0]]).T), axis=1)\n",
    "    z_axis = np.concatenate((z_axis, np.array([[0, 0, val]]).T), axis=1)\n",
    "\n",
    "theta_rad = np.deg2rad(theta)\n",
    "phi_rad = np.deg2rad(phi)\n",
    "\n",
    "A = np.sin(theta_rad) * np.cos(phi_rad)\n",
    "B = np.sin(theta_rad) * np.sin(phi_rad)\n",
    "C = np.cos(theta_rad)\n",
    "\n",
    "def t0(x, y, z):\n",
    "    t0 = (-1) * (A*x + B*y + C*z) / (A**2 + B**2 + C**2)\n",
    "    return t0\n",
    "\n",
    "def P_func(x, y, z):\n",
    "    t1 = t0(x, y, z)\n",
    "    x1 = x + A * t1\n",
    "    y1 = y + B * t1\n",
    "    z1 = z + C * t1\n",
    "    return np.array([x1, y1, z1])\n",
    "\n",
    "# Create the projection matrix\n",
    "P_matrix = np.zeros((3, 3))\n",
    "P_matrix[0, :] = P_func(1, 0, 0)\n",
    "P_matrix[1, :] = P_func(0, 1, 0)\n",
    "P_matrix[2, :] = P_func(0, 0, 1)\n",
    "print(\"Projection Matrix:\")\n",
    "print(P_matrix)\n",
    "\n",
    "proj_x = P_matrix @ x_axis\n",
    "proj_y = P_matrix @ y_axis\n",
    "proj_z = P_matrix @ z_axis\n",
    "\n",
    "print(\"x-axis projection:\")\n",
    "print(proj_x)\n",
    "print(\"y-axis projection:\")\n",
    "print(proj_y)\n",
    "print(\"z-axis projection:\")\n",
    "print(proj_z)\n",
    "\n",
    "# Bundle the results for comparison with Approach 2 later\n",
    "approach1 = [P_matrix.copy(), proj_x.copy(), proj_y.copy(), proj_z.copy()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Approach 2\n",
    "\n",
    "The space orthogonal to the space of interest ($W$ being the plane we're projecting onto) can be defined as $W^{\\perp} = \\textrm{span}(\\hat{n})$; anything parallel to the normal vector obviously has no projection onto the plane defined by said normal vector.\n",
    "\n",
    "Using the same values defined in Approach 1, the equation of our plane is\n",
    "\n",
    "$$\n",
    "A (x + A t_0) + B (y + B t_0) + C (z + C t_0) = 0\n",
    "$$\n",
    "\n",
    "Again, our normal vector is \n",
    "\n",
    "$$\n",
    "\\hat{n} = \\langle A, B, C \\rangle\n",
    "$$\n",
    "\n",
    "The definition of a projection of some vector $\\vec{a}$ onto our normal $\\hat{n}$ is\n",
    "\n",
    "$$\n",
    "proj_\\hat{n} \\vec{a} = \\frac{\\hat{n} \\cdot \\vec{a}}{\\| \\hat{n} \\|^2} \\hat{n}\n",
    "$$\n",
    "\n",
    "which is equivalent to\n",
    "\n",
    "$$\n",
    "proj_\\hat{n} \\vec{a} = \\left( \\frac{\\hat{n}}{\\| \\hat{n} \\|^2} \\hat{n} \\right) \\cdot \\vec{a}\n",
    "$$\n",
    "\n",
    "We can therefore set up our projection matrix (onto $W^\\perp$) as\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "    P_{W^\\perp} &= \\frac{1}{A^2 + B^2 + C^2} \n",
    "    \\left[ {\\begin{array}{c}\n",
    "    A  \\\\\n",
    "    B  \\\\\n",
    "    C  \\\\\n",
    "    \\end{array} } \\right]\n",
    "    \\left[ {\\begin{array}{c}\n",
    "    A & B & C \\\\\n",
    "    \\end{array} } \\right] \\\\\n",
    "    &= \\frac{1}{A^2 + B^2 + C^2} \n",
    "    \\left[ {\\begin{array}{c}\n",
    "    AA & AB & AC \\\\\n",
    "    BA & BB & BC \\\\\n",
    "    CA & CB & CC \\\\\n",
    "    \\end{array} } \\right]\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Remembering that $W$ and $W^\\perp$ are distjoint, that means that the projection onto $W$ is simply $P_W = I - P_{W^\\perp}$.\n",
    "\n",
    "See code implementation for full example. The results are compared with Approach 1 in the code block following.\n",
    "\n",
    "*[Contributing Source](https://math.stackexchange.com/a/1140397)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.75       -0.25       -0.35355339]\n",
      " [-0.25        0.75       -0.35355339]\n",
      " [-0.35355339 -0.35355339  0.5       ]]\n"
     ]
    }
   ],
   "source": [
    "coords = np.array([1, 2, 3, 4, 5])\n",
    "\n",
    "x_axis = np.zeros((3, 1))\n",
    "y_axis = np.zeros((3, 1))\n",
    "z_axis = np.zeros((3, 1))\n",
    "for val in coords:\n",
    "    x_axis = np.concatenate((x_axis, np.array([[val, 0, 0]]).T), axis=1)    \n",
    "    y_axis = np.concatenate((y_axis, np.array([[0, val, 0]]).T), axis=1)\n",
    "    z_axis = np.concatenate((z_axis, np.array([[0, 0, val]]).T), axis=1)\n",
    "\n",
    "theta_rad = np.deg2rad(theta)\n",
    "phi_rad = np.deg2rad(phi)\n",
    "\n",
    "A = np.sin(theta_rad) * np.cos(phi_rad)\n",
    "B = np.sin(theta_rad) * np.sin(phi_rad)\n",
    "C = np.cos(theta_rad)\n",
    "\n",
    "P_perp = (1 / (A**2 + B**2 + C**2)) * np.array([[A], [B], [C]]) @ np.array([[A, B, C]])\n",
    "P_matrix = np.identity(3) - P_perp\n",
    "print(P_matrix)\n",
    "\n",
    "proj_x = P_matrix @ x_axis\n",
    "proj_y = P_matrix @ y_axis\n",
    "proj_z = P_matrix @ z_axis\n",
    "\n",
    "# Bundle results for comparison with Approach 1 later.\n",
    "approach2 = [P_matrix.copy(), proj_x.copy(), proj_y.copy(), proj_z.copy()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test the results of Approach 1 and Approach 2 by comparing the projection matrices and subsequent projections of each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection matrices and projections were equal.\n"
     ]
    }
   ],
   "source": [
    "for a1, a2 in zip(approach1, approach2):\n",
    "    if not np.array_equal(a1, a2):\n",
    "        raise Exception(\"Projection matrices or projections do not match.\")\n",
    "print(\"Projection matrices and projections were equal.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Couldn't figure out the rotation matrices to view the axis \n",
    "# projected onto the plane on a plot, but left my attempt here \n",
    "# in case someone in the future would like to try.\n",
    "\n",
    "# def Rz(alpha):\n",
    "#     R = np.array([[np.cos(alpha), -np.sin(alpha), 0],\n",
    "#                   [np.sin(alpha), np.cos(alpha), 0],\n",
    "#                   [0, 0, 1]])\n",
    "#     return R\n",
    "\n",
    "# def Ry(beta):\n",
    "#     R = np.array([[np.cos(beta), 0, np.sin(beta)],\n",
    "#                   [0, 1, 0],\n",
    "#                   [-np.sin(beta), 0, np.cos(beta)]])\n",
    "#     return R\n",
    "\n",
    "# def Rx(gamma):\n",
    "#     R = np.array([[1, 0, 0],\n",
    "#                   [0, np.cos(gamma), -np.sin(gamma)],\n",
    "#                   [0, np.sin(gamma), np.cos(gamma)]])\n",
    "#     return R\n",
    "\n",
    "# R = Rz(-phi_rad) @ Ry(-(np.pi/2 - theta_rad)) @ Rx(0)\n",
    "\n",
    "# # Rotate our x-y viewing plane\n",
    "# P_flatten = np.array([[1, 0, 0],\n",
    "#                       [0, 1, 0]])\n",
    "# P_viewing = P_flatten @ R\n",
    "\n",
    "# proj_x_2d = P_viewing @ proj_x\n",
    "# proj_y_2d = P_viewing @ proj_y\n",
    "# proj_z_2d = P_viewing @ proj_z\n",
    "\n",
    "# print(proj_x_2d)\n",
    "# print(proj_y_2d)\n",
    "# print(proj_z_2d)\n",
    "\n",
    "# plt.plot(proj_x_2d[0,:], proj_x_2d[1,:], label='x-axis')\n",
    "# plt.plot(proj_y_2d[0,:], proj_y_2d[1,:], label='y-axis')\n",
    "# plt.plot(proj_z_2d[0,:], proj_z_2d[1,:], label='z-axis')\n",
    "# plt.legend()\n",
    "# plt.show()"
   ]
  }
 ],
 "metadata": {
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}