psets/pset1sol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "08cd437b",
   "metadata": {},
   "source": [
    "# 18.065 Pset 1 Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d53b289",
   "metadata": {},
   "source": [
    "## Problem 1 (4+4+4+4+4 points)\n",
    "\n",
    "Recall from class that multiplying an $m \\times p$ by a $p \\times n$ matrix costs $mnp$ scalar multiplications (and a similar number of additions) by the standard (practical) algorithms.\n",
    "\n",
    "Matrix multiplication is **not commutative** ($AB \\ne BA$ in general), but it **is associative**: $(AB)C=A(BC)$.   It turns out that where you put the parentheses (i.e. in *what order* you do the multiplications) can make a *huge* difference in computational cost.\n",
    "\n",
    "**(a)** If $x \\in \\mathbb{R}^n$ and $A,B$ are $n \\times n$ matrices, compare the scalar multiplication counts of $(AB)x$ vs. $A(Bx)$, i.e. if we do the multiplications in the order indicated by the parentheses.\n",
    "\n",
    "**(b)** If $x, b\\in \\mathbb{R}^n$, **how many scalar multiplications** does the computation $$p = (I - (xx^T)/(x^T x)) b$$ take if we *do it in the order indicated by the parentheses*?  (Note that dividing by a scalar $\\alpha$ is equivalent to multiplying by $\\alpha^{-1}$ at the negligible cost of one scalar division.)\n",
    "\n",
    "**(c)** Explain how to compute the *same* $p$ as in part (b) using as *few multiplications as possible*.  Outline the sequence of computational steps, and give the count of multiplications.\n",
    "\n",
    "**(d)** $p^T x = $ what?\n",
    "\n",
    "**(e)** Implement your algorithm from (c) in Julia, filling in the code below, and time it for $n=1000$ using the `@btime` macro from the [BenchmarkTools package](https://github.com/JuliaCI/BenchmarkTools.jl), along with the algorithm from part (b), following the outline below.  How does the ratio of the two times compare to your ratio of multiplication counts?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5d43b59",
   "metadata": {},
   "source": [
    "### Solution:\n",
    "\n",
    "**(a)** $AB$ requires $n^3$ scalar multiplications ($n^2$ row×col dot products), and $Bx$ requires $n^2$ multiplications ($n$ dot products).  So, $(AB)x$ requires $\\boxed{n^3 + n^2}$ multiplications, while $A(Bx)$ requires $\\boxed{2n^2}$ multiplications (a *lot* fewer for large $n$!).\n",
    "\n",
    "**(b)** $$p = (I - (xx^T)/(x^T x)) b$$, computed in this order, requires (1) $n$ multiplies for $x^T x$, (2) $n^2$ multiplies for $xx^T$, (3) $n^2$ multiplies to scale $xx^T$ by $(x^T x)^{-1}$ (+ 1 division), (4) no multiplies to subtract from $I$, and (5) $n^2$ multiplies to multiply by $b$.  Total: $\\boxed{3n^2 + n}$ scalar multiplications.\n",
    "\n",
    "**(c)** If we re-arrange it as $$p = b - x((x^T b)/(x^T x))$$, then it requires only $\\boxed{3n+1}$ multiplications: $2n$ for the two dot products, 1 for the scalar division (counting it as a scalar multiplication by the inverse as instructed), and $n$ to scale $x$.\n",
    "\n",
    "**(d)** By remembering your linear algebra, you might recognize $P = (xx^T)/(x^T x)$ as **orthogonal projection** onto the direction of $x$, and hence $I-P$ is projection **perpendicular** to $x$.  So, we must have $\\boxed{p^T x = 0}$.  But we can also compute this explicitly:\n",
    "$$\n",
    "p^T x = x^T p = x^T \\left[ b - x((x^T b)/(x^T x)) \\right] = \\\\\n",
    "x^T b - (x^T x)(x^T b)/(x^T x) = x^T b - x^T b = 0 \\, .\n",
    "$$\n",
    "\n",
    "**(e)** See the following code.  On my laptop, my part (c) solution is about **3400× faster** than part (b), which is actually bigger than the ratio of $(3n^2 + n) / (3n+1) = n = 1000$ you might predict just by counting multiplications, but it is the right order of magnitude.\n",
    "\n",
    "Of course, computers are *much* more complicated than simple arithmetic counts than capture, and you will almost never be able to predict performance ratios simply by arithmetic counts.   To get a sense of this complexity, it's interesting to plot $\\frac{1}{n}\\text{(time for b)}/\\text{(time for c)}$ versus $n$.  Arithmetic counts would predict this ratio to be a constant, but in fact it is much more complicated. (You were *not* required to do this — I show it below only for interest.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b149871c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hooray, part (c) and part (b) agree!\n",
      "\n",
      "part (b): \n",
      "  4.143 ms (7 allocations: 22.90 MiB)\n",
      "\n",
      "part (c): \n",
      "  1.279 μs (2 allocations: 15.88 KiB)\n"
     ]
    }
   ],
   "source": [
    "using LinearAlgebra, BenchmarkTools\n",
    "\n",
    "# algorithm from part (b)\n",
    "function part_b(x, b)\n",
    "    return (I - (x*x')*(x'*x)^-1) * b\n",
    "end\n",
    "\n",
    "# algorithm from part (c)\n",
    "function part_c(x, b)\n",
    "    # CHANGE THIS:\n",
    "    return b - x * ((x'*b)/(x'*x))\n",
    "end\n",
    "\n",
    "# test and benchmark on random vectors:\n",
    "n = 1000\n",
    "x, b = rand(n), rand(n)\n",
    "\n",
    "# test it first — should give same answer up to roundoff error\n",
    "if part_c(x, b) ≈ part_b(x, b)\n",
    "    println(\"Hooray, part (c) and part (b) agree!\")\n",
    "else\n",
    "    error(\"You made a mistake: part (c) and part (b) do not agree!\")\n",
    "end\n",
    "\n",
    "# benchmark it:\n",
    "\n",
    "println(\"\\npart (b): \")\n",
    "@btime part_b($x, $b);\n",
    "\n",
    "println(\"\\npart (c): \")\n",
    "@btime part_c($x, $b);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "27552024",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3450.351837372948"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4.413e-3 / 1.279e-6 # ratio of times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bf2cee16",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 10\n",
      "n = 14\n",
      "n = 21\n",
      "n = 30\n",
      "n = 43\n",
      "n = 62\n",
      "n = 89\n",
      "n = 127\n",
      "n = 183\n",
      "n = 264\n",
      "n = 379\n",
      "n = 546\n",
      "n = 785\n",
      "n = 1129\n",
      "n = 1624\n",
      "n = 2336\n",
      "n = 3360\n",
      "n = 4833\n",
      "n = 6952\n",
      "n = 10000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "20-element Vector{Float64}:\n",
       " 0.4198663315479747\n",
       " 0.4024561123451314\n",
       " 0.5016266576985883\n",
       " 0.5217398342398343\n",
       " 0.5499016471424656\n",
       " 0.5297567217200534\n",
       " 0.7579014598586505\n",
       " 1.0874537900252532\n",
       " 1.0783124610718169\n",
       " 0.7767890148134123\n",
       " 0.9204295359856696\n",
       " 1.2168404617858037\n",
       " 2.2617932818717397\n",
       " 3.025106985400225\n",
       " 3.060177735357155\n",
       " 3.1275025975637285\n",
       " 3.0167103840430975\n",
       " 7.520125153516695\n",
       " 8.25342810247817\n",
       " 8.673908026056296"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# for interest, let's plot (time for b) / (n*(time for c)) versus n,\n",
    "# and see how different it is from a constant:\n",
    "\n",
    "ns = round.(Int, 10 .^ range(1, 4, length=20)) # log-spaced points\n",
    "\n",
    "t_bc = map(ns) do n\n",
    "    @show n\n",
    "    x, b = rand(n), rand(n)\n",
    "    @belapsed(part_b($x, $b)), @belapsed(part_c($x, $b))\n",
    "end\n",
    "\n",
    "t_ratio = first.(t_bc) ./ (ns .* last.(t_bc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d6fa1641",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHJCAYAAACMppPqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbO0lEQVR4nO3deVxUVf8H8M+wLyKCgqyCqYiouJK5oJL7gibaY2nmUplmubboY2VqT2a2uISV9ZTZo6Ypmqllbii55JLSopYLKiqmogGKIgzn98f5DTAMywzMzJ3l83695sVw58693xkuzJdzvucclRBCgIiIiMgGOSgdABEREZGpMNEhIiIim8VEh4iIiGwWEx0iIiKyWUx0iIiIyGYx0SEiIiKbxUSHiIiIbBYTHSIiIrJZTHSIiIjIZjHRMbGcnBy8/PLL6NmzJ/z8/KBSqfDGG2+Uu39+fj4iIyPx9ttvG3yuv/76Cy4uLvjll1+qEXGx7777DvHx8ahbty5cXFzg6+uLbt26YeXKlcjPzzfKOZS2devWCn8eROZW+m9EcnIyVCoVkpOTTXK+K1eu4I033sDx48eNfuzly5dDpVLh/PnzFe5X0e+hSqXC888/b/TYbMWqVauwcOFCne3nz5+HSqXCu+++W+kx3njjDahUqiqdX/MzPnLkSJWebw5MdEwsMzMTy5YtQ15eHh555JFK91+6dClu3bqFF154weBzRUREYPjw4ZgyZUoVIi0mhMDo0aMxYMAAFBYW4v3338eOHTvw5ZdfokWLFnjuueewdOnSap3DUmzduhWzZ89WOgyicrVu3RoHDhxA69atTXL8K1euYPbs2SZJdPTF38OqKy/RMcTTTz+NAwcOGCcgC+SkdAC2LiwsDLdu3YJKpcKNGzfw2WeflbtvQUEBFixYgDFjxsDT07NK53v++efRtm1b7N+/Hx06dKjSMRYsWIDly5dj9uzZeP3117Uei4+Px8svv4wzZ85U6dj2Ijc3Fx4eHnZ3bmPIz8+HSqWCkxP/PAFAzZo18dBDDykdBtmwkJAQhISEKB2GybBFx8RUKpXeTYKbNm3C5cuXMWLECJ3HTp8+jWHDhsHf3x+urq5o0qQJEhMTdfZr06YNmjRpgo8//rhK8ebn52P+/PmIjIzEa6+9VuY+AQEB6NSpU9H3N2/exHPPPYfg4GC4uLjggQcewMyZM5GXl6f1PE0T9BdffIHGjRvD3d0dbdu2xcGDByGEwIIFC1C/fn3UqFEDDz/8sE4y1bVrVzRr1gwpKSl46KGH4O7ujuDgYLz22mtQq9VF+5XX1K9pyl2+fDkAYNSoUUXvoebnVLKZXQiBpUuXomXLlnB3d4ePjw+GDBmCc+fOlRnX3r170aFDB3h4eGDMmDEAgF27dqFr166oXbs23N3dUa9ePQwePBi5ubkV/hzWrFmDnj17IjAwEO7u7mjSpAmmT5+OO3fuaO03atQo1KhRA7/99ht69uwJLy8vdOvWDQBw//59vPnmm4iMjISrqyv8/PwwevRoXL9+vcJzlzzuqVOn0KtXL3h6eiIwMLCoS/XgwYPo1KkTPD09ERERgS+//FLnGL///jsGDhwIHx8fuLm5oWXLljr7aX5WX331FaZNm4bg4GC4uroW/ex37NiBbt26oWbNmvDw8EDHjh2xc+fOSuMHgH/++QfTpk3DAw88AFdXV/j7+6Nv3744depU0T6Wdu2Wpbzr+eeff0Z8fDxq164NNzc3NGjQAJMnTy56/MyZMxg9ejQaNWoEDw8PBAcHIz4+Hr/99pvWsWNiYgAAo0ePLvodKNmNdOTIEQwYMAC+vr5wc3NDq1atsHbtWp04Dx48iI4dO8LNzQ1BQUGYMWOGXl3clf0eanz11Vdo0qQJPDw80KJFC2zevFnrcX1eb8n3c/Xq1Zg5cyaCgoJQs2ZNdO/eHX/++WeFsaakpBQ9t7QVK1ZApVLh8OHDRds2bdqE9u3bw8PDA15eXujRo4dOy8moUaMQHh6uczx9upO6du2KLVu24MKFC1rvXWnvv/9+0fXZvn17HDx4sNJzhYeHo3///vjhhx/QunVruLu7IzIyEp9//nmFMQFARkYG2rRpg0aNGuH06dOV7m9ygszm+vXrAoCYNWtWmY+PGTNG+Pv762z/448/hLe3t2jevLlYsWKF+PHHH8W0adOEg4ODeOONN3T2Hz9+vKhTp44oLCws2rZ79+4Kz62xf/9+AUC88sorer2mu3fviujoaOHp6Sneffdd8eOPP4rXXntNODk5ib59+2rtC0CEhYWJDh06iKSkJLFhwwYREREhfH19xZQpU8TAgQPF5s2bxcqVK0XdunVFdHS01mvo0qWLqF27tggKChKLFy8W27ZtExMnThQAxIQJE3Re6+7du7XOn5aWJgCIL774QgghxJkzZ8SQIUMEAHHgwIGi271794QQQjzzzDPC2dlZTJs2Tfzwww9i1apVIjIyUtStW1dcvXpVKy5fX18RGhoqlixZInbv3i327Nkj0tLShJubm+jRo4fYuHGjSE5OFitXrhQjRowQt27dqvB9nTt3rvjggw/Eli1bRHJysvj4449F/fr1RVxcnNZ+I0eOFM7OziI8PFzMmzdP7Ny5U2zbtk2o1WrRu3dv4enpKWbPni22b98uPvvsMxEcHCyioqJEbm5uhecfOXKkcHFxEU2aNBGLFi0S27dvF6NHjxYAxIwZM0RERIT473//K7Zt2yb69+8vAIgjR44UPf/UqVPCy8tLNGjQQKxYsUJs2bJFPP744wKAmD9/vs7PKjg4WAwZMkRs2rRJbN68WWRmZoqvvvpKqFQq8cgjj4ikpCTx3Xffif79+wtHR0exY8eOCuPPzs4WTZs2FZ6enmLOnDli27ZtYv369WLSpEli165dQgjLvHY15yr5e1rW9fzDDz8IZ2dnER0dLZYvXy527dolPv/8c/HYY48V7bNnzx4xbdo0sW7dOrFnzx6xYcMG8cgjjwh3d3dx6tQpIYQQWVlZ4osvvhAAxKuvvlr0O5Ceni6EEGLXrl3CxcVFxMbGijVr1ogffvhBjBo1Suv3SAj5N8rDw0NERUWJ1atXi2+//Vb06tVL1KtXTwAQaWlp5f6sKvs9BCDCw8PFgw8+KNauXSu2bt0qunbtKpycnMTZs2cNer0l38/w8HAxfPhwsWXLFrF69WpRr1490ahRI1FQUFBurEII0apVK9GxY0ed7TExMSImJqbo+5UrVwoAomfPnmLjxo1izZo1ok2bNsLFxUWkpKQU7Tdy5EgRFhamc7xZs2aJyj6i//jjD9GxY0cREBCg9d4JUfz3Ljw8XPTu3Vts3LhRbNy4UTRv3lz4+PiIf/75p8JzhYWFiZCQEBEVFSVWrFghtm3bJh599FEBQOzZs6doP831c/jwYSGEEL/99psIDQ0V7du3F9evX68wfnNhomNGlSU6TZo0Eb1799bZ3qtXLxESEiKysrK0tj///PPCzc1N3Lx5U2v7p59+KgCIkydPFm1LTk4Wjo6OYvbs2RXG+PXXXwsA4uOPP9brNX388ccCgFi7dq3W9vnz5wsA4scffyzaBkAEBASI27dvF23buHGjACBatmyp9cGwcOFCAUD8+uuvRdu6dOkiAIhvv/1W61zPPPOMcHBwEBcuXBBC6J/oCCHEhAkTyvxjcuDAAQFAvPfee1rb09PThbu7u3j55Zd14tq5c6fWvuvWrRMAxPHjx3WOb4jCwkKRn58v9uzZIwCI1NTUosdGjhwpAIjPP/9c6zmrV68WAMT69eu1th8+fFgAEEuXLq3wnJrjlnx+fn6+8PPzEwDEL7/8UrQ9MzNTODo6iqlTpxZte+yxx4Srq6u4ePGi1nH79OkjPDw8iv7Ian5WnTt31trvzp07wtfXV8THx2ttV6vVokWLFuLBBx+sMP45c+YIAGL79u3l7mOJ167mXJUlOg0aNBANGjQQd+/erfB9KKmgoEDcv39fNGrUSEyZMqVou+aaKPl7oREZGSlatWol8vPztbb3799fBAYGCrVaLYQQYujQocLd3V3rH4CCggIRGRlZaaIjRPm/h0LI96Nu3boiOzu7aNvVq1eFg4ODmDdvnsGvV/N+lk5m165dW5RsVUTzwX7s2LGibYcOHRIAxJdffimEkNdpUFCQaN68edF7JIQQOTk5wt/fX3To0KFoW3USHSGE6NevX5nP1/y9a968uVbypol19erVFZ4rLCxMuLm5aV2bd+/eFb6+vuLZZ5/VeT8OHz4stm/fLmrWrCmGDBli0LVpauy6siBXrlyBv7+/1rZ79+5h586dGDRoEDw8PFBQUFB069u3L+7du6fTDKk5xuXLl4u2denSBQUFBTo1N9W1a9cueHp6YsiQIVrbR40aBQA63QxxcXFa9UdNmjQBAPTp00er6VSz/cKFC1rP9/LywoABA7S2DRs2DIWFhdi7d2/1XkwJmzdvhkqlwhNPPKH1ngcEBKBFixY63Qg+Pj54+OGHtba1bNkSLi4uGDt2LL788kudLq+KnDt3DsOGDUNAQAAcHR3h7OyMLl26AABOnjyps//gwYN14q9Vqxbi4+O14m/ZsiUCAgL0GsGjUqnQt2/fou+dnJzQsGFDBAYGolWrVkXbfX194e/vr/Wz2rVrF7p164bQ0FCtY44aNQq5ubk6zfel49+/fz9u3ryJkSNHasVfWFiI3r174/DhwzrdeCV9//33iIiIQPfu3cvdx1qv3b/++gtnz57FU089BTc3t3L3KygowFtvvYWoqCi4uLjAyckJLi4uOH36dJnXUGlnzpzBqVOnMHz48KLjlfzbk5GRUdTVs3v3bnTr1g1169Yter6joyOGDh2q9+uqSFxcHLy8vIq+r1u3rs41Z+jrLf2ziI6OBqD7cyvt8ccfh7+/v1bpwJIlS+Dn51f0ev/8809cuXIFI0aMgIND8cdsjRo1MHjwYBw8eLDS7mtj6devHxwdHYu+1/d1AvJvWL169Yq+d3NzQ0RERJnP/fLLL9G3b188/fTTWLt2bYXXprmx2s+C3L17V+fiyMzMREFBAZYsWYIlS5aU+bwbN25ofa85xt27dw2OQXNRp6Wl6bV/ZmYmAgICdPp3/f394eTkhMzMTK3tvr6+Wt+7uLhUuP3evXta20v+IdUICAgoisVY/v77bwghyjwfADzwwANa3wcGBurs06BBA+zYsQPvvPMOJkyYgDt37uCBBx7AxIkTMWnSpHLPffv2bcTGxsLNzQ1vvvkmIiIi4OHhgfT0dCQkJOj8XD08PFCzZk2d+P/555+i97G00tdMWTw8PHSuR800A6W5uLho/awyMzPLfE+CgoKKHi+p9L5///03AOgkISXdvHmz3KL969eva/2BLou1XruaGqvKikenTp2KxMREvPLKK+jSpQt8fHzg4OCAp59+Wq+/DZqfwYsvvogXX3yxzH0015HmvSytrG1VUbt2bZ1trq6uWq/D0Ndb+piurq4AKv+76erqimeffRbvvfceFixYgPz8fKxduxZTp04tOobm51ne70BhYSFu3bpllkEDVX2dZT1X8/yynvv111/D3d0dTz/9dJWHqpsKEx0LUqdOHdy8eVNrm4+PDxwdHTFixAhMmDChzOfVr19f63vNMerUqWNwDG3btoWvry++/fZbzJs3r9ILtnbt2vj5558hhNDa99q1aygoKKhSDBXR/PEt6erVq0WxAMWJXumCUn0+3DXq1KkDlUqFlJSUoj8MJZXeVt77FBsbi9jYWKjVahw5cgRLlizB5MmTUbduXTz22GNlPmfXrl24cuUKkpOTi1pxAFlcW5ayzl2nTh3Url0bP/zwQ5nPKfnfsSnUrl0bGRkZOtuvXLlSFF9JpV+D5vElS5aUO+KovCQUAPz8/HDp0qVKY7S0a1cffn5+AFDp6/vf//6HJ598Em+99ZbW9hs3bqBWrVqVnkfz+mfMmIGEhIQy92ncuDEAGb/mtZRU1jZTqe7rNcT48ePx9ttv4/PPP8e9e/dQUFCAcePGFT2u+XmW9zvg4OAAHx8fAPLvVem/VZq4rcnKlSvx2muvoUuXLvjxxx/RsmVLpUMqwq4rCxIZGYmzZ89qbfPw8EBcXByOHTuG6OhotG3bVudW+o/kuXPn4ODgUPRHyBDOzs545ZVXcOrUKcydO7fMfa5du4Z9+/YBALp164bbt29j48aNWvusWLGi6HFjysnJwaZNm7S2rVq1Cg4ODujcuTMAFI1g+PXXX7X2K/08oPz/bvr37w8hBC5fvlzme968eXOD4nZ0dES7du2KmrsrmtRR86FbOpn65JNP9D5f//79kZmZCbVaXWb8Vbk2DNGtW7eihK2kFStWwMPDo9Lh0h07dkStWrVw4sSJMuNv27Ztua1VgOxO+uuvv7Br164KY7S0a1cfERERaNCgAT7//PMyPyA1VCqVzjW0ZcsWrS5toPzfgcaNG6NRo0ZITU0t92egSZjj4uKwc+dOrWROrVZjzZo1er0mQ1oZyqPv6zWGwMBAPProo1i6dCk+/vhjxMfHa7UgNm7cGMHBwVi1ahWEEEXb79y5g/Xr1xeNxALk36tr165pvXf379/Htm3b9IqlvBYWc/P19cWOHTvQpEkTxMXF6ZRUKIktOmbw/fff486dO8jJyQEAnDhxAuvWrQMA9O3bt+iC79q1K+bMmaMzD8qiRYvQqVMnxMbGYvz48QgPD0dOTg7OnDmD7777TueP+cGDB9GyZcui/xgAYM+ePejWrRtef/31Sut0XnrpJZw8eRKzZs3CoUOHMGzYMISGhiIrKwt79+7FsmXLMHv2bHTs2BFPPvkkEhMTMXLkSJw/fx7NmzfHTz/9hLfeegt9+/atsEaiKmrXro3x48fj4sWLiIiIwNatW/Hpp59i/PjxRX9oAgIC0L17d8ybNw8+Pj4ICwvDzp07kZSUpHM8TcIyf/589OnTB46OjoiOjkbHjh0xduxYjB49GkeOHEHnzp3h6emJjIwM/PTTT2jevDnGjx9fYawff/wxdu3ahX79+qFevXq4d+9e0dDMit6XDh06wMfHB+PGjcOsWbPg7OyMlStXIjU1Ve/36bHHHsPKlSvRt29fTJo0CQ8++CCcnZ1x6dIl7N69GwMHDsSgQYP0Pp6hZs2ahc2bNyMuLg6vv/46fH19sXLlSmzZsgXvvPMOvL29K3x+jRo1sGTJEowcORI3b97EkCFD4O/vj+vXryM1NRXXr1/HRx99VO7zJ0+ejDVr1mDgwIGYPn06HnzwQdy9exd79uxB//79ERcXZ5HXrr4SExMRHx+Phx56CFOmTEG9evVw8eJFbNu2DStXrgQgk93ly5cjMjIS0dHROHr0KBYsWKDT5dWgQQO4u7tj5cqVaNKkCWrUqIGgoCAEBQXhk08+QZ8+fdCrVy+MGjUKwcHBuHnzJk6ePIlffvkF33zzDQDg1VdfxaZNm/Dwww/j9ddfh4eHBxITEyusoyqpvN/DipLZ0vR9vcYyadIktGvXDgDwxRdfaD3m4OCAd955B8OHD0f//v3x7LPPIi8vDwsWLMA///yjNfP90KFD8frrr+Oxxx7DSy+9hHv37mHx4sWVTjug0bx5cyQlJeGjjz5CmzZt4ODggLZt2xrvhRrAy8sLP/zwAxISEtCjRw9s2rQJcXFxisSiRdFSaDsRFhYmAJR5Kzka4cyZM0KlUumMAhFCVtCPGTNGBAcHC2dnZ+Hn5yc6dOgg3nzzTa39cnJyhIeHh85oIX2Hl5f07bffin79+gk/Pz/h5OQkfHx8RFxcnPj4449FXl5e0X6ZmZli3LhxIjAwUDg5OYmwsDAxY8aMouGhGihjKK1mZMCCBQvKjPebb74p2talSxfRtGlTkZycLNq2bStcXV1FYGCg+Pe//60zKiQjI0MMGTJE+Pr6Cm9vb/HEE0+II0eO6IwuycvLE08//bTw8/MTKpVK52fy+eefi3bt2glPT0/h7u4uGjRoIJ588kmtodSauEo7cOCAGDRokAgLCxOurq6idu3aokuXLmLTpk2Vvvf79+8X7du3Fx4eHsLPz088/fTT4pdfftGJf+TIkcLT07PMY+Tn54t3331XtGjRQri5uYkaNWqIyMhI8eyzz4rTp09XeP7yjlveaw0LCxP9+vXT2vbbb7+J+Ph44e3tLVxcXESLFi10RvaU9XMuac+ePaJfv37C19dXODs7i+DgYNGvX79y9y/p1q1bYtKkSaJevXrC2dlZ+Pv7i379+mkNNbbEa7f072l5owgPHDgg+vTpI7y9vYWrq6to0KCB1uiiW7duiaeeekr4+/sLDw8P0alTJ5GSkiK6dOkiunTponWs1atXi8jISOHs7Kxz/tTUVPGvf/1L+Pv7C2dnZxEQECAefvhhnZGZ+/btEw899JBwdXUVAQEB4qWXXhLLli3Ta9RVRb+HZb33QshrbuTIkQa/3vKuubJGZVYmPDxcNGnSpNzHN27cKNq1ayfc3NyEp6en6Natm9i3b5/Oflu3bhUtW7YU7u7u4oEHHhAffvih3qOubt68KYYMGSJq1apV9N6VfD2lr08hdK+x8kZdlf6dFkLovJ+lh5cLIX+egwcPFm5ubmLLli2VvgZTUwlRol2NFKcZJfP9999X6fn//e9/MWnSJKSnp2u16NiCrl274saNG/j999+VDoXIILx2bc+vv/6KFi1aIDExEc8995zS4VAFWKNjYebNm4cdO3Zoza6pr4KCAsyfPx8zZsywuSSHiMgSnD17Frt27cLYsWMRGBhYNB0BWS4mOhamWbNm+OKLL6o0WiE9PR1PPPEEpk2bZoLIiIho7ty56NGjB27fvo1vvvnGqteVsxc20XU1aNAgJCcno1u3bkVFvkREREQ20aIzceLEoiGhRERERBo2keiUnh6ciIiICLCARGfv3r2Ij49HUFAQVCqVzuRdALB06VLUr18fbm5uaNOmDVJSUswfKBEREVkdxScMvHPnDlq0aIHRo0frLOwHAGvWrMHkyZOxdOlSdOzYsWgCqxMnThg8yVZZCgsLceXKFXh5eVnc+hxERERUNiEEcnJyEBQUpLV4alk7WgwAYsOGDVrbHnzwQTFu3DitbZGRkWL69Ola23bv3i0GDx5c6Tnu3bsnsrKyim4nTpwodzI/3njjjTfeeOPNsm/p6ekVfu4r3qJTkfv37+Po0aOYPn261vaePXti//79VTrmvHnzMHv2bJ3t6enpOitAExERkWXKzs5GaGhopTW6Fp3o3LhxA2q1WmeV4rp162rNM9OrVy/88ssvuHPnDkJCQrBhwwbExMSUecwZM2Zg6tSpRd9r3qiaNWsy0SEiIrIylZWdWHSio1H6RQghtLbpu8orIFd6Lb3CLREREdkmxUddVaROnTpwdHTUmSX42rVrOq08RERERKVZdKLj4uKCNm3aYPv27Vrbt2/fjg4dOigUFREREVkLxbuubt++jTNnzhR9n5aWhuPHj8PX1xf16tXD1KlTMWLECLRt2xbt27fHsmXLcPHiRYwbN07BqImIiMgaKJ7oHDlyBHFxcUXfawqFR44cieXLl2Po0KHIzMzEnDlzkJGRgWbNmmHr1q0ICwtTKmQiIiKyEjaxqGd1ZGdnw9vbG1lZWRx1RUREZCX0/fy26BodIiIioupgokNEREQ2S/EaHSIiIrJBajWQkgJkZACBgUBsLODoaPYw7LZFJzExEVFRUeXOoExERERVlJQEhIcDcXHAsGHya3i43G5mLEZmMTIREZHxJCUBQ4YApdMLzYoG69YBCQnVPg2LkYmIiMi81Gpg0iTdJAco3jZ5stzPTJjoEBERkXGkpACXLpX/uBBAerrcz0xYjExERERVl5sLHDoE7Nsnu6X0kZFh2phKYKJDRERkj6o6KuryZWD/fpnY7NsHHD8OFBQYdu7AwCqFXBVMdIiIiOxNUpKspSnZzRQSAixapF0orFYDv/2mndhcuKB7vOBgoGNHoH17YN484Pr1sut0VCp5nthY47+mcjDRISIisifljYq6fFluf+01mZDs3w8cPAjk5Gjv5+AAtGgBdOggk5uOHYHQ0OJRVfXqyeOoVNrn0Dy+cKFZ59Ph8HIOLyciInuhVsv5bCoqGC7Ny0u21HTsKJObdu3ktoqU1WIUGiqTHCMMLQf0//xmiw4REZG9qGxUlEb37jIh6dgRaNrU8BaYhARg4ECLmBmZiQ4REZE9UKuB777Tb98xY4DHH6/e+Rwdga5dq3cMI2CiQ0REZMtu3QL++18gMRE4f16/55hxVJSpMdEhIiKyRb//DixZAnz1FXD3rtzm6wvk5wO3b1vMqChTs9uZkbmoJxER2Ry1Gvj2W6BbN6B5c2DZMpnktGghW3UuXQKWL5f7akZBaSg0KsrUOOqKo66IiMjaldU95eAADBoETJwoW2hKJjZmGBVlahx1RUREZO0qm734jz+Ku6dyc+U2X1/gmWeA556Tc9qUxYJGRZkaEx0iIiJLVN7sxR98ALi4AIsXAzt3Fj8WHS1bb4YNA9zdKz++hYyKMjUmOkRERJamvNmLL10CHn20+HsHB+CRR2SC07mzbt0NMdEhIiKyKGq1bMmpqIRWpQJefBGYMAEICzNfbFbIbkddERERWSR9Zi8WAujbl0mOHpjoEBERWZKMDOPuZ+eY6BAREVkSfWcltqHZi02JiQ4REZEliY2Vo6vKKyxWqeScNzY0e7EpMdEhIiKyJI6OwKJF5S/RANjc7MWmxESHiIjI0iQkAH366G4PCQHWrbOa2YstAYeXExERWSJNsfGsWUDjxjY9e7EpMdEhIiKyNFlZQGqqvD9uHBAQoGw8VoxdV0RERJZm/35Zo9OwIZOcarLbRCcxMRFRUVGIiYlROhQiIiJtKSnyK0dWVZvdJjoTJkzAiRMncPjwYaVDISIi0qZJdDp1UjYOG2C3iQ4REZFFuncPOHRI3meLTrUx0SEiIrIkR44A9+8DdevKGh2qFiY6RERElqRkfU55syOT3pjoEBERWRLW5xgVEx0iIiJLoVYD+/bJ+6zPMQomOkRERJbi99+B7GzAywto0ULpaGwCEx0iIiJLoem26tCBSz0YCRMdIiIiS8GJAo2OiQ4REZElEIKFyCbARIeIiMgSpKXJFcudnYEHH1Q6GpvBRIeIiMgSaFpzYmIAd3dlY7EhTHSIiIgsAetzTIKJDhERkSVgfY5JMNEhIiJS2t9/A3/9JZd86NhR6WhsChMdIiIipWlmQ27WDPDxUTYWG2O3iU5iYiKioqIQExOjdChERGTvWJ9jMnab6EyYMAEnTpzA4cOHlQ6FiIjsHetzTMZuEx0iIiKLkJMDHDsm77NFx+iY6BARESnp4EGgsBAIDwdCQpSOxuYw0SEiIlIS63NMiokOERGRklifY1JMdIiIiJRy/77sugLYomMiTHSIiIiUcvQocO8eUKcOEBmpdDQ2iYkOERGRUn76SX7t1EnOikxGx0SHiIhIKSxENjkmOkREREooLNRu0SGTYKJDRESkhBMngFu3AA8PoFUrpaOxWUx0iIiIlKDptmrfHnB2VjYWG8ZEh4iISAmabivW55gUEx0iIiIlcKJAs2CiQ0REZG4XLgDp6YCTE/DQQ0pHY9OY6BAREZmbpjWndWvA01PZWGwcEx0iIiJzY32O2TDRISIiMjfW55gNEx0iIiJzysyUc+gATHTMwG4TncTERERFRSEmJkbpUIiIyJ5ouq2aNJGLeZJJ2W2iM2HCBJw4cQKHDx9WOhQiIrInXN/KrOw20SEiIlIEC5HNiokOERGRudy5Axw9Ku+zPscsmOgQERGZy88/AwUFQEgIEBamdDR2gYkOERGRuZSsz1GplI3FTjDRISIiMhcWIpsdEx0iIiJzyM8HDh6U91mfYzZMdIiIiMzh+HFZjOzjAzRtqnQ0doOJDhERkTlouq06dgQc+PFrLnyniYiIzIH1OYpgokNERGRqQhRPFMj6HLNiokNERGRqf/4J3LgBuLkBbdsqHY1dYaJDRERkappuq3btABcXZWOxM0x0iIiITI31OYpxMvQJ58+fR0pKCs6fP4/c3Fz4+fmhVatWaN++Pdzc3EwRIxERkXVjoqMYvROdVatWYfHixTh06BD8/f0RHBwMd3d33Lx5E2fPnoWbmxuGDx+OV155BWFcv4OIiEi6dAk4f14OKX/oIaWjsTt6JTqtW7eGg4MDRo0ahbVr16JevXpaj+fl5eHAgQP4+uuv0bZtWyxduhSPPvqoSQImIiKyKprRVi1bAjVrKhqKPdIr0Zk7dy769etX7uOurq7o2rUrunbtijfffBNpaWlGC5CIiMiqsdtKUXolOhUlOaXVqVMHderUqXJARERENoWJjqIMHnW1detWbNu2TWf7tm3b8P333xslKCIiIptw6xbw++/yPicKVITBic706dOhVqt1tgshMH36dKMERUREZBP275ezIjdqBNStq3Q0dsngROf06dOIiorS2R4ZGYkzZ84YJSgiIiKbwG4rxRmc6Hh7e+PcuXM628+cOQNPT0+jBEVERGQTmOgozuBEZ8CAAZg8eTLOnj1btO3MmTOYNm0aBgwYYNTgTCkxMRFRUVGIiYlROhQiIrJFd+8Chw/L+6zPUYxKCCEMeUJWVhZ69+6NI0eOICQkBABw6dIlxMbGIikpCbVq1TJFnCaTnZ0Nb29vZGVloSbnNyAiImPZuxfo0gUICACuXAFUKqUjsin6fn4bvASEt7c39u/fj+3btyM1NRXu7u6Ijo5G586dqxUwERGRTSnZbcUkRzEGJzoAoFKp0LNnT/Ts2dPY8RAREdkG1udYBK5eTkREZGxqtRxaDrA+R2FMdIiIiIwtNRXIyZFrW0VHKx2NXWOiQ0REZGyahTw7dAAcHZWNxc4ZlOgUFBTgyy+/xNWrV00VDxERkfVjfY7FMCjRcXJywvjx45GXl2eqeIiIiKybEEx0LIjBXVft2rXD8ePHTRAKERGRDThzBvj7b8DFBeCktIozeHj5c889h6lTpyI9PR1t2rTRWfYhmkVXRERkzzT1OTExgJubsrGQ4YnO0KFDAQATJ04s2qZSqSCEgEqlKnNlcyIiIrvBbiuLYnCik5aWZoo4iIiIbAMTHYticKITFhZmijiIiIis39WrskZHpZJDy0lxVVoC4uzZs1i4cCFOnjwJlUqFJk2aYNKkSWjQoIGx4yMiIrIemtac5s0BK1vk2lYZPOpq27ZtiIqKwqFDhxAdHY1mzZrh559/RtOmTbF9+3ZTxEhERGQdNIXI7LayGAa36EyfPh1TpkzB22+/rbP9lVdeQY8ePYwWHBERkVVhfY7FUQkhhCFPcHNzw2+//YZGjRppbf/rr78QHR2Ne/fuGTVAU8vOzoa3tzeysrJQs2ZNpcMhIiJrlZ0N+PgAhYXApUtAcLDSEdk0fT+/De668vPzK3PCwOPHj8Pf39/QwxEREdmG/ftlklO/PpMcC2Jw19UzzzyDsWPH4ty5c+jQoQNUKhV++uknzJ8/H9OmTTNFjERERJaP9TkWyeBE57XXXoOXlxfee+89zJgxAwAQFBSEN954Q2sSQSIiIrvC+hyLpFeNzqZNm9CnTx84Oztrbc/JyQEAeHl5mSY6M2CNDhERVVteHuDtLb+eOgU0bqx0RDbPqDU6gwYNwj///AMAcHR0xLVr1wDIBMeakxwiIrJxajWQnAysXi2/mmqZoiNHZJLj5wdERJjmHFQleiU6fn5+OHjwIAAUrWlFRERk0ZKSgPBwIC4OGDZMfg0Pl9uNTdNt1amTnBWZLIZeNTrjxo3DwIEDoVKpoFKpEBAQUO6+XNSTiIj0plbLJCEjAwgMlPUtjo7VP25SEjBkCFC6OuPyZbl93TogIaH659FgIbLF0nsenVOnTuHMmTMYMGAAvvjiC9QqZ2rrgQMHGjM+k2ONDhGRQpKSgEmT5JwzGiEhwKJF1UtC1GrZclPyuCWpVPI8aWnGSaoKC4HatYF//gEOHwbatq3+MalS+n5+6z3qKjIyEpGRkZg1axYeffRReHh4GCVQIiKycKZodTFGi8vduzKmK1fkTXP/l1/KT3IAec70dOCTT4AnngCq80+uWg0sXy6THDc3ucYVWRSDZ0a2NWzRISKqgClaXfRpcQkIAL7+Gvj777KTmStXZHJhDKGhQNOmQFRU8deoqMoTIFO1SJFe9P38ZqLDRIeIqGzltbpoim3LanURQo4+ysmRSyJovpa8f/w48PHHxonR3R0ICpItTUFB8pabCyxbVvlzfX2BmzfLfzwkRDsBatoUaNJEDiOvyntDRsVER09MdIiIylBZqwsgu2patwZu39ZOZvLzjRODvz/QqFFxAlMymdHc9/bWHeWkif3yZd1EBNCu0cnKAk6eBP74AzhxQn794w/ZclSe4GDgxg2Z0JXF2DVAVCYmOnpiokNEVIbkZDkcuzo8PWX3j5eX/Kq5n5sLbN9e+fN37wa6dq3auTUtLoB2sqNvi8utWzLx0SQ/mq9XrugfQ3Xip0oZvRi5LPfu3YObm1t1DkFERJaoohaNkiZPBvr21U1matQovzVD3xaX6gzVTkiQyUxZNTQLF1bereTjA3TsKG8l/fOPrMF5443KY9D3PSSTMnj18sLCQsydOxfBwcGoUaMGzp07B0CugfXf//7X6AGaSmJiIqKiohATE6N0KERElicwUL/9Bg4EevQAHnpI1rKEhMjupIq6bBwdZbIA6HY7ab5fuLD63T4JCcD587JlZdUq+TUtrXq1M7VqAV266Levvu8hmZTBic6bb76J5cuX45133oGLi0vR9ubNm+Ozzz4zanCmNGHCBJw4cQKHDx9WOhQiIsvTrh1Qan1DLSqVHK1U1VYXTYtLcLD29pAQ4xbyOjrK7qPHH5dfjVEzExsr4yxvBuTqvjdkVAYnOitWrMCyZcswfPhwOJa4YKKjo3Hq1CmjBkdERAqZObO4qNhUrS6maHExB3O1SJFRGJzoXL58GQ0bNtTZXlhYiHxjVdoTEZFyNm0CPvhA3n/lFdO2upiixcUczNUiRdVmcDFy06ZNkZKSgrCwMK3t33zzDVq1amW0wIiISAEXLwKjRsn7kycDb78N/Oc/plmPytolJMgaJb43Fs3gRGfWrFkYMWIELl++jMLCQiQlJeHPP//EihUrsHnzZlPESERE5pCfDzz2mBxaHRMDzJ8vt2taXUgX3xuLZ3DXVXx8PNasWYOtW7dCpVLh9ddfx8mTJ/Hdd9+hR48epoiRiIjM4dVXgQMH5KipNWuAEgNOiKwVJwzkhIFERMDWrUC/fvL++vWsMSGLZ5YJA2/fvo3CwkKtbUwWiIiszKVLwJNPyvvPP88kh2yKwV1XaWlp6NevHzw9PeHt7Q0fHx/4+PigVq1a8PHxMUWMRERkKgUFcsRTZqZct+rdd5WOiMioDG7RGT58OADg888/R926daEqb8IkIiKyfLNmAT/9JJdtWLMGcHVVOiIiozI40fn1119x9OhRNG7c2BTxEBGRufz4IzBvnrz/6adAGXOkEVk7g7uuYmJikJ6ebopYiIjIXDIygCeekItqPvssMHSo0hERmYTBLTqfffYZxo0bh8uXL6NZs2ZwLrUWSnR0tNGCIyIiE1CrgeHDgevXgejo4lmQiWyQwYnO9evXcfbsWYwePbpom0qlghACKpUKarXaqAESEZGRzZ0r15Xy9ATWrgXc3ZWOiMhkDE50xowZg1atWmH16tUsRiYisja7dgFz5sj7n3wCsN6SbJzBic6FCxewadOmMhf2JCIiC/b337LLSgjgqafkfSIbZ3Ax8sMPP4zU1FRTxEJERKaiVsvi46tXgaZNgcWLlY6IyCwMbtGJj4/HlClT8Ntvv6F58+Y6xcgDBgwwWnBERGQkb78N7NgBeHjIuhwPD6UjIjILg9e6cnAovxHIGouRudYVEdm8vXuBuDigsBD44gtg1CilIyKqNpOtdVV6bSsiIrJg16/LJR4KC+V6VkxyyM4YXKNDRERWQpPcXLkCREYCiYlKR0Rkdnq16CxevBhjx46Fm5sbFldSwDZx4kSjBEZERNW0YAHwww+Am5usy6lRQ+mIiMxOrxqd+vXr48iRI6hduzbq169f/sFUKpw7d86oAZoaa3SIyCbt2wd06SJHW336KfD000pHRGRURq3RSUtLK/M+ERFZoMxMWZejVsuvTz2ldEREijG4RmfOnDnIzc3V2X737l3M0cy2SURE5qNWA8nJwOrVcmmHUaOA9HSgUSM5+zFnsCc7ZvDwckdHR2RkZMDf319re2ZmJvz9/Tm8nIjInJKSgEmTgEuXtLc7OQGHDwMtWyoSFpGp6fv5bXCLjmbxztJSU1Ph6+tr6OGIiKiqkpKAIUN0kxwAKCgArKxmksgU9J5Hx8fHByqVCiqVChEREVrJjlqtxu3btzFu3DiTBElERKWo1bIlp7xGeZUKmDwZGDgQcHQ0a2hElkTvRGfhwoUQQmDMmDGYPXs2vL29ix5zcXFBeHg42rdvb5IgiYiolJSUsltyNISQdTopKUDXrmYLi8jS6J3ojBw5EoAcat6xY0c4ORk8qTIRERlLRoZx9yOyUXrV6Ny5c6fofpcuXSpNckruT0REJhAYaNz9iGyUXolOw4YN8dZbb+HKlSvl7iOEwPbt29GnT59KZ08mIqJqio0FQkLKHzquUgGhoXI/IjumV/9TcnIyXn31VcyePRstW7ZE27ZtERQUBDc3N9y6dQsnTpzAgQMH4OzsjBkzZmDs2LGmjpuIyL45OgKLFgGDB+s+pkl+Fi5kITLZPYPm0bl06RK++eYb7N27F+fPn8fdu3dRp04dtGrVCr169ULfvn3h4GBd64RyHh0ismo9ewLbt2tvCw2VSU5CgiIhEZmDvp/fBk8YaGuY6BCR1bp/X9bg3LwJvPsuEBQkv4+NZUsO2TyjrnVFREQWaPt2meTUrSvnzGFyQ6TDuvqZiIio2OrV8uvQoUxyiMphUKJz9OhR9O/fH7dv39bafvfuXQwaNAj79+83anBERFSO3Fxg40Z5//HHFQ2FyJIZlOhERkYiOzsb3bt3R1ZWFgAgNzcXffr0QUZGBpo2bWqSIImIqJTvvgPu3AHq1wfatVM6GiKLZVCi4+npiR9++AGenp6Ii4vD2bNn0b17d+Tn5+PHH3/UWhaCiIhMaNUq+fXxx8ufS4eIDK/R8fDwwJYtW+Dn54cmTZrA2dkZ27Zt44glIiJzuXUL+P57eX/YMGVjIbJwVSpGzs3NxfXr11GjRg1kZmYiNzfX2HEREVF5kpKA/HygeXOAJQNEFTI40bl+/Tq6dOmCgIAAXLhwAZGRkYiNja1weQgiIjKikt1WRFQhgxKdjIwMdO7cGQ0aNMDGjRvh5eWFNWvWoE2bNujcuTMuXrxoqjiNLjExEVFRUYiJiVE6FCIi/WVkALt3y/uPPaZsLERWwKBE5/bt2+jRowfWr18PFxcXAICjoyNWrlyJgQMHIicnxyRBmsKECRNw4sQJHD58WOlQiIj0t2YNIATQvr0ccUVEFeISEFwCgoisSbt2wKFDwJIlwPPPKx0NkWL0/fzmzMhERNbizBmZ5Dg4AI8+qnQ0RFahymtd3b9/H9euXUNhYaHW9nr16lU7KCIiKsPXX8uv3bvL9a2IqFIGJzqnT5/GmDFjdJZ7EEJApVJBrVYbLTgiIvp/QnC0FVEVGJzojBo1Ck5OTti8eTMCAwOh4oycRESm9+uvwMmTgKsrMGiQ0tEQWQ2DE53jx4/j6NGjiIyMNEU8RERUFs1K5f36AVxuh0hvBhcjR0VF4caNG6aIhYiIylJYWJzosNuKyCAGJzrz58/Hyy+/jOTkZGRmZiI7O1vrRkRERnbgAHDxIuDlJVt0iEhvBnddde/eHQDQrVs3re0sRiYiMhFNEfKgQYC7u7KxEFkZgxOd3Zqpx4mIyPQKCoBvvpH3uVI5kcEMTnS6dOliijiIiKgsO3cC168Dfn5AqZZ0IqocZ0YmIrJkmm6rf/0LcKryHK9EdouJDhGRpbp7F9iwQd7naCuiKmGiQ0RkqbZsAXJygHr15GrlRGQwvROd27dvmzIOIiIqreTcOQ78v5SoKvT+zalTpw769OmDjz76CJcvXzZlTERElJUlW3QAdlsRVYPeic6ff/6Jvn37Yv369XjggQcQExODuXPn4tdffzVlfERE9mnDBiAvD4iKAqKjlY6GyGrpneiEhYXhhRdewI4dO3Dt2jVMnToVf/zxBzp37oz69etj0qRJ2LVrFycMJCIyhpLdVlw8majKVEIIUZ0DFBQUYNeuXfjuu++wadMm5OTkYMmSJRg+fLixYjSp7OxseHt7IysrCzVr1lQ6HCIi4O+/gaAgucbVmTNAgwZKR0RkcfT9/K72pAxOTk7o2bMnevbsiSVLluDYsWMoKCio7mGJiOzXN9/IJOfBB5nkEFWT0WefatWqlbEPSURkXzSTBLIImajaOF6RiMiSpKXJ1cpVKmDoUKWjIbJ6THSIiCzJ11/Lr3FxQGCgsrEQ2QAmOkRElkQz2oorlRMZBRMdIiJL8fvvwG+/Ac7OQEKC0tEQ2QSDipH//PNPrF69GikpKTh//jxyc3Ph5+eHVq1aoVevXhg8eDBcXV1NFSsRkW3TtOb07Qv4+CgbC5GN0KtF59ixY+jRowdatGiBvXv3IiYmBpMnT8bcuXPxxBNPQAiBmTNnIigoCPPnz0deXp6p4yYisi1CaE8SSERGodeEgWFhYXjppZcwbNgw+Pr6lrvfgQMH8MEHH6Bly5b497//bdRATYUTBhKRRTh4UK5Q7ukJXLsGeHgoHRGRRTPqhIGnT5+Gi4tLpfu1b98e7du3x/379/WPlIiIiltzHnmESQ6REenVdaVPklOd/YmI7FpBAbBmjbzPbisio6rSqKudO3eif//+aNCgARo2bIj+/ftjx44dxo6NiMg+JCfL9a1q1wZ69lQ6GiKbYnCi8+GHH6J3797w8vLCpEmTMHHiRNSsWRN9+/bFhx9+aIoYiYhsm6bbasgQObSciIzG4NXLg4ODMWPGDDz//PNa2xMTE/Gf//wHV65cMWqApsZiZCJSVF4eULcukJUF7NkDdO6sdEREVkHfz2+DW3Sys7PRu3dvne09e/ZEdna2oYcjIrJv338vk5yQEKBTJ6WjIbI5Bic6AwYMwIYNG3S2f/vtt4iPjzdKUEREdkOzUvnQoYADJ6snMja9hpcvXry46H6TJk3wn//8B8nJyWjfvj0A4ODBg9i3bx+mTZtmmiiJiGxRTg7w3XfyPte2IjIJvWp06tevr9/BVCqcO3eu2kGZE2t0iEgxX30FPPkkEBEBnDoFqFRKR0RkNYw6YWBaWprRAiMiov9XcqVyJjlEJsEOYSIiJVy/Dvz4o7zPSQKJTIaJDhGREtatA9RqoE0b2XVFRCbBRIeISAlcqZzILJjoEBGZ28WLQEqKrMsZOlTpaIhsmkGJTkFBAWbPno309HRTxUNEZPs0C3h27iwnCiQikzEo0XFycsKCBQugVqtNFQ8Rke3TTBLIbisikzO466p79+5ITk42QShERHbg5Eng+HHAyUku4klEJqXXPDol9enTBzNmzMDvv/+ONm3awNPTU+vxAQMGGC04IiKboylC7tULqF1b2ViI7IDBq5c7VLAWi0qlsrpuLc6MTERmI4QcSn7mDLByJZd9IKoGo86MXFJhYWG1AiMisjtqtRxltW+fTHLc3AC2fhOZBYeXExGZUlISEB4OxMUBr74qt6lUxbMiE5FJVSnR2bNnD+Lj49GwYUM0atQIAwYMQEpKirFjIyKybklJsuD40iXt7Xfvyu1JScrERWRHDE50/ve//6F79+7w8PDAxIkT8fzzz8Pd3R3dunXDKs2QSSIie6dWA5Mmybqc8kyeLPcjIpMxuBi5SZMmGDt2LKZMmaK1/f3338enn36KkydPGjVAU2MxMhGZRHKy7K6qzO7dQNeupo6GyObo+/ltcIvOuXPnEB8fr7N9wIABSEtLM/RwRES2KSPDuPsRUZUYnOiEhoZi586dOtt37tyJ0NBQowRFRGT1AgONux8RVYnBw8unTZuGiRMn4vjx4+jQoQNUKhV++uknLF++HIsWLTJFjERE1ic2FggOBi5fLvtxlUqucxUba964iOyMwYnO+PHjERAQgPfeew9r164FIOt21qxZg4EDBxo9QCIiq+TgICcHLCvRUank14ULAUdHs4ZFZG/0SnQWL16MsWPHws3NDRcvXsQjjzyCQYMGmTo2IiLr9emnstBYpQLq1AGuXy9+LCREJjkJCYqFR2Qv9Bp15eTkhCtXrsDf3x+Ojo7IyMiAv7+/OeIzOY66IiKjO3wY6NQJuH8fmD8fmDZNzoyckSFrcmJj2ZJDVE1GXQIiKCgI69evR9++fSGEwKVLl3Dv3r0y961Xr17VIiYisgU3bgCDB8skZ9Ag4KWXZKsOh5ATKUKvFp1ly5bhhRdeQEFBQbn7CCG4qCcR2Te1GujdG9ixQ9bnHDoEeHsrHRWRTTJqi87YsWPx+OOP48KFC4iOjsaOHTtQu3ZtowVLRGagWVjSVN0npj6+NXj9dZnkeHjI5R2Y5BApTu9RV15eXmjWrBm++OILdOzYEa6urqaMi4iMKSlJLkdQcs2lkBBg0SLjFMSa+vjWYNMm4K235P3//hdo2lTZeIgIQBWWgLA17Loim6dZWLL0r7pmiPO6ddVLRkx9fGtw+jTQti2QnS0TvoULlY6IyObp+/nNRIeJDtkytRoID9ddPbuk2rWB5csBLy/Z5VL65u4OOJXT+FvZ8TWT4qWl2W431p07QPv2wG+/AR07yiHlzs5KR0Vk84xao0NEViolpeIkBwAyM4Ey1q/T4uJSdhJ0717FxxcCSE+XcdjiqCMhgGeflUlOQADwzTdMcogsDBMdIlum74KR4eGAmxuQm1t8u3u3uDvq/n15++cf08ZhbZYuBVaulK1Va9dy3SoiC1TlROf+/ftIS0tDgwYN4FReszYRKUvfD94vvtBtcRECyMvTTn5K344cAebONV4c1mT/fmDyZHl/wQKuWUVkoQyu0cnNzcULL7yAL7/8EgDw119/4YEHHsDEiRMRFBSE6dOnmyRQU2GNDtk0tRoIC6t8Ycmq1tBoanQuX9YtRjbG8S3V338DrVsDV64A//oX8PXXxcXXRGQW+n5+Oxh64BkzZiA1NRXJyclwc3Mr2t69e3esWbOmatFW0+bNm9G4cWM0atQIn332mSIxEFkkR0e5FEFZjLGwpKOjHEJe8nil2drClQUFwGOPySSnSRM5lJxJDpHFMjjR2bhxIz788EN06tQJqhK/3FFRUTh79qxRg9NHQUEBpk6dil27duGXX37B/PnzcfPmTbPHQWSRLlwAvv1W3vf11X4sJMQ4Q78TEuRxgoN1HwsPr7zQ2dr8+99AcjJQo4YcWl+jhtIREVEFDE50rl+/XuaCnnfu3NFKfMzl0KFDaNq0KYKDg+Hl5YW+ffti27ZtZo+DyCK99JIcGdWli+xu2b0bWLVKfk1LM978NgkJwPnzxcdfvx7w8ZHnePtt45zDEqxfL+txADkkPzJS0XCIqHIGJzoxMTHYsmVL0fea5ObTTz9F+/btDQ5g7969iI+PR1BQEFQqFTZu3Kizz9KlS1G/fn24ubmhTZs2SElJKXrsypUrCC7xn2RISAgul1ePQGRPdu2Sw50dHIDFi+VcOF27Ao8/Lr8auzvJ0bH4+AkJwIcfyu1z5wK//mrccynh1Clg1Ch5/8UX5cKdRGTxDE505s2bh5kzZ2L8+PEoKCjAokWL0KNHDyxfvhz/+c9/DA7gzp07aNGiBT7U/FEsZc2aNZg8eTJmzpyJY8eOITY2Fn369MHFixcByMVES1OiZYnIohQUABMnyvvjxwPR0eaP4fHHgYEDgfx8mSDk55s/BmO5fVsmb7dvy2Ru3jylIyIiPRmc6HTo0AH79u1Dbm4uGjRogB9//BF169bFgQMH0KZNG4MD6NOnD958800klNOE/v777+Opp57C008/jSZNmmDhwoUIDQ3FRx99BAAIDg7WasG5dOkSAisYypqXl4fs7GytG5HN+egj4I8/5KzHc+YoE4NKBXz8sawNOnbMeruwhACeego4eRIICpIjrDilBpH1EBYEgNiwYUPR93l5ecLR0VEkJSVp7Tdx4kTRuXNnIYQQ+fn5omHDhuLSpUsiOztbNGzYUNy4caPcc8yaNUsA0LllZWWZ5DURmd21a0LUqiUEIMTHHysdjRArV8pYnJ2FSE1VOhrDvf9+cfz79ysdDRH9v6ysLL0+v6v8b8m1a9dw7do1FBYWam2PNmIT+Y0bN6BWq1G3bl2t7XXr1sXVq1cBAE5OTnjvvfcQFxeHwsJCvPzyy6hdu3a5x5wxYwamTp1a9H12djZCQ0ONFjOR4mbOlDMYt2oFPP200tHILqy1a+Xor1GjgJ9/tp5lEvbulQXdAPD++3JNKyKyKgYnOkePHsXIkSNx8uRJnfoYlUoFtVpttOBKHrckIYTWtgEDBmDAgAF6HcvV1RWurq5GjY/IYhw5AmjmklqyxDLmr9F0YaWkFHdhvfaa0lFVLiMDGDpUToo4fDgwYYLSERFRFRhcozN69GhERERg//79OHfuHNLS0opu586dM2pwderUgaOjY1Hrjca1a9d0WnmI7F5hoSxAFkJ+MHfsqHRExQICZOIFWMcorPx84NFHgatXgebNgU8+4aSARFbK4BadtLQ0JCUloWHDhqaIR4uLiwvatGmD7du3Y9CgQUXbt2/fjoEDB5r8/ERWZeVK4MABwNMTmD9f6Wh0WXIXllotW5wyMuS6XElJwL59QM2acu4cT0+lIySiKjI40enWrRtSU1ONlujcvn0bZ86cKfo+LS0Nx48fh6+vL+rVq4epU6dixIgRaNu2Ldq3b49ly5bh4sWLGDdunFHOT2QTsrOBl1+W9197rexZipWm6cLau9eyurCSkoBJk4BLl3QfW7ECaNTI/DERkdEYnOh89tlnGDlyJH7//Xc0a9YMzqX+I9O3VkbjyJEjiIuLK/peUyg8cuRILF++HEOHDkVmZibmzJmDjIwMNGvWDFu3bkVYWJihoRPZrjfflN0sjRoVr6htiTRdWE88IbuwBg5UZo4fjaQkYMiQshckBWRLDxFZNYNXL9+0aRNGjBiBnJwc3YOZqBjZlLh6OVm9P/+UdST5+cCWLUDfvkpHVDEhgEGDZBdWq1bKdWFpVl4vqyUHsN2V14lshMlWL584cSJGjBiBjIwMFBYWat2sLckhsnpCyG6X/HygXz/LT3KA4i4sHx9lJxJMSSk/yQHke5ueLvcjIqtlcKKTmZmJKVOmcNQTkSXYvBnYtg1wcQE++EDpaPRnCaOwMjKMux8RWSSDE52EhATs3r3bFLEQkSHu3Suux5k61fqKZocNU3YtrIAA/farYEkZIrJ8BhcjR0REYMaMGfjpp5/QvHlznWLkiZqFBC1cYmIiEhMT2d1G1uv994Fz5+T6SzNnKh2N4ZQchZWTAyQmVh5fSAgQG2uemIjIJAwuRq5fv375B1OpjD5poKmxGJmsUno6EBkJ5ObK+XOGDVM6oqpbuVKOwnJ2ljM7m3oU1p9/ymLokydlkbFaLZOakn8KNZMDrlsnVy0nIouj7+d3lSYMJCKFvfyyTHI6dpQT8VmzYcOAb74xz0SCGzYAI0fKFp2gIJnIZGTozqMTEgIsXMgkh8gGGNyiY2vYokNWZ+9eoEsX2epw9Kgcom3trl4FoqKAW7eAOXOM34WlVstjzpsnv+/cGVizprhOp/TMyLGxHFJOZOH0/fzWK9GZOnUq5s6dC09PT62Vv8vy/vvvGx6tgpjokFUpKADatJGjlMaNAz76SOmIjMdUXVg3bshWo+3b5feTJwPvvGM5y08QUZUYtevq2LFjyP//ERHHjh0zToREZLhly2SS4+MjZ0O2JabowvrlF9n9dOEC4OEhV3a39q4+IjIIu67YokPWIjNTDiG/dUuOGHruOaUjMr6MDKBpU/ka584FXn216sf64gtg/HggLw9o0EDW5zRvbrxYiUhRJpsZecyYMWUu/3Dnzh2MGTPG0MMRkb5efVUmANHRwNixSkdjGoGBxRMJzplTtYkE8/JkgjNmjLzfv7/sCmOSQ2SXDG7RcXR0REZGBvz9/bW237hxAwEBASgoKDBqgKbGFh2yCseOydocIYDkZFmMbKuqsxbWpUtykc6ff5bF2m+8IRNEB4P/pyMiC2f0Fp3s7GxkZWVBCIGcnBxkZ2cX3W7duoWtW7fqJD9EZARCABMnyq+PPWbbSQ4gE5SPPipeC2v+fP2el5wsk8GffwZq1ZLLY7z+OpMcIjun9zw6tWrVgkqlgkqlQkREhM7jKpUKs2fPNmpwRARg9Wrgp59kMe2CBUpHYx6aLqwnnpBdWAMGlD8KSwg5581LL8lh4tHRQFKSrMshIrund6Kze/duCCHw8MMPY/369fD19S16zMXFBWFhYQgKCjJJkER26/Zt+QEOyGUeQkKUjcecSo/C2r8fOHhQe66bu3eBp5+Wc+IAMjH65BOZFBIRoQo1OhcuXEC9evWg0kyRbuVYo0MWbcYMuQbUAw8Af/wBuLkpHZF5lRyFVbMmkJ1d/FhAgKzdSU8HnJzk6u0TJhQv30BENs2oNToXL14suh8WFlZpknP58mU9w1ROYmIioqKiEBMTo3QoRGU7fVou3AnID3F7S3IA2XLz5JPyfskkB5CzKaeny3qc3buB559nkkNEOvRKdGJiYvDMM8/g0KFD5e6TlZWFTz/9FM2aNUNSUpLRAjSVCRMm4MSJEzh8+LDSoRAVU6tlUe3q1fID/v59oHdvID5e6ciUoVYD69dXvI+HB9C+vXniISKro1eNzsmTJ/HWW2+hd+/ecHZ2Rtu2bREUFAQ3NzfcunULJ06cwB9//IG2bdtiwYIF6NOnj6njJrI9SUm6i0sCQN++9ttSkZKi+36UduWK3K9rV7OERETWRa8WHV9fX7z77ru4cuUKPvroI0RERODGjRs4ffo0AGD48OE4evQo9u3bxySHqCqSkuT8L2V9qE+aJB+3RxkZxt2PiOwOl4BgMTIpTa0GwsPLb7lQqeRoq7Q0+1tROzkZiIurfL/du9miQ2RnTLYEBBEZWWXdM0LIotuUFPPFZCliY2WSV17XnUoFhIbK/YiIysBEh0hp7J4pn6MjsGiRvF862dF8v3Ch/bV0EZHemOgQKS0w0Lj72ZqEBGDdOiA4WHt7SIjcnpCgTFxEZBVYo8MaHVKaWg2EhQHlzT9lzzU6JanVsvuu5MzI9vx+ENk5fT+/9V4CgohMxNERePhh4KuvdB9j90wxR0cWHBORwdh1RaS0S5eKh4/7+Gg/xu4ZIqJqYYsOkdKmTgXu3AE6dpTDpPftY/cMEZGRMNEhUtL27XKFbgcHIDFRLlLJ7hkiIqNh1xWRUvLy5EKUgPzaooWy8RAR2SAmOkRK+eAD4K+/gLp1gTlzlI6GiMgm2W2ik5iYiKioKMTExCgdCtmjixeBuXPl/XffBby9lY2HiMhGcR4dzqNDShgyBFi/XhYb79ljv6uTExFVEde6IrJU27bJJMfRURYgM8khIjIZJjpE5pSXB7zwgrw/cSLQvLmy8RAR2TgmOkTm9O67wOnTQEAA8MYbSkdDRGTzmOgQmcv588B//iPvv/cewJowIiKTY6JDZC5TpgB378oJAR9/XOloiIjsAhMdInPYuhXYuBFwcgI+/JAFyEREZsJEh8jU7t2ThccAMHky0LSpouEQEdkTJjpEprZgAXD2LBAUBLz+utLREBHZFSY6RKaUlga89Za8//77gJeXsvEQEdkZJjpEpjRpkuy6evhh4F//UjoaIiK7w0SHyFS++07eWIBMRKQYJjpEpnD3rmzNAYCpU4EmTZSNh4jITjHRITKF+fNlfU5ICPDaa0pHQ0Rkt5joEBnb2bPA22/L+x98ANSooWw8RER2jIkOkTEJIefMycsDevQABg9WOiIiIrtmt4lOYmIioqKiEBMTo3QoZEu++07OguzsDCxZwgJkIiKFqYQQQukglJSdnQ1vb29kZWWhJhdZpOrIzQWiooALF4AZM4rnzyEiIqPT9/Pbblt0iIxu3jyZ5ISGAjNnKh0NERGBiQ6RcZw+Dbzzjry/cCHg6aloOEREJDkpHQCRUanVQEoKkJEBBAYCsbGAo6NpzykE8MILwP37QK9ewKBBpj0fERHpjYkO2Y6kJDlJ36VLxdtCQoBFi4CEBNOdd+NGYNs2wMWFBchERBaGXVdkG5KSgCFDtJMcALh8WW5PSjLNee/cASZPlvdffhlo1Mg05yEioiphokPWT62WLTllDSDUbJs8We5nbG+9BVy8CISFyZFWRERkUdh1RdYvJUW3JackIYD0dGDnTqBnz+qdq2QNUEFBcQHyokWAh0f1jk1EREbHRIesU14ecPgwsGcP8PXX+j2nd28gIgJo1kzemjaVXxs1kiuMV6asGiAAaN0aGDDA8NdAREQmxwkDOWGg+VVlZFRuLnDwILB3r0xuDh4E7t0zTjwuLkBkpHby06wZEB4OOPx/766mBqisXxeVCli3zrQFz0REpEXfz28mOkx0zEvfkVE5OcD+/cWJzaFDQH6+9rH8/YEuXYBOnWStzLVr5SciISHyeCdPAr//Xnz74w9ZUFwWDw8503FUFPDtt0BWVtn7aY6flmb6oexERASAiY7emOiYUXmtIprh2NOny2Rmzx7gl190i4eDg2Vi06UL0Lkz0Lhx8XM1xwa0j695vLwWl8JCWUxcMvH5/XeZEOXlGfb6du8GunY17DlERFQlTHT0xETHTNRq2RVUUdFwafXry4RGk9zUr1/xHDVltRaFhsqZig3tViooAM6elYnP6tUyUarMqlXA448bdh4iIqoSfT+/WYxMZTPWDMNqtVweYdUq/ZKcfv2Axx6TiU1oqGHnSkgABg40TtxOTrLFqHFjwNdXv0QnMNDw8xARkUkx0SFdVZ1hODsb+PVXIDVV3o4fl91Ad+/qf+7hw6vXKuLoaPzuo9hY+fovX664Big21rjnJSKiamOiQ9rKq6PRzDC8bp1cy+nCheJkRpPYnDtX9jE9PYF69WTdS2UssVXE0VEmeUOGyKSmrBqghQtZiExEZIFYo8ManWL61NG4uspbdnbZj4eEAC1bAi1ayFvLlkCDBjI5CA+vvFXEkkcuGbMGiIiIqoU1OmS4ymYYBuRIpLw8wNlZzjlTMqGJjgZq1y7/udbeKmLMGiAiIjILJjpULCNDv/3efhuYMkVOtGeIhATZ9VVW/Y+1tIqYogaIiIhMxm4TncTERCQmJkJtioUerZW+9THt2hme5GiwVYSIiMyINTqs0ZGEAN54A5gzp/x9rKGOhoiI7AJrdEh/ajUwcSKwdGnxNmutoyEiIirBQekASGF5eXLemqVLZTLz4YfA+vVyuYWSQkK4cCUREVkdtujYs5wcOSfOzp1yFNX//gf861/yMdbREBGRDWCiY6+uXQP69JGLZ9aoAWzYAHTvXvw4RxcREZENYKJjj9LSgJ49gTNnAD8/YOtWoG1bpaMiIiIyOiY69ubXX4FevYCrV4GwMODHH4GICKWjIiIiMgkWI9uTvXuBzp1lktO8ObB/P5McIiKyaUx07MW338ruqqwsoFMnmfQEBSkdFRERkUkx0bEHn30mh4Xn5QEDBsjuqlq1lI6KiIjI5Jjo2DIhgLfeAp55BigsBMaMkXPkuLsrHRkREZFZMNGxVYWFwOTJwMyZ8vsZM2TLjhPrz4mIyH7wU88W3b8PjBoFrF4tv1+4UK4YTkREZGeY6JiCWq3crMK3bwODB8s6HCcn4MsvgWHDzHNuIiIiC8NEx9iSkmTryaVLxdtCQoBFi0y/TtSNG0C/fsChQ4CHh6zH6d3btOckIiKyYKzRMaakJGDIEO0kBwAuX5bbk5JMd+4LF+Sw8UOHgNq1gV27mOQQEZHdY6JjLGq1bMkRQvcxzbbJk+V+xvbHH0DHjsCffwKhocBPPwHt2hn/PERERFaGiY6xpKTotuSUJASQng4sWSLraKpCrQaSk2WRcXKy/H7/flkDdPkyEBUlv4+MrNrxiYiIbAxrdIwlI0O//aZMAaZNk0nJgw8W35o1A5ydy39eWbU/deoA2dlylFX79sDmzYCvb/VeBxERkQ1homMsgYH67efnB1y/Dvz+u7x9/rnc7uYGtG6tnfw88ACgUhXX/pTuFrtxQ35t3RrYsUMWIBMREVERJjrGEhsrR1ddvlx2nY5KJR9PSwOuXQMOH5aFw5pbVpbsdtq/v/g5vr5ATAxw4EDZx9S4dg1wdTX+ayIiIrJyKiEq+gS1fdnZ2fD29kZWVhZq1qxZvYNpWl4A7cREpZJf160re4h5YSFw5ox24nPsmOyS0tfu3UDXrlUOnYiIyJro+/nNYmRjSkiQyUxwsPb2kJDykxwAcHAAIiKAJ54AFi8GDh4EcnJkq8/o0fqdW98aISIiIjtit4lOYmIioqKiEBMTY9wDJyQA58/LFpZVq+TXtDTDJwt0cQHatgWefFK//fWtESIiIrIj7LoyZteVKajVQHi4frU/5lpmgoiISGHsurIVjo5y+QiguNZHQ/P9woVMcoiIiMrARMcaVLX2h4iIyM5xeLm1SEgABg5UblV0IiIiK8REx5o4OnIIORERkQHYdUVEREQ2i4kOERER2SwmOkRERGSzmOgQERGRzWKiQ0RERDaLiQ4RERHZLCY6REREZLOY6BAREZHNYqJDRERENsvuZ0bWLN6enZ2tcCRERESkL83ntuZzvDx2n+jk5OQAAEJDQxWOhIiIiAyVk5MDb2/vch9XicpSIRtXWFiIK1euwMvLCyqVqsx9YmJicPjwYYMey87ORmhoKNLT01GzZk2jxmxKFb1WSz5XVY9l6PP03V+f/Srbp6zHrfW6Asx3bfG6Mvy6Aqz32rKn66oqzzXWtWWJ15UQAjk5OQgKCoKDQ/mVOHbfouPg4ICQkJAK93F0dCz3B1TRYwBQs2ZNq/qjUdnrsdRzVfVYhj5P3/312a+yfSp63NquK8B81xavq6pfV4D1XVv2dF1V5bnGurYs9bqqqCVHg8XIepgwYUKVHrNG5nw9xjxXVY9l6PP03V+f/Srbh9eW8ufhdWX57Om6qspzjXVtWfN1ZfddV6aSnZ0Nb29vZGVlWdV/R2TZeF2RqfDaIlOwhOuKLTom4urqilmzZsHV1VXpUMiG8LoiU+G1RaZgCdcVW3SIiIjIZrFFh4iIiGwWEx0iIiKyWUx0iIiIyGYx0SEiIiKbxUSHiIiIbBYTHYUMGjQIPj4+GDJkiNKhkI1IT09H165dERUVhejoaHzzzTdKh0Q2ICcnBzExMWjZsiWaN2+OTz/9VOmQyIbk5uYiLCwML774osnOweHlCtm9ezdu376NL7/8EuvWrVM6HLIBGRkZ+Pvvv9GyZUtcu3YNrVu3xp9//glPT0+lQyMrplarkZeXBw8PD+Tm5qJZs2Y4fPgwateurXRoZANmzpyJ06dPo169enj33XdNcg626CgkLi4OXl5eSodBNiQwMBAtW7YEAPj7+8PX1xc3b95UNiiyeo6OjvDw8AAA3Lt3D2q1Gvz/mIzh9OnTOHXqFPr27WvS8zDRqYK9e/ciPj4eQUFBUKlU2Lhxo84+S5cuRf369eHm5oY2bdogJSXF/IGSVTHmdXXkyBEUFhYiNDTUxFGTpTPGdfXPP/+gRYsWCAkJwcsvv4w6deqYKXqyVMa4rl588UXMmzfP5LEy0amCO3fuoEWLFvjwww/LfHzNmjWYPHkyZs6ciWPHjiE2NhZ9+vTBxYsXzRwpWRNjXVeZmZl48sknsWzZMnOETRbOGNdVrVq1kJqairS0NKxatQp///23ucInC1Xd6+rbb79FREQEIiIiTB+soGoBIDZs2KC17cEHHxTjxo3T2hYZGSmmT5+utW337t1i8ODBpg6RrFBVr6t79+6J2NhYsWLFCnOESVamOn+vNMaNGyfWrl1rqhDJClXlupo+fboICQkRYWFhonbt2qJmzZpi9uzZJomPLTpGdv/+fRw9ehQ9e/bU2t6zZ0/s379foajI2ulzXQkhMGrUKDz88MMYMWKEEmGSldHnuvr777+RnZ0NQK5EvXfvXjRu3NjssZL10Oe6mjdvHtLT03H+/Hm8++67eOaZZ/D666+bJB4nkxzVjt24cQNqtRp169bV2l63bl1cvXq16PtevXrhl19+wZ07dxASEoINGzYgJibG3OGSldDnutq3bx/WrFmD6Ojoov7yr776Cs2bNzd3uGQl9LmuLl26hKeeegpCCAgh8PzzzyM6OlqJcMlK6Ps5aC5MdExEpVJpfS+E0Nq2bds2c4dENqCi66pTp04oLCxUIiyychVdV23atMHx48cViIqsXWWfgxqjRo0yaRzsujKyOnXqwNHRUSdrvXbtmk52S6QvXldkCryuyBQs7bpiomNkLi4uaNOmDbZv3661ffv27ejQoYNCUZG143VFpsDrikzB0q4rdl1Vwe3bt3HmzJmi79PS0nD8+HH4+vqiXr16mDp1KkaMGIG2bduiffv2WLZsGS5evIhx48YpGDVZOl5XZAq8rsgUrOq6MslYLhu3e/duAUDnNnLkyKJ9EhMTRVhYmHBxcRGtW7cWe/bsUS5gsgq8rsgUeF2RKVjTdcW1roiIiMhmsUaHiIiIbBYTHSIiIrJZTHSIiIjIZjHRISIiIpvFRIeIiIhsFhMdIiIisllMdIiIiMhmMdEhIiIim8VEh4iIiGwWEx0iIiKyWUx0iIiIyGYx0SEiIiKbxUSHiGzOmTNnoFKpsGXLFnTr1g0eHh5o3Lgxfv75Z6VDIyIzY6JDRDYnNTUVKpUK7733Hl599VWkpqaiXr16mD59utKhEZGZMdEhIpuTmpoKb29vrFmzBnFxcWjUqBEeeeQRXL9+XenQiMjMmOgQkc1JTU1FfHw8/Pz8iradO3cODRs2VDAqIlICEx0isjmpqalo37691rZjx46hZcuWygRERIphokNENiUrKwsXLlxAq1attLYfP36ciQ6RHWKiQ0Q2JTU1FY6OjmjRokXRtgsXLuDWrVtMdIjsEBMdIrIpqampiIyMhLu7e9G2Y8eOoVatWggPD1cuMCJShEoIIZQOgoiIiMgU2KJDRERENouJDhEREdksJjpERERks5joEBERkc1iokNEREQ2i4kOERER2SwmOkRERGSzmOgQERGRzWKiQ0RERDaLiQ4RERHZLCY6REREZLP+D1xceXoCGlR4AAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject Text(24.000000000000007, 0.5, '(time for b) / n×(time for c)')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "loglog(ns, t_ratio, \"ro-\")\n",
    "title(\"1(e): Computers are more complicated than you think\")\n",
    "xlabel(L\"n\")\n",
    "ylabel(\"(time for b) / n×(time for c)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef31bc89",
   "metadata": {},
   "source": [
    "The ratio (the speedup from part c) actually **grows with n**, mainly because part (b) also uses more memory ($\\sim n^2$ rather than $\\sim n$) and [using more memory is slower](https://en.wikipedia.org/wiki/Memory_hierarchy).\n",
    "\n",
    "Interestingly, part (b) is actually *faster* than (c) for very small matrices, possibly because matrix–vector products are so heavily optimized in [linear-algebra software](https://www.openblas.net/).   We could get a bit faster in (c) by optimizing `b - x * ((x'*b)/(x'*x))` in order to [fuse the loops](https://julialang.org/blog/2017/01/moredots/), or maybe by calling an optimized [BLAS](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms) operation for [axpy](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.axpy!)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfa41bf1",
   "metadata": {},
   "source": [
    "## Problem 2 (8+4 points)\n",
    "\n",
    "**(a)** Describe the four fundamental subspaces of the **rank-1 matrix** $A = uv^T$ where $u \\in \\mathbb{R}^m$ and $n \\in \\mathbb{R}^n$.\n",
    "\n",
    "**(b)** For any column vectors $u,v \\in \\mathbb{R}^3$, the matrix $uv^T$ is rank 1, except when \\_\\_\\_\\_\\_\\_\\_\\_, in which case $uv^T$ has rank \\_\\_\\_\\_."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54223d31",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "**(a)** Since the matrix is given to be \"rank-1\", we can assume $u, v \\ne 0$ (see part b).  In that case, the four subspaces are:\n",
    "\n",
    "1. Column space $C(uv^T) \\subseteq \\mathbb{R}^m$ is 1-dimensional and is **spanned by u**.\n",
    "2. Row space $C(vu^T) \\subseteq \\mathbb{R}^n$ is 1-dimensional and is **spanned by v**.\n",
    "3. Nullspace $N(uv^T)$ is $n-1$ dimensional and is **every vector orthogonal to v**.\n",
    "3. Left nullspace $N(vu^T)$ is $m-1$ dimensional and is **every vector orthogonal to u**.\n",
    "\n",
    "**(b)** $uv^T$ is always rank 1 exept when **u and/or v is 0**, in which case it has **rank 0**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82b2c169",
   "metadata": {},
   "source": [
    "## Problem 3 (5+4+4+4 points)\n",
    "\n",
    "**(a)** Pick the choices that makes this statement correct for arbitrary matrices A and B:  $C(AB)$ (*contains / is contained in*) the column space of (*A / B*).  Briefly justify your answer.\n",
    "\n",
    "**(b)** Suppose that $A$ is a $1000\\times 1000$ matrix of rank $< 10$.  Suppose we multiply it by 10 random vectors $x_1, x_2, \\ldots, x_{10}$, e.g. generated by `randn(1000)`.  How could we use the results to get a $10\\times 10$ matrix $C$ whose rank (almost certainly) matches $A$'s?\n",
    "\n",
    "**(c)** Suppose we instead make $1000\\times 10$ matrix $X$ whose columns are $x_1, x_2, \\ldots, x_{10}$.  Give a formula for the *same* matrix $C$ in terms of matrix products involving $A$ and $X$.\n",
    "\n",
    "**(d)** Fill in the code for $C$ below, and compare the biggest 10 singular values of $A$ (chosen to be rank ≈ 4 in this case) to the corresponding 10 singular values of $C$.  Does it match what you expect?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dc3782a",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "**(a)** C(AB) **is contained in** the column space of **A**.   The reason is that any $y \\in C(AB)$, we have $y = ABx = A(Bx)$, so $y \\in C(A)$.   (The converse is not necessarily true, e.g. if $B = 0$.)\n",
    "\n",
    "**(b and c)** I will answer these together.  The only reason for asking (c) separately was to make sure you know how to express your answer in matrix notation.\n",
    "\n",
    "First, realize that $AX$ *almost certainly* (probability 1) has the **same rank as A**, i.e. we have a $1000 \\times 10$ matrix with the same rank as A.\n",
    "* Proof: Suppose $C(AX)$ were lower-dimensional than $C(A)$, in which case consider the orthogonal complement of $C(AX)$ in $C(A)$: there must be some direction $y\\ne 0$ in $C(A) = N(A^T)^\\perp$ perpendicular to everything in $C(AX)$, i.e. $y^TAX = 0$.  But that would require every one of our \"random\" vectors $x_k$ to be orthogonal to $A^T y$ (≠ 0 since $y \\notin  N(A^T)$), which happens with probability 0.\n",
    "\n",
    "Now, how can we reduce this to a $10\\times 10$ matrix?  If we are to use only elementary matrix operations (no fair using something like an SVD that you don't personally know how to compute), this suggests multipling $AX$ on the left by some $10 \\times 1000$ matrix.  But which one?  There are several valid choices:\n",
    "\n",
    "1. You could use $X^T AX$.  (This could only reduce the rank if $X$ was $\\perp$ to some portion of $C(A)$, which happens with probability 0.)\n",
    "2. You could use $Y^T AX$, where $Y$ is another $1000 \\times 10$ random matrix like `randn(1000,10)`.  (Again almost surely will have the same rank.)\n",
    "3. You could use $(AX)^T (AX) = X^T A^T A X$.   In 18.06 and similar intro linear-algebra courses we prove that $B^T B$ has the same rank as $B$ for any $B$.  (Proof by showing that $N(B) = N(B^T B)$:  $Bx = 0 \\implies B^TBx = 0$, and $B^T Bx = 0 \\implies x^T B^T B x = 0 = \\Vert Bx \\Vert^2 \\implies Bx = 0$.)\n",
    "\n",
    "I think these are the easiest possibilities.\n",
    "\n",
    "**(d)** Let's try all three possibilities from above in Julia:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e7b32901",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(svdvals(A))[1:10] = [1069.3195877770895, 1032.123298970407, 998.5420983487572, 963.6168764929388, 1.197027719541634e-12, 9.397013502391205e-13, 8.816489090342173e-13, 8.368786422681453e-13, 7.163119064816855e-13, 6.558305069443912e-13]\n",
      "svdvals(X' * A * X) = [17355.334466001634, 10725.604519449384, 8360.035414898422, 2369.8158461371777, 7.68456896867752e-12, 5.9484402119558325e-12, 2.9567758084067988e-12, 2.0377094428213978e-12, 1.0361892803083682e-12, 1.7189758977956257e-13]\n",
      "svdvals(Y' * A * X) = [11903.189430493283, 11182.390116386241, 6413.062698992587, 3644.0437181289303, 4.420883286604323e-12, 3.670046234886016e-12, 3.053253554336238e-12, 2.327401156370041e-12, 2.2006586809078315e-12, 1.4777079446259352e-12]\n",
      "svdvals((A * X)' * A * X) = [1.8887124303351242e7, 1.2865209490083521e7, 7.204356096870981e6, 3.603491360462107e6, 4.0543542578357826e-9, 3.7280041982774055e-9, 2.130547526629226e-9, 1.2473174612561205e-9, 1.113501138973212e-9, 8.507451167986952e-12]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "10-element Vector{Float64}:\n",
       " 1.8887124303351242e7\n",
       " 1.2865209490083521e7\n",
       " 7.204356096870981e6\n",
       " 3.603491360462107e6\n",
       " 4.0543542578357826e-9\n",
       " 3.7280041982774055e-9\n",
       " 2.130547526629226e-9\n",
       " 1.2473174612561205e-9\n",
       " 1.113501138973212e-9\n",
       " 8.507451167986952e-12"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "# random 1000x1000 matrix of rank 4\n",
    "A = randn(1000, 4) * randn(4, 1000)\n",
    "@show svdvals(A)[1:10]\n",
    "\n",
    "# possibility 1: XᵀAX\n",
    "X = randn(1000, 10)\n",
    "@show svdvals(X'*A*X)\n",
    "\n",
    "# possibility 2: YᵀAX\n",
    "Y = randn(1000, 10)\n",
    "@show svdvals(Y'*A*X)\n",
    "\n",
    "# possibility 3: (AX)ᵀAX\n",
    "@show svdvals((A*X)'*A*X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "626f50e3",
   "metadata": {},
   "source": [
    "A has rank 4, so it should only have 4 nonzero singular values in theory.  In fact, *all* of its singular values are nonzero, but **only 4** of them (≈ 1069, 1032, 998, 963) are **non-negligible** — the remaining singular values are on the order of [machine precision 10⁻¹⁶](https://en.wikipedia.org/wiki/Machine_epsilon) times $\\sigma_1 \\approx 1000$, i.e. they are $\\sim 10^{-12}$.  These are due to [roundoff errors](https://en.wikipedia.org/wiki/Round-off_error) in computer arithmetic.\n",
    "\n",
    "For all 3 of our $10\\times 10$ matrices, we indeed see the same thing: **4 non-negligible** singular values, and the rest are roundoff errors.  Indeed, 4 singular values are much *larger* than those of $A$, which is a consequence of our lack of normalization of $X$ and $Y$.  The biggest are the singular values of $(AX)^T AX = X^T A^T A X$, which shouldn't be surprising — this *squares* the singular values of $AX$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee67e16b",
   "metadata": {},
   "source": [
    "## Problem 4 (4+5+5 points)\n",
    "\n",
    "The famous Hadamard matrices are filled with ±1 and have orthogonal columns (orthonormal if we divide $H_n$ by $1/\\sqrt{n}$).  The first few are:\n",
    "\\begin{align}\n",
    " H_1 &= \\begin{pmatrix} 1 \\end{pmatrix}, \\\\\n",
    "  H_2 &= \\begin{pmatrix}\n",
    "    1 &  1 \\\\\n",
    "    1 & -1\n",
    "  \\end{pmatrix}, \\\\\n",
    "  H_4 &= \\begin{pmatrix} H_2 & H_2 \\\\ H_2 & -H_2 \\end{pmatrix} = \\begin{pmatrix}\n",
    "    1 &  1 &  1 &  1\\\\\n",
    "    1 & -1 &  1 & -1\\\\\n",
    "    1 &  1 & -1 & -1\\\\\n",
    "    1 & -1 & -1 &  1\n",
    "  \\end{pmatrix} \\, .\n",
    "\\end{align}\n",
    "Notice that (for power-of-2 sizes), they are built up \"recursively\" out of smaller Hadamard matrices.  Multiplying a vector by a Hadamard matrix requires no multiplications at all, only additions/subtractions.\n",
    "\n",
    "**(a)** If you multiply $H_4 x$ for some $x \\in \\mathbb{R}^4$ by the normal \"rows-times-columns\" method (*without* exploiting any special patterns), exactly how many scalar additions/subtractions are required?\n",
    "\n",
    "**(b)** Let's break $x$ into two blocks: $x = \\begin{pmatrix} x_1 \\\\ x_2 \\end{pmatrix}$ for $x_1, x_2 \\in \\mathbb{R}^2$.  Write out $H_4 x$ in terms of a sequence of $2\\times 2$ block multiplications with $\\pm H_2$.  You'll notice that some of these $2\\times 2$ multiplications are repeated.  If we re-use these repeated multiplications rather than doing them twice, we can save a bunch of arithmetic — what is the new count of scalar additions/subtractions if you do this?\n",
    "\n",
    "**(c)** Similarly, the $8\\times 8$ Hadamard matrix $H_8 = \\begin{pmatrix} H_4 & H_4 \\\\ H_4 & -H_4 \\end{pmatrix}$ is made out of $H_4$ matrices.   To multiply it by a vector $y \\in \\mathbb{R}^8$, the naive rows-times columns method would require \\_\\_\\_\\_ scalar additions/subtractions, whereas if you broke them up first into blocks of 4, used your solution from (b), and then re-used any repeated $H_4$ products, it would only require \\_\\_\\_\\_ scalar additions/subtractions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e419713b",
   "metadata": {},
   "source": [
    "### Solutions:\n",
    "\n",
    "**(a)** $H_4 x$ requires 4 dot products (with the rows of $H_4$), each of which requires 3 scalar additions or subtractions, so there are **12 additions/subtractions** overall.\n",
    "\n",
    "**(b)** We have:\n",
    "$$\n",
    "H_4 x = \\begin{pmatrix} H_2 x_1 + H_2 x_2 \\\\ H_2 x_1 - H_2 x_2 \\end{pmatrix}\n",
    "$$\n",
    "so we see that we need to compute $H_2 x_1$ and $H_2 x_2$ ($=2 \\times 2 \\times 1 = 4$ additions/subtractions) and then add/subtract the resulting 2-component vectors ($= 2\\times 2 = 4$ additions/subtractions), for **8 additions/subtractions overall**.\n",
    "\n",
    "**(c)** In the same way, a naive $H_8 y$ would require $8 \\times 7 = \\boxed{56}$ additions/subtractions, but breaking it into blocks would require 2 multiplications by $H_4$ (= 8 adds/subs each) and 2 additions/subtractions of 4-component vectors (= 8 add/subs), for $\\boxed{8 \\times 3 = 24}$ additions/subtractions.\n",
    "\n",
    "If you continue this process recursively, you obtain the [fast Walsh–Hadamard transform](https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform), which reduces the cost to multiply by $H_n$ from $\\Theta(n^2)$ to $\\Theta(n \\log n)$ for $n = 2^k$!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2a2b8eb",
   "metadata": {},
   "source": [
    "## Problem 5 (5+5 points)\n",
    "\n",
    "The famous \"discrete Fourier transform\" matrix $F$ has columns that are actually eigenvectors of the (unitary) permutation matrix:\n",
    "$$\n",
    "P = \\begin{pmatrix}  & 1 & & \\\\ & &1 &  \\\\ & &  &1 \\\\ 1 &  & & \\end{pmatrix}\n",
    "$$\n",
    "for the $4\\times 4$ case, and similarly for larger matrices.\n",
    "\n",
    "**(a)** One way of saying way Fourier transforms are practically important is that they *diagonalize* (are eigenvectors of) matrices that *commute* with $P$.  If $A$ is a $4\\times 4$ matrix whose first row is $(a\\, b\\, c\\, d)$\n",
    "$$\n",
    "A = \\begin{pmatrix} a & b & c & d \\\\ ? & ? &? & ? \\\\ ? &? & ? &? \\\\ ? & ? & ? &? \\end{pmatrix}\n",
    "$$\n",
    "that commutes with $P$ (i.e. $AP=PA$), what must be true of the other (\"?\") entries of $A$?\n",
    "\n",
    "**(b)** Fill in the matrix `A` in Julia below and fill in and run the code to check that it commutes with $P$ and is diagonalized by $F$:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14b558f8",
   "metadata": {},
   "source": [
    "### Solution:\n",
    "\n",
    "**(a)** By inspection:\n",
    "$$\n",
    "AP = \\begin{pmatrix} d & a & b & c \\\\ ? & ? &? & ? \\\\ ? &? & ? &? \\\\ ? & ? & ? &? \\end{pmatrix},\\; PA = \\begin{pmatrix} ? & ? &? & ? \\\\ ? &? & ? &? \\\\ ? & ? & ? &? \\\\ a & b & c & d\\end{pmatrix}\n",
    "$$\n",
    "i.e. $AP$ shifts the columns, and $PA$ shifts the rows.  For these two to be equal, the second row of $A$ must be \"d a b c\".  Continuing in this process, the third row of $A$ must be \"c d a b\", and the fourth row must be \"b c d a\", i.e. A must be of the form\n",
    "$$\n",
    "\\boxed{A = \\begin{pmatrix} a & b & c & d \\\\ d & a &b & c \\\\ c &d & a &b \\\\ b & c & d &a \\end{pmatrix}}\\, ,\n",
    "$$\n",
    "which is known as a [circulant matrix](https://en.wikipedia.org/wiki/Circulant_matrix) — it expresses a linear operation known as a [circular discrete convolution](https://en.wikipedia.org/wiki/Convolution#Circular_discrete_convolution).\n",
    "\n",
    "**(b)** See below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e7222137",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Int64}:\n",
       " 1  7  3  2\n",
       " 2  1  7  3\n",
       " 3  2  1  7\n",
       " 7  3  2  1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a, b, c, d = 1, 7, 3, 2    # 4 arbitrarily chosen values\n",
    "\n",
    "P = [0 1 0 0\n",
    "     0 0 1 0\n",
    "     0 0 0 1\n",
    "     1 0 0 0]\n",
    "\n",
    "F = im .^ ((0:3) .* (0:3)') # the 4×4 Fourier matrix\n",
    "\n",
    "# fill in:\n",
    "A = [a b c d\n",
    "     d a b c\n",
    "     c d a b\n",
    "     b c d a]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "25dc7673",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check:\n",
    "P * A == A * P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9344787",
   "metadata": {},
   "source": [
    "Yay, it commutes!\n",
    "\n",
    "Now we want to check that $F$ \"diagonalizes\" $A$, which was defined (in the problem) as (columns of) $F$ being eigenvectors of $A$.\n",
    "\n",
    "There are several ways to check this, the easiest being to check that $A$ is diagonal in the **basis of F**, i.e. that $F^{-1} A F$ is diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0ae85697",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{ComplexF64}:\n",
       " 13.0+0.0im   0.0+0.0im   0.0+0.0im   0.0+0.0im\n",
       "  0.0+0.0im  -2.0+5.0im   0.0+0.0im   0.0+0.0im\n",
       "  0.0+0.0im   0.0+0.0im  -5.0+0.0im   0.0+0.0im\n",
       "  0.0+0.0im   0.0+0.0im   0.0+0.0im  -2.0-5.0im"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A is diagonal in the basis of F:\n",
    "\n",
    "F^-1 * A * F   # or F \\ A * F, or inv(F) * A * F, or F'/4 * A * F"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "897b25b0",
   "metadata": {},
   "source": [
    "(From class, note also that $F^{-1} = \\overline{F^T} / 4$: F is unitary up to a scale factor.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e6492277",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(F) ≈ F'/4   # note that F' is the conjugate-transpose"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "714942d0",
   "metadata": {},
   "source": [
    "We could also, of course, compute eigenvectors of $A$ and compare them to F, but this is a lot more work (both for Julia and for us).   Julia will give us eigenvectors, but it might give them in a different order and will give them with a different normalization (normalized to unit length), and might give them with different phases or signs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "dce982d1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{ComplexF64}:\n",
       "  0.5+0.0im  -3.67437e-17+0.5im          …  -0.5+0.0im\n",
       " -0.5+0.0im           0.5-0.0im             -0.5+0.0im\n",
       "  0.5+0.0im   1.27272e-16-0.5im             -0.5+0.0im\n",
       " -0.5+0.0im          -0.5-4.19352e-16im     -0.5+0.0im"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "X = eigvecs(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acebc0b5",
   "metadata": {},
   "source": [
    "Also, it doesn't show the whole matrix to save space.  To make it look more like $F$, we can divide the columns of $X$ so that their first row is \"1 1 1 1\", and then call a lower-level `show` function to make sure we show the whole matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c81825b9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4×4 Matrix{ComplexF64}:\n",
      "  1.0+0.0im           1.0-0.0im                   1.0+0.0im          1.0-0.0im\n",
      " -1.0+0.0im  -7.34873e-17-1.0im          -7.34873e-17+1.0im          1.0-0.0im\n",
      "  1.0+0.0im          -1.0-1.81057e-16im          -1.0+1.81057e-16im  1.0-0.0im\n",
      " -1.0+0.0im  -7.65216e-16+1.0im          -7.65216e-16-1.0im          1.0-0.0im"
     ]
    }
   ],
   "source": [
    "show(stdout, \"text/plain\", X ./ transpose(X[1,:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72946b10",
   "metadata": {},
   "source": [
    "If you stare at this long enough, you'll see that they are indeed the columns of $F$, but in a permuted order:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "cb1467ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X ./ transpose(X[1,:]) ≈ F[:, [3,4,2,1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd09f649",
   "metadata": {},
   "source": [
    "But that is a lot of effort!   It's a bit silly to check whether something is an eigenvector by anything other than multiplying it by $A$.   (Analogously, if you wanted to check whether $x$ solves $Ax=b$, you wouldn't re-do Gaussian elimination and then compare the result to $x$, you would just check whether $Ax \\approx b$.)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}