jstac/pcr.ipynb

## pcr.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "f1d2cabd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'%.3f'"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%precision 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "e76d30a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.matrices import Matrix\n",
    "import numpy as np\n",
    "from numpy.random import rand, randn\n",
    "from numpy.linalg import svd, inv\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77707729",
   "metadata": {},
   "source": [
    "The model is \n",
    "\n",
    "$$ y = \\beta_0 x_0 + \\beta_1 x_1 + \\beta_2 x_2 + \\epsilon $$\n",
    "\n",
    "The column vectors $x_0$ and $x_1$ are independent but $x_2 = x_1 + \\eta$ where $\\eta$ has low variance.\n",
    "\n",
    "The error vector obeys $\\epsilon \\sim N(0, \\sigma^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 372,
   "id": "b5914a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 16\n",
    "k = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 373,
   "id": "2d738776",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.empty((n, k))\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 374,
   "id": "80418a7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "X[:,0] = rand(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 375,
   "id": "de4a1cba",
   "metadata": {},
   "outputs": [],
   "source": [
    "X[:,1] = rand(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 376,
   "id": "4b607ed8",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.15  # small value means x_2 \\approx x_1\n",
    "η = a * randn(n)\n",
    "X[:,2] = X[:,1] + η"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 377,
   "id": "2c05ba4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.696,  0.182,  0.144],\n",
       "       [ 0.286,  0.175, -0.244],\n",
       "       [ 0.227,  0.532,  0.266],\n",
       "       [ 0.551,  0.532,  0.427],\n",
       "       [ 0.719,  0.634,  0.774],\n",
       "       [ 0.423,  0.849,  0.823],\n",
       "       [ 0.981,  0.724,  0.725],\n",
       "       [ 0.685,  0.611,  0.714],\n",
       "       [ 0.481,  0.722,  0.591],\n",
       "       [ 0.392,  0.323,  0.366],\n",
       "       [ 0.343,  0.362,  0.241],\n",
       "       [ 0.729,  0.228, -0.031],\n",
       "       [ 0.439,  0.294,  0.235],\n",
       "       [ 0.06 ,  0.631,  0.717],\n",
       "       [ 0.398,  0.092,  0.143],\n",
       "       [ 0.738,  0.434,  0.432]])"
      ]
     },
     "execution_count": 377,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "id": "b885c366",
   "metadata": {},
   "outputs": [],
   "source": [
    "β = -1, 1, 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "id": "552f8927",
   "metadata": {},
   "outputs": [],
   "source": [
    "σ = 1.2\n",
    "y = X @ β + σ * randn(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 380,
   "id": "206a0072",
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, Vt = np.linalg.svd(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "id": "be1f0c47",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.433, 1.064, 0.411])"
      ]
     },
     "execution_count": 392,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 393,
   "id": "9fae9640",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAja0lEQVR4nO3deXRUdZ738fe3lpBAAoEkCIQdArLJYoC44dYe1NHW9nFaUVxZRtve7Bm6x35OO2OfcbpPM2132446bOIGbR/Htnt8dNSxG3CZRAKyI2GVHcISIBAgqfo9f1SZjiFABarqJlWf1zl1rKp7b9XH8vrh8qvfvWXOOUREpPXzeR1ARETiQ4UuIpIiVOgiIilChS4ikiJU6CIiKSLg1Rvn5+e73r17e/X2IiKt0pIlS/Y55wqaWuZZoffu3Zvy8nKv3l5EpFUysy9Ot0xDLiIiKeKshW5mmWb2qZktN7PVZvZEE+tcZWaHzGxZ9PZ4YuKKiMjpxDLkcgK4xjlXbWZB4CMze8c5V9povQ+dczfFP6KIiMTirIXuItcGqI4+DEZvul6AiEgLE9MYupn5zWwZsBd43zlX1sRql0SHZd4xsyGneZ2pZlZuZuWVlZXnnlpERE4R0ywX51wIGGFmucAfzGyoc25Vg1WWAr2iwzI3Am8CRU28zgxgBkBxcXGzj/Lf/GwH099dx86qGrrlZjFt/EBuHVnY3JcREUlJzZrl4pyrAhYA1zd6/rBzrjp6/20gaGb5ccoIRMr8sTdWsqOqBgfsqKrhsTdW8uZnO+L5NiIirVYss1wKokfmmFkW8DXg80brdDEzi94fE33d/fEMOv3dddTUhr7yXE1tiOnvrovn24iItFqxDLl0BV40Mz+Rov69c+4tM3sIwDn3PHA78LCZ1QE1wJ0uzhda31lV06znRUTSTSyzXFYAI5t4/vkG958BnolvtK/qlpvFjibKu2tuZiLfVkSk1Wg1Z4pOGz+QrKD/lOfzszOoDYU9SCQi0rK0mkK/dWQhP7ttGIW5WRhQmJvFrSO6sWL7YR59bRl1KnURSXOeXZzrXNw6svCUaYqDurbnZ+98TsBn/PKbI/D7zKN0IiLealWF3pS/u7IfdWHH9HfX4ff5mH77RfhU6iKShlp9oQM8cnV/6kKOX/1PBQGf8bPbhqnURSTtpEShA3zva0XUhcP89s8bCPiNf7l1KNGp8SIiaSFlCh3gB9cNoC7seG7BRgI+45+/PkSlLiJpI6UK3cz44fiB1IXCzPxwM36fj5/cNEilLiJpIaUKHSKl/uMbB1EXdsz5eDMBv/HYDReq1EUk5aVcoUOk1B+/aTB1IceMRZsI+Ixp4weq1EUkpaVkoUOk1J/4+hDqwo5nF2wk4Pfxg+sGeB1LRCRhUrbQAXw+48lbhxIKh3n6g/UEfMZ3rz3lMu0iIikhpQsdIqX+89suoi7seOr9CgJ+41tX9fc6lohI3KV8oUOk1KffPpxQ2PGL/15H0Odjyri+XscSEYmrtCh0AL/P+OXfDqcu7Hjy7bX4fcaDl/fxOpaISNykTaEDBPw+fn3HCMJhx0/fWkPAb9x7SW+vY4mIxEWruXxuvAT9Pn5z50i+NugCHv/jauaVbfU6kohIXKRdoQNkBHz8+90juebCzvz4Dyv5/eJtXkcSETlvaVnoAG0Cfp69exTjBhTwozdW8PqS7V5HEhE5L2lb6ACZQT8z7rmYy/rlM+315fxx2Q6vI4mInLO0LnSIlPrMe4sZ26cTj762jP9avtPrSCIi5yTtCx0gK8PPnPtHU9yrE99/bRnvrNzldSQRkWZToUe1zQgw54HRjOiRy3fmf8Z7q3d7HUlEpFnOWuhmlmlmn5rZcjNbbWZPNLGOmdnTZrbBzFaY2ajExE2s7DYB5j4wmqGFHXhk3lI+WLvH60giIjGL5Qj9BHCNc244MAK43sxKGq1zA1AUvU0FnotnyGTKyQzy4oNjGNS1PQ+/spQF6/Z6HUlEJCZnLXQXUR19GIzeXKPVbgFeiq5bCuSaWdf4Rk2eDllBXnpwDP07ZzP15SV8tH6f15FERM4qpjF0M/Ob2TJgL/C+c66s0SqFQMOzc7ZHn2u1cttm8OrksfTNb8fklxbzyUaVuoi0bDEVunMu5JwbAXQHxpjZ0EarNPVTQI2P4jGzqWZWbmbllZWVzQ6bbB3bRUq9Z6e2TJpbTtmm/V5HEhE5rWbNcnHOVQELgOsbLdoO9GjwuDtwyoRu59wM51yxc664oKCgeUk9kpfdhlcnl9AtN5MH5i5myRcHvI4kItKkWGa5FJhZbvR+FvA14PNGq/0JuDc626UEOOScS5nJ3AU5bZg/pYQu7TO5b85iPtt60OtIIiKniOUIvSvwFzNbASwmMob+lpk9ZGYPRdd5G9gEbABmAt9KSFoPdW6fybwpJeRlZ3Dv7E9Zsb3K60giIl9hzp0y1J0UxcXFrry83JP3Ph87q2q4Y8b/cuhYLfOmlDC0sIPXkUQkjZjZEudccVPLdKZoM3XLzWLe5BJyMoNMnF3Gmp2HvY4kIgKo0M9Jj05tmT+lhKygn4mzy1i3+4jXkUREVOjnqmdeW+ZNKSHoN+6eVcqGvSp1EfGWCv089Mlvx7wpJZgZE2aWsbGy+uwbiYgkiAr9PPUryGbe5LGEw467ZpayZd9RryOJSJpSocdB0QU5zJtSQm3IMWFmKVv3H/M6koikIRV6nAzsksMrk8ZSUxtiwsxSth9UqYtIcqnQ42hwt/a8MmksR47XMmFmKTuraryOJCJpRIUeZ0MLO/DypLFUHa3lrpml7D503OtIIpImVOgJMLxHLi9OGsO+6pPcNbOUvYdV6iKSeCr0BBnVsyNzHxjN7sPHuWtWGZVHTngdSURSnAo9gYp7d+KF+0ez42ANE2eVsb9apS4iiaNCT7CxffOYfV8xW/Yf5e5ZZRw8etLrSCKSolToSXBp/3xm3VfMpn1HmTi7jEPHar2OJCIpSIWeJFcUFfAf91zM+j3V3DOnjEM1KnURiS8VehJdPbAzz00cxdpdh7lvzqccOa5SF5H4UaEn2bWDLuCZu0axaschHnhhMUdP1HkdSURShArdA+OHdOG3E0by2bYqHpi7mGMnVeoicv5U6B65YVhXfn3HCMq3HGDS3HJqToa8jiQirZwK3UM3D+/GU98cQenm/Ux9uZzjtSp1ETl3KnSP3TqykOm3D+ejDfv4u5eXcKJOpS4i50aF3gLcfnF3fn7bMBZWVPLwK0s5WRf2OpKItEIq9BbijtE9efIbQ/nz53v59ryl1IZU6iLSPCr0FuTusb346S1DeG/NHr47/zOVuog0y1kL3cx6mNlfzGytma02s+81sc5VZnbIzJZFb48nJm7qu/eS3vzkpsG8s2o3j762jDqVuojEKBDDOnXA3zvnlppZDrDEzN53zq1ptN6Hzrmb4h8x/Uy6vA+hcJh/fftzAj7jl98cgd9nXscSkRburIXunNsF7IreP2Jma4FCoHGhSxxNHdeP2pBj+rvr8Pt8TL/9InwqdRE5g1iO0OuZWW9gJFDWxOJLzGw5sBP4B+fc6ia2nwpMBejZs2ezw6abR67uT13I8av/qSDgM3522zCVuoicVsyFbmbZwH8C33fOHW60eCnQyzlXbWY3Am8CRY1fwzk3A5gBUFxc7M41dDr53teKCIXDPP3nDfj9xpO3DsVMpS4ip4qp0M0sSKTMX3XOvdF4ecOCd869bWbPmlm+c25f/KKmr0evG0Bt2PHcgo0EfMYTXx+iUheRU5y10C3SHLOBtc65p06zThdgj3POmdkYIrNn9sc1aRozM344fiB1oTAzP9xMwOfjJzcNUqmLyFfEcoR+GXAPsNLMlkWf+zHQE8A59zxwO/CwmdUBNcCdzjkNqcSRmfHjGwdRF3bM+XgzAb/x2A0XqtRFpF4ss1w+As7YGs65Z4Bn4hVKmmZmPH7TYEJhx4xFmwj4jGnjB6rURQRo5iwX8Z6Z8c83D6E25Hh2wUYCfh8/uG6A17FEpAVQobdCPl9ktksoHObpD9YT8BnfvfaUSUUikmZU6K2Uz2f8/LaLqAs7nnq/goDf+NZV/b2OJSIeUqG3Yj6fMf324YTCjl/89zqCPh9TxvX1OpaIeESF3sr5fcYv/zZS6k++vRa/z3jw8j5exxIRD6jQU0DA7+NXd4wgFHb89K01BPzGvZf09jqWiCSZroeeIoJ+H7+5cyTXDb6Ax/+4mnllW72OJCJJpkJPIRkBH8/cNZJrLuzMj/+wkt8v3uZ1JBFJIhV6imkT8PPs3aMYN6CAH72xgteXbPc6kogkiQo9BWUG/cy452Iu65fPtNeX88dlO7yOJCJJoEJPUZlBPzPvLaakTx6PvraM/1q+0+tIIpJgKvQUlpXhZ/b9xRT36sT3X1vGOyt3eR1JRBJIhZ7i2mYEmPPAaEb0yOU78z/jvdW7vY4kIgmiQk8D2W0CzH1gNEMLO/DIvKV8sHaP15FEJAFU6GkiJzPIiw+OYVDX9jz8ylIWrNvrdSQRiTMVehrpkBXkpQfHUHRBNlNfXsJH6/ULgSKpRIWeZnLbZvDKpLH0zW/HpBcX88lGlbpIqlChp6GO7TJ4dfJYeuW1ZdLccso26edfRVKBCj1N5WW34dXJJXTLzeSBuYsp33LA60gicp5U6GmsIKcN86eU0KV9Jve/sJilWw96HUlEzoMKPc11bp/JvCkl5GVncN/sT1mxvcrrSCJyjlToQpcOmcyfUkJuuyATZ5WxaschryOJyDlQoQsA3XKzmDe5hJzMIBNnl7Fm52GvI4lIM6nQpV6PTm2ZP6WErKCfibPLWLf7iNeRRKQZzlroZtbDzP5iZmvNbLWZfa+JdczMnjazDWa2wsxGJSauJFrPvEipB/3G3bNK2bBXpS7SWsRyhF4H/L1zbhBQAjxiZoMbrXMDUBS9TQWei2tKSare+e2YN6UEM2PCzDI2VlZ7HUlEYnDWQnfO7XLOLY3ePwKsBQobrXYL8JKLKAVyzaxr3NNK0vQryGbe5LE457hrZilb9h31OpKInEWzxtDNrDcwEihrtKgQaPgDlts5tfQxs6lmVm5m5ZWVlc2MKslWdEEOr04uoTbkmDCzlK37j3kdSUTOIOZCN7Ns4D+B7zvnGk+BsCY2cac84dwM51yxc664oKCgeUnFEwO75PDKpLHU1IaYMLOU7QdV6iItVUyFbmZBImX+qnPujSZW2Q70aPC4O6DfPEsRg7u155VJYzlyvJYJM0vZWVXjdSQRaUIss1wMmA2sdc49dZrV/gTcG53tUgIccs7p985SyNDCDrw8aSxVR2u5a2Ypuw8d9zqSiDQSyxH6ZcA9wDVmtix6u9HMHjKzh6LrvA1sAjYAM4FvJSaueGl4j1xenDSGfdUnuWtmKXsPq9RFWhJz7pSh7qQoLi525eXlnry3nJ/yLQe4d86ndMvNYv6UEgpy2ngdSSRtmNkS51xxU8t0pqg0W3HvTrxw/2h2HKxh4qwy9lef8DqSiKBCl3M0tm8es+8rZsv+o9w9q4yDR096HUkk7anQ5Zxd2j+fWfcVs2nfUSbOLuPQsVqvI4mkNRW6nJcrigqYcc/FrN9TzT1zyjhUo1IX8YoKXc7bVQM789zEUazddZj75nzKkeMqdREvqNAlLq4ddAHP3DWKVTsOcf8Li6k+Ued1JJG0o0KXuBk/pAu/nTCSZduqePCFxRw7qVIXSSYVusTVDcO68us7RlD+xQEmzS2n5mTI60giaUOFLnF38/BuPPXNEZRu3s+Ul8o5XqtSF0kGFbokxK0jC5l++3A+3riPv3t5iUpdJAlU6JIwt1/cnZ/fNoyFFZV869WlnKwLex1JJKWp0CWh7hjdkye/MZQ/f76XR+YtpTakUhdJFBW6JNzdY3vx01uG8P6aPXx3/mcqdZEEUaFLUtx7SW9+ctNg3lm1m0dfW0adSl0k7gJeB5D0MenyPoTCYf717c8J+IxffnMEfl9Tv14oIudChS5JNXVcP2pDjunvrsPv8zH99ovwqdRF4kKFLkn3yNX9CYUdT71fQcBn/Oy2YSp1kThQoYsnvnttEXWhME//eQN+v/HkrUOJ/HytiJwrFbp45tHrBlAbdjy3YCMBn/HE14eo1EXOgwpdPGNm/HD8QEJhx4xFmwj4fPzkpkEqdZFzpEIXT5kZj91wIbWhMHM+3kzAH3msUhdpPhW6eM7MePymwQ2O1I1p4weq1EWaSYUuLYKZ8c83D6Eu7Hh2wUYCfh8/uG6A17FEWhUVurQYPp/xL7cMJRRyPP3BegI+47vXFnkdS6TVOOup/2Y2x8z2mtmq0yy/yswOmdmy6O3x+MeUdOGLzkv/P6O689T7FTy7YIPXkURajViO0OcCzwAvnWGdD51zN8UlkaQ9n8/4xe0XEQqH+cV/ryPo8zFlXF+vY4m0eGctdOfcIjPrnYQsIvX8PuPf/nY4dWHHk2+vxe8zHry8j9exRFq0eI2hX2Jmy4GdwD8451Y3tZKZTQWmAvTs2TNOby2pKuD38as7RhAKO3761hoCfuPeS3p7HUukxYrH5XOXAr2cc8OB3wJvnm5F59wM51yxc664oKAgDm8tqS7o9/H0hJFcN/gCHv/jauaVbfU6kkiLdd6F7pw77Jyrjt5/GwiaWf55JxOJCvp9PHPXSK65sDM//sNKfr94m9eRRFqk8y50M+ti0TNAzGxM9DX3n+/rijTUJuDn2btHceWAAn70xgpeX7Ld60giLc5Zx9DNbD5wFZBvZtuBfwKCAM6554HbgYfNrA6oAe50zrmEJZa0lRn08x/3XMzkF8uZ9vpyAj7j1pGFXscSaTHMq+4tLi525eXlnry3tG41J0M8OHcxZZv385s7R3Lz8G5eRxJJGjNb4pwrbmqZflNUWp2sDD+z7y+muFcnvv/aMt5ZucvrSCItggpdWqW2GQHmPDCaET1y+c78z3hv9W6vI4l4ToUurVZ2mwBzHxjN0MIOPDJvKR+s3eN1JBFPqdClVcvJDPLig2MY1LU9D7+ylAXr9nodScQzKnRp9TpkBXn5wbEUXZDN1JeX8NH6fV5HEvGECl1SQoe2QV6ZNJa++e2Y9OJiPtmoUpf0o0KXlNGxXQavTh5Lr7y2TJpbTtkmnd8m6UWFLiklL7sNr04uoVtuJg/MXUz5lgNeRxJJGhW6pJyCnDbMn1JCl/aZ3P/CYpZuPeh1JJGkUKFLSurcPpN5U0rIy87gvtmfsmJ7ldeRRBJOhS4pq0uHTOZPKSG3XZCJs8pYteOQ15FEEkqFLimtW24W86eUkJMZZOLsMtbsPOx1JJGEUaFLyuvesS3zp5SQFfQzcXYZ63Yf8TqSSEKo0CUt9MyLlHrQb9w9q5QNe1XqknpU6JI2eue3Y96UEsyMCTPL2FhZ7XUkkbhSoUta6VeQzfwpY3HOcdfMUrbsO+p1JJG4UaFL2unfOYdXJ5dQG3JMmFnK1v3HvI4kEhcqdElLA7vk8MqksdTUhpgws5TtB1Xq0vqp0CVtDe7WnlcmjeXI8VomzCxlZ1WN15FEzot+U1TS3vJtVUycVUZedgb3X9abmYs2s7Oqhm65WUwbP1A/RC0tin5TVOQMhvfI5cVJY9hVVcMTf1rDjqoaHLCjqobH3ljJm5/t8DqiSExU6CLAqJ4dyckK0vjvqzW1Iaa/u86TTCLNpUIXidpffbLJ53dU1bB+zxG8Gp4UiVXgbCuY2RzgJmCvc25oE8sN+A1wI3AMuN85tzTeQUUSrVtuFjtO88Xodb9aRLcOmYwbUMCVAwq4tH8+HbKCSU4ocmZnLXRgLvAM8NJplt8AFEVvY4Hnov8UaVWmjR/IY2+spKY2VP9cVtDPD68fSGbQz8J1lfy/Fbv43eJt+H3GqJ65XDmggHEDChjarQM+n3mYXiSGQnfOLTKz3mdY5RbgJRf5+2ipmeWaWVfn3K54hRRJhi9ns0x/d12Ts1wmjOlJbSjMsm1VLFxXycKKSv7tvQr+7b0K8tplcEVRPlcOLOCKogLys9t4+a8iaSqmaYvRQn/rNEMubwE/d859FH38AfAj59wpcxLNbCowFaBnz54Xf/HFF+eXXsRj+6pP8OH6ShZV7GNRRSX7j0bG4YcWtufKAQVcOaAzI3vmEvTr6yqJjzNNW4xlyOWsr9/Ec03+KeGcmwHMgMg89Di8t4in8rPb8I2R3fnGyO6Ew47VOw+zsGIviyr28fzCTfz7XzaS0ybApf3zuHJAZ8YNyKd7x7Zex5YUFY9C3w70aPC4O7AzDq8r0qr4fMaw7h0Y1r0D376miMPHa/lkwz4WVlSycF0l767eA0D/ztnRo/cCxvTpRGbQ73FySRXxKPQ/Ad82s98R+TL0kMbPRaB9ZpDrh3bl+qFdcc6xsbKaBdGx95dLv2D2R5tpE/BR0jcvUvADC+ib347IxDGR5jvrGLqZzQeuAvKBPcA/AUEA59zz0WmLzwDXE5m2+EBT4+eN6dR/SWc1J0OUbt7PoopIwW+qjFzGtzA3iysHRqdG9ssjJ1NTI+WrzjSGrmu5iLQA2w4cY2FFJYsqKvl4wz6OngwR8BkX9+pYP/d9cNf2mhopKnSR1uRkXZilWw/Wj72v2RX5Yev87DaMG5DPlQMiUyM7tcvwOKl4QYUu0ortPXKcDysiX65+uL6Sg8dqMYOLCjvUj70P755LQFMj04IKXSRFhMKOVTsORY7eKyr5bOtBwg7aZwa4vCi//szVrh2yvI4qCaJCF0lRh47V8tGGffVfru4+fByAgRfkRIdnOjO6T0faBDQ1MlWo0EXSgHOOij3V9Sc2fbr5ACdDYbKCfi7pl8e4onyuHNiZ3nltNTWyFVOhi6ShYyfrKN20v/66M1uiP4bds1Pb+hObLumXR7s28TgdRZJFhS4ifLH/aP3QzCcb93PsZIig3yju1al+7vuFXXJ09N7CqdBF5CtO1IVYsuUgC9dHpkZ+vvsIAJ1z2tR/sXpFUT65bTU1sqVRoYvIGe05fLz+xKYP1+/jUE0tPov83uqXBT+8ey5+ndjkORW6iMQsFHYs3/7Xa74v316Fc5DbNsjl/fPrx987t8/0OmpaUqGLyDk7ePQkH3151ciKSiqPnADgwi459WPvxb06kRHQiU3JoEIXkbhwzvH57iP1lyUo/+IAtSFH2ww/l/bLq/9Rj555uuZ7oqjQRSQhqk/U8b8bI1eNXFCxl20HIj+y3Se/XXTsPZ+Svnm0zdDUyHhRoYtIwjnn2LL/GAvX7WXR+n18snEfx2vDZPh9jOnTqf7L1QEXZGtq5HlQoYtI0h2vDVG+5SALK/aysKKSij3VAHRpn1l/UbHL+ufTIUvXfG8OFbqIeG7XoZr6E5s+XL+PI8fr8PuMkT1y66/5Pqywg675fhYqdBFpUepCYZZtq6ov+BU7DuEcdGqXwRVFf73me0FOG6+jtjgqdBFp0fZXn4hMjVxXyaL1leyrPgnAkG7t6+e9j+rVkaCu+a5CF5HWIxx2rNl1uH7e+9IvDlIXdmS3CUSmRg4sYFxRAT06pefUSBW6iLRaR47X8snG/fVz33dURaZG9itoVz/2XtI3j8xgelzzXYUuIinBOcfGyr9eNbJ0035O1IVpE/Axtu+XJzbl068gdadGqtBFJCUdrw1RtvlA/dj7hr2RqZGFuVn1R++X9s+jfWbqTI1UoYtIWth+8BiLKvaxsGIvH2/YT/WJyNTIi3t2rL/uzOCu7Vv11EgVuoikndpQmM+2VtWf2LRqx2EA8rMzuKKoIDo1Mp+87NY1NfK8C93Mrgd+A/iBWc65nzdafhXwR2Bz9Kk3nHM/PdNrqtBFJJkqj5zgow2V0eGZfRw4ehIzGFbYgXFFkTNXR/bIJdDCp0aeV6GbmR+oAK4DtgOLgQnOuTUN1rkK+Afn3E2xhlKhi4hXwmHHqp2H6sfel26tIhR25GQG6q/5Pm5AAd1ys7yOeoozFXosl0AbA2xwzm2KvtjvgFuANWfcSkSkhfL5jIu653JR91y+c20Rh2pq+aTBNd/fWbUbgKLO2fXXnRndu1OLnxoZS6EXAtsaPN4OjG1ivUvMbDmwk8jR+urGK5jZVGAqQM+ePZufVkQkATpkBblhWFduGNYV5xwb9lbXl/tLpV8w66PNZAZ9XNI3r372TJ/8di1uamQshd5U4sbjNEuBXs65ajO7EXgTKDplI+dmADMgMuTSvKgiIolnZhRdkEPRBTlMvqIvNSdDlG7eHxmeqajkif+KDE706JQVGZopKuDS/vlkt/H+mu+xJNgO9GjwuDuRo/B6zrnDDe6/bWbPmlm+c25ffGKKiHgjK8PP1QM7c/XAzgBsO3Cs/uj9D0t38ErpVoJ+4+JeHblyQGeuHFDAoK45nhy9x/KlaIDIl6LXAjuIfCl6V8MhFTPrAuxxzjkzGwO8TuSI/bQvri9FRaS1O1kXZskXB1lYETl6X7MrcmxbkNOmfubMFf3z6dguA4A3P9vB9HfXsbOqhm65WUwbP5BbRxY26z3jMW3xRuDXRKYtznHOPWlmDwE45543s28DDwN1QA3wA+fcJ2d6TRW6iKSavYePs2j9vug13yupOlaLGVzUPZcuOW1YUFHJibpw/fpZQT8/u21Ys0pdJxaJiCRZKOxYuSMyNXJhxV6Wbq1qcr3C3Cw+/sdrYn7dMxV6y55BLyLSSvl9xogeuXzva0W88a3LmpxdArAzevXIeFChi4gkwelOUornyUsqdBGRJJg2fiBZjU5Mygr6mTZ+YNzew/uJkyIiaeDLLz7Pd5bLmajQRUSS5NaRhXEt8MY05CIikiJU6CIiKUKFLiKSIlToIiIpQoUuIpIiPDv138wqgS/OcfN8oCVeybGl5oKWm025mke5micVc/VyzhU0tcCzQj8fZlZ+umsZeKml5oKWm025mke5mifdcmnIRUQkRajQRURSRGst9BleBziNlpoLWm425Woe5WqetMrVKsfQRUTkVK31CF1ERBpRoYuIpIgWV+hmdr2ZrTOzDWb2j00sNzN7Orp8hZmNinXbBOe6O5pnhZl9YmbDGyzbYmYrzWyZmcX1d/diyHWVmR2KvvcyM3s81m0TnGtag0yrzCxkZp2iyxL5ec0xs71mtuo0y73av86Wy6v962y5vNq/zpYr6fuXmfUws7+Y2VozW21m32tincTuX865FnMj8iPUG4G+QAawHBjcaJ0bgXcAA0qAsli3TXCuS4GO0fs3fJkr+ngLkO/R53UV8Na5bJvIXI3Wvxn4c6I/r+hrjwNGAatOszzp+1eMuZK+f8WYK+n7Vyy5vNi/gK7AqOj9HKAi2f3V0o7QxwAbnHObnHMngd8BtzRa5xbgJRdRCuSaWdcYt01YLufcJ865g9GHpUD3OL33eeVK0Lbxfu0JwPw4vfcZOecWAQfOsIoX+9dZc3m0f8XyeZ2Op59XI0nZv5xzu5xzS6P3jwBrgcYXP0/o/tXSCr0Q2Nbg8XZO/UBOt04s2yYyV0OTiPwp/CUHvGdmS8xsapwyNSfXJWa23MzeMbMhzdw2kbkws7bA9cB/Nng6UZ9XLLzYv5orWftXrJK9f8XMq/3LzHoDI4GyRosSun+1tF8sauqHsRvPqzzdOrFse65ifm0zu5rI/3CXN3j6MufcTjPrDLxvZp9HjzCSkWspkWs/VJvZjcCbQFGM2yYy15duBj52zjU82krU5xULL/avmCV5/4qFF/tXcyR9/zKzbCJ/gHzfOXe48eImNonb/tXSjtC3Az0aPO4O7IxxnVi2TWQuzOwiYBZwi3Nu/5fPO+d2Rv+5F/gDkb9eJSWXc+6wc646ev9tIGhm+bFsm8hcDdxJo78OJ/DzioUX+1dMPNi/zsqj/as5krp/mVmQSJm/6px7o4lVErt/xfuLgfO5EfkbwyagD3/9YmBIo3X+hq9+qfBprNsmOFdPYANwaaPn2wE5De5/AlyfxFxd+OsJZGOArdHPztPPK7peByLjoO2S8Xk1eI/enP5LvqTvXzHmSvr+FWOupO9fseTyYv+K/nu/BPz6DOskdP+K24cbx/9INxL5dngj8H+jzz0EPNTgQ/v36PKVQPGZtk1irlnAQWBZ9FYefb5v9D/OcmC1B7m+HX3f5US+TLv0TNsmK1f08f3A7xptl+jPaz6wC6glclQ0qYXsX2fL5dX+dbZcXu1fZ8zlxf5FZBjMASsa/He6MZn7l079FxFJES1tDF1ERM6RCl1EJEWo0EVEUoQKXUQkRajQRURShApdRCRFqNBFRFLE/wdrMcXkwdnP+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(S, 'o-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "id": "6ad77f17",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.603,  0.784,  0.149],\n",
       "       [-0.578, -0.301, -0.758],\n",
       "       [-0.55 , -0.543,  0.635]])"
      ]
     },
     "execution_count": 394,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = Vt.transpose()\n",
    "V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "id": "719a8399",
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = X @ inv(X.T @ X) @ X.T @ y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "id": "c33cb0ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.176,  0.388,  0.138],\n",
       "       [-0.041,  0.286, -0.597],\n",
       "       [-0.172, -0.119, -0.488],\n",
       "       [-0.255,  0.038, -0.122],\n",
       "       [-0.357, -0.044,  0.284],\n",
       "       [-0.349, -0.349, -0.142],\n",
       "       [-0.41 ,  0.148,  0.138],\n",
       "       [-0.338, -0.033,  0.223],\n",
       "       [-0.301, -0.152, -0.247],\n",
       "       [-0.182,  0.011,  0.11 ],\n",
       "       [-0.16 ,  0.027, -0.171],\n",
       "       [-0.162,  0.488, -0.205],\n",
       "       [-0.164,  0.12 , -0.02 ],\n",
       "       [-0.232, -0.501, -0.035],\n",
       "       [-0.108,  0.194,  0.195],\n",
       "       [-0.272,  0.201,  0.134]])"
      ]
     },
     "execution_count": 402,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U3 = U[:, [0, 1, 2]]\n",
    "U3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "id": "45a614c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 403,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(yhat, U3 @ U3.T @ y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 404,
   "id": "368cd505",
   "metadata": {},
   "outputs": [],
   "source": [
    "U2 = U[:, [0, 1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 406,
   "id": "3edf12e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.176,  0.388],\n",
       "       [-0.041,  0.286],\n",
       "       [-0.172, -0.119],\n",
       "       [-0.255,  0.038],\n",
       "       [-0.357, -0.044],\n",
       "       [-0.349, -0.349],\n",
       "       [-0.41 ,  0.148],\n",
       "       [-0.338, -0.033],\n",
       "       [-0.301, -0.152],\n",
       "       [-0.182,  0.011],\n",
       "       [-0.16 ,  0.027],\n",
       "       [-0.162,  0.488],\n",
       "       [-0.164,  0.12 ],\n",
       "       [-0.232, -0.501],\n",
       "       [-0.108,  0.194],\n",
       "       [-0.272,  0.201]])"
      ]
     },
     "execution_count": 406,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 407,
   "id": "de9a4a0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat_2 = U2 @ U2.transpose() @ y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 408,
   "id": "8394e0c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABj/klEQVR4nO2dd3zV9fX/n+87cm/WTUL2AAJhhIQtMp1sxIHWhbXV2mptHbi11fqj+q3aqlVqbWtrrRsHKo7IEBAHIEM2JAEyCNl75yZ3vH9/3NyQcbPvSvJ5Ph48Qj7zZJ3P+Zz3OeclpJQoKCgoKHgvKk8boKCgoKDQNYqjVlBQUPByFEetoKCg4OUojlpBQUHBy1EctYKCgoKXozhqBQUFBS9H05ODhBDZQA1gAcxSyhmuNEpBQUFB4Sw9ctTNXCylLHWZJQoKCgoKDlFSHwoKCgpejuhJZ6IQIguoACTwipTy310dHxYWJuPj451ioIKCgsJQ4McffyyVUoY72tfT1Mc8KWW+ECIC+EoIkSal/Lb1AUKI24DbAEaMGMG+ffv6ZbSCgoLCUEIIcbqzfT1KfUgp85s/FgOfADMdHPNvKeUMKeWM8HCHDwUFBQUFhT7QraMWQvgLIQLt/wcWA0ddbZiCgoKCgo2epD4igU+EEPbj35VSbnSpVQoKCgoKLXTrqKWUmcAUN9iioKCgoOAApTxPQUFBwcvpTcOLgkK/WX8gj2c3pZNf2UBMsC8PLhnPimmxnjZLQcGrURy1gttYfyCP3318hAaTBYC8ygZ+9/ERAMVZK3TJUH/AK6kPBbfx7Kb0Fidtp8Fk4dlN6R6ySGEgYH/A51U2IDn7gF9/IM/TprkNxVEruI38yoZebVdQAOUBD4qjVnAjMcG+vdquoADKAx4UR63gRh5cMh5bOf5ZfLVqHlwy3jMGKQwIlAe84qgV3MjCpEjbWK9mYoP1PH3VpCG1KKTQex5YNK7DtqH2gFcctYLb2JddjgQWJ0UC8OHtcxUnrdAt46ID23weG+w75B7wSnmegtvYk1WORiW4cfZINh8v4lh+9ZB6fVXoG7syygCYFBuEyWJl4z0XeNgi96NE1ApuY3dWOZPigjhnZAhCwLH8Kk+bpDAA2HGqlNFh/kwZHkRRtdHT5ngExVEruIWGJguHcyuZNSoUf52GUaH+HM+v9rRZCl6OyWJlT1Y5c8eEEhmop6LeRKPZ0v2JgwzFUSu4hf05FZgsklmjhgGQFGPgmOKoFbrhcG4ldU0W5iaEEWnQA1Bc3ehhq9yP4qgV3MLurHJUAmbEhwA2R51X2UBlfZOHLVPwZnaesuWn54wOJTLI5qiHYvpDcdQKbmF3ZhnJMUEE6rUAJMcEASjpD4Uu2ZFRSlK0gRB/HyINOgCKlIhawZ2sP5DHvGe2MeqRFOY9s23Qzi4wmiwcOFPZkvYASIo2AHC8QHHUCo4xmizsP13JvDGhAEQG2iLqwiEYUSvleR5iKE2SO5xbRZPZysxWjjo8UEdEoE7JUyt0yr7sCposVuYmhAEQ7KfFR6OieAg6aiWi9hBDadDM7swyhKCNowZIjjEoJXoKnbIzoxSNSnBu8++NEIJIg07JUSu4j6E0aGZ3VjnjIwMJ9vNpsz05JoiMkjqMpqFXbqXQPTsyypgyPJgA3dkX/8hA/ZBMfSiO2kMMlUEzJouVH09XtMlP20mOMWCxStILazxgmYI3U200cSS3knkJoW22RwbplfI8Bfdx1/yEDtsG46CZI3lVNJgszBod2mFfUoxtQVHJUyu0Z3dmOVYJc5rz03YiA/VK6kPBfRh8bWkAfx81YFtcG4yDZnZnlgMd89MAw0P8CNRplDy1Qgd2ZpSi06iYPjK4zfZIg466Jgs1RpNnDPMQiqP2EFuOFxHsp2X9HfMABq0G3J6sMhLC/QkL0HXYp1IJJsQYlBI9hQ7sPFXGufHD0GnUbbZHtTS9DK30h+KoPYDZYuXr9GIuHh9BQngAfj7qQdn4YbFK9mVXOEx72EmOMZBWUIPFKjs9RmFoUVLTSHpRDXPHdPy9iQgcmt2JiqP2APtzKqmoN7FwQqQtqow2DEpHfTy/mppGs8OFRDtJ0QYaTBaySmvdaJmCN7Mr09Y2Prddfhpo1Z2oOGoFF7M1tQitWnDBONsvYlK07fXfOsiiyt1Ztj+4WaO6iqhtreTKgqKCnV0ZpQTqNUxsXmxujX0wk5L6UHA5X6UWMXt0aKu5FwZqG82cqaj3sGXOZXdWOSND/Vryio4YGxmAj1o1KN8oFPrGjlNlzBoVikbd0T356zQE6jRKRK3gWrJK68gsqWNBYkTLtsFYpma1SvZml3eZ9gDQqlWMiwoYVF+7Qt85U15PTnl9y3wPR0QMwe7EHjtqIYRaCHFACPGFKw0a7GxNLQJgwYTIlm3jIgNRq8SgiipPFNdQWW9iZhdpDztJ0bZWcikHV+pHoffYZbcc5aftRAUNvVrq3kTUq4BUVxkyVPjqeBGJUYEMH+bXsk2vVTM2ImBQ1RPb66e7i6jBlqeuqDcNydZghbbszCglLMCHcZEBnR5ja3pRctQdEELEAcuBV11rzuCmsr6JfacrWNgqmrZjX1AcLOzOKiM22LfNAwkgJTOFxesWM/mNySxet5iUzBSS7amfvMHz9Sv0HiklOzLKmJMQhhCi0+MiDHqKa4yDbvG9K3oaUb8IPARYOztACHGbEGKfEGJfSUmJM2wbdGxPL8FilSyYENFhX1KMgaLqRkprB36kIKVkT1bH/HRKZgqrd66moK4AiaSgroDVO1dzuvH7ZrFbxVEPZTJKaimpaeww36M9UQYdJoukYgipA3XrqIUQlwLFUsofuzpOSvlvKeUMKeWM8PBwpxk4mNiSWkRYgI4pccEd9g2mBcWMkjpKa5vatI0bzUae3vM0Rkvb9IbRYuSVI38nPtR/UKV+FHrPjlPd56fhbIneUEqV9SSingdcLoTIBt4D5gsh3napVYOQJrOVb9JLWJAYgUrV8bUuOdoLpakOfwAvTITVwbaPhz/o0Wkt9dOjQ8mrzeOFH19g0bpFVDU6dsSFdYUkKa3kQ54dp0qJC/FlRKhfl8fZtROH0hS9bh21lPJ3Uso4KWU8cD2wTUp5o8stG2TszS6nptHsMO0BEOSnJS7E13uiysMfwOd3Q9UZQNo+fn53j5z1DxmlhIZl89dDj7Dso2W8fux1ZkTOYJje8cJilH8UyTEGcisaqKofWsN2FGxYrJIfMsuY203aA1o3vQydiFqR4nITXx0vQqdRcd7Yzl/rvGpBcesTYGonYmBqgA2PgM4AukDQBdg++gSCLpAaaxOfZnzGN/WvIcNLOFw6jF9N+hXXjLuG6IBoUrb/gdVZn2Bs/UYhJRf5xp3tUCyo6vbVV2HwcSy/imqjmXljuv/ZhzcP+BpKqY9eOWop5XZgu0ssGcRIKdmaVsS8MWH4+XT+LU+OCeKr1CLqGs346zz8DK3Kdby9oQzWXtdm00mtlvcMAXwe4E+DSkWixcLVDTquNPvgc/wbOPUj6AJZfuwT0FpZExJMoUZNpMWCzmrlw5I9TBp3GLClfhRHPfTY2Vw/PaeLAV52fDQqwgJ8hlSJnhJRu4ETRbWcKW/gNxeO6fK45BgDUkJaYTXnjOy+/tilBERAbZGD7VGw8l1MDZV8XbSHtQXfs68uBx+hYpnvcObVBlGfW8WSMX74YIT6cqjMgcYaaKpleRMsrzvbKl8jBL+MjuSPux8mNPRXHM8ffKNeFbpnx6lSxkYEEGHofNxAayIC9UNK5FZx1G5gS0s3ouP8tB175cfxfA876sIj0FhLir8/a0KCKNSoiTJbWFVdz6wLH2Zd6V4+PPEhxfXFxAbEcu8593LlmCsJ0Ydw/weH+NqnmKt+vhDa18K+MLE5532WQCl5pV7LLyLiyAp/lQNFfsBUt32pCp6nyWxlb3Y5180Y3uNzIg06JfWh4Fy2pBYxOS6oZRGkM6KD9IT4aT1bold4BN64nJSgEFYbfDBKMwAFWg2/Dw1CHn0RiWRuzFz+MPsPnB97PmrV2eHuu7PKmBk/zHHDwoLHbQuSrXPfah9C5j/Ov8dcxBUf30BJwN85WjKTieGJrv5KFbyEAzkVGE1W5vYgP20nKkjPkSHUIKUMZXIxpbWNHDxTyYLEjt2I7RFCkBRj8JyjLjgMb1wGWj/WRES2OGk7ViR+Gj8+X/E5ryx6hYuGX9TGSedVNpBb0cCs0Z28DUy+Fi77GwQNBwSodWC1wrAEIvwiuCPxL0irltu3/Jqc6hwXfqEK3sTOjDJUAmb3YC6MnYhAPWV1jZgsnfbgDSoUR+1itqUVIyUsTOo67WEnOSaI9KIa9/8CFhyi5K3L+TwggN9Nnk9Bg+Pu0npzPfFB8Q737Wmun3akj9jC5Gvh3qOwuhLuOw5BMfDeDVCdz/mjxtOQ8yuaLGZu3XwrhXWF/fyiFAYCOzNKmRgbRJCftsfnRBr0SGlTgxkKKI7axWw5XkRMkJ6k6I5D0B2RHGOgyWwlo8T1iidGs5Gd+Tt5fvvD/CRlJfMjA/l9gGBn2WH0asdpmij/qE6vtzuzHINeQ2JUz75W/MNg5fvQVAtrVzI8APxFDDN9H6GmqYZbN99KWUNZX740hQFCfZOZAzmVzOlB/XRrooKGltKL4qhdiNFk4buTpSyYENnlkJnW2B26KzoUpZSkl6fz+tHX+fVXv+a8987j11/9mneyUwhBzb0TbubDyz7k62u/ZvXc1R2ctV6tZ9X0VZ1ef3dWOTNHDUPtoPOyUyKT4Cf/hYJDqD67gwnRgeQVDePlhS9TVF/Er7/6dacdjQoDnz1Z5Zitknm9LMk8q504NCJqZTHRhezKKKPBZOm22qM1o8MD0GtVHMuv5qrpvbtfSmYKa/avobCukCj/KFZNX8Ws6Fnsyt9l+1ewi9KGUgASghK4JuZC5h76hHOEP343fQEhI1uutXz0coAO17Nvb09xtZGs0jpWzuz5yn0L45fCoj/CV4/zm+hQfluwiMlhc3nx4he5c+ud/Hbrb/nPov/gp+26tVhh4LErowytWjAjPqRX5w217kTFUbuQLalF+Puoe/Vap1YJxkcZet1Kbp9MZx96VFBXwO+++x0S2yjIEF0Is2NmMyd6DnNi5hBVmQdvXgm+wXDzFxA8osM1l49e3qljbs/uLPv86d69wrYw924oTuXiQ69yscWPrNLzmBszl2cveJb7v7mfu7fdzcsLX0an1vXt+gpeyY6MUqaNCOmyEcwRof4+aFRiyDhqJfXhIqSUbE0t5vyx4eg06u5PaEVyjE2VvDeKJ2v2r+kwmU4iCdQG8v6l77P9uu385YK/cOXYK4mqaHbSfiFwc4pDJ91bdmeVEaDTtMyW7jVCwGVrqI+cwfPaf5F7fCcAC0Yu4Ml5T7K7cDcPbH8Ak1WZBTJYqKxv4lh+dY/me7RHpRJEBOqGTOpDcdQu4lh+NYXVxl6lPewkxxioNprJrWjo/uBmOquQqDXVkhSahEo0/6hz98FbK1o56T6kKhywO7Occ0aGOBQk7TEaHZob3qUcA9N3/BZqbF/TZQmX8disx9ieu51Hv38Ui9XiFJsVPMsPmWVISY/mezgiwjB0JLkUR+0ivjpehBAwP7H3jrplQbEXA5o6q8Zos/3MXnjrSvALhZu/hKC4XtvmiLLaRk4W13ZdltdDfIIieTrocXzMNbB2ZUtzzHWJ13HvOfeyIWsDT/7wpKKvOAjYmVGGr1btcD57T4hSHLVCf9maVsT0ESGEBvQ+p5oYZUDVS8WTC+Iu6LCtTZXGmT02J+0fZoukg5w3U2Nvti0/PbuzRpde4jdiKr8Xd0P+fvj0Tmh2yrdMvIVbJ93KRyc/4vl9zyvOeoCz41QpM0cNw0fTNzcUOYTUyBVH7QIKqho4mlftUBuxJ/j6qBkdHsDxHi4omq1mfij4gUi/SKL9oxEIov2jWT13tW0xMGc3vHWVbdCSk500wA+Z5ei1KibFBjvleskxQXxcP5Xqeb+Ho+vgu+da9t017S5WJq7kjeNv8MrhV5xyPwX3U1RtJKOkrk/5aTsRBj3VRjMNTYM/FaZUfbiAranFACzsQ37aTnKMgb3NlRTd8XnG55yuPs2LF7/IghEL2u7M+QHe/gkERsFNn4Mhps82dcaerHKmjwjpc2TUHvuC5J6Ym1g4OQO2/R+EJ8KEyxBC8MjMR6gz1fHywZc5XXWaH4t/7FEJoYL3sDPDViba1/w02FIfYHP68WH+TrHLW1EiahewJbWIkaF+jInoXPK+O5JjDORXGamo61rA02Qx8crhV0gKTWL+8Pltd57eZYukA6Pgpi9c4qSr6k2kFlb3vSzPAYnRBoSA44U1ttkgsTPg49tss0gAlVDxx7l/ZGLoRL7I+qKDWG5KZorTbFFwDTtPlRHkq2VCDzt2HTGUaqkVR+1k6pvM7MwoY0Fiz7sRHZFk11DsZkHxk1OfkFebx51T70Qc+fCsxuGzY+GNy23O+eYUMET32Zau2JtdjpR0PoipDwToNGfFbrV6uP5d8A2xLS7W2EbGalQaSo2lHc41Woys2b+m7zfvo06kQs+RUrIzo4w5o0N718XajkjD0FF6URy1k/nuZClNZmuPhzB1xllV8s7z1I2WRl45/ApTw6dyXll+W43DumKwmmDWr20RtYvYk12Oj1rF1OHBTr1umymCgZGwci00lMP7PwWT7Q+zqM6BsAGdlyp2Sz90IhV6Tk55PXmVDcwd07+3sKEkcqs4aiez5XgRgXoN58b3L8Ic5u9DdJC+y8qPD9Ntw/vvmnYXYtuTHTUOkbCjH9FlD9idWcbU4cHotb1r6umOFrHbhuYGl+gpcOUrkLsXPl8FUnZakigQvHnsTRrMPa9Dx9wIm37vWCdy6xN9/CoUHLHjlG3QVn8l1wJ1Gny1aiWiVugdFqtkW1oxF4+PQNufxo9m7B2Kjqg31fOfI/9hZtRMZkbP7FzjsLPtTqC20czR/Gqnpj3sOBxOlXQ5XPwYHH4PdrzIqumrOgyO8lH5MCpoFM/ue5ZlHy3r2mHXFML+N+G9n8JfRkOd49GurvweDkV2ZpQSEagjIbx/C4BCiCFToqc4aidy8EwlZXVNfepGdERSTBAZJbUOy4/eS3+PcmM5d06707ahs4VCJzW1OOLH0xVYrNIpjS7taVElb5/6ueABmPgT2PJHljcJVs9d3aYk8Yl5T7B+xXpeX/o6Y0LGtDjsN469QUNTHeT9CF8/Ba9cCM+Ph8/ugvyDtjnZfp1EeC78Hg41rFbJrowy5o0J69cajp1Ig35IpD6U8jwnsjW1CI1KcNE4JznqaAPWZrHbaSPOTherbarltaOvcV7seUyLmGZrCAmIhOq8thfQ+trkr1zE7swyNCrBOSN7N/msJ4QH6ogI1HV8oxACrngZyrPg41tZ/svNLL96c4fzz4k8h1cXv8qPZ77ln/ue57l9z/Hanme5paKSa2vr8Y091/a9GbcUIpJs1x0xp6NUmIu/h0ONE8U1lNU19Xr+dGdEGvQcPFPplGt5M4qjdiJbUos4N35Yr5QqusJeT3y8oK2jfiv1Laoaq7hzanM0/f0Lti6+iVfDmd22V/WgOJuDmXytU2xxxO6scibFBfV68llPSY4xOK560fraKkH+Mx9ev9RWGVJdcPZrjj0HTmyEE5s45/ROXrWa2B84jH9GxPBcqOC1mNHcMukGrhl3TdvRqfbv1dYnzorwXvCgS7+HQ42z+WlnOWpb6kNK6ZQI3VtRHLWTyCmr50RRLY8td86QI4C4EF8Mek2bBcWqxirePPYm84fPJzksGU5tsTmWiT+Bn7zaUfnbRTQ0WTicW8kvzxvtsnskxRj49mQpRpOl42KlIRpm/AK+/hPYA+CqM7Z66+bRroQnwpzfwtglTB8+i/+oNewv2s8/D/3TFmEffY1bJt7S1mFPvtb2r64UnhsHja5X2hlK7MooZWSoH3EhzpktHmnQ02i2Ut1gdlqA5I0oOWonsSXVViq2KKlvbeOOsIvdtn79f+PYG9SZ6rhj2h221/91v7S9ul/+ktucNNiUo00WySwX5KftJMcEYbFKThTVOD5g/5sONkrQB8OqQ3DHblj0BMTPA7UtJpkeOZ3/LP4Pby57k/Eh43lu33Ms+3gZrx99nXpTPSmZKSxet5jJ6+azOH4kKSfWtcwaUegfZouV3Znl/a72aI296WWwV34ojtpJbE0rYkxEACNDndvKmhwTRFphNRarpNxYztupb7Mkfgnj/GPh/RsBCde/DT7ubaH9IasclaDXyhy9IbmllryTEsXOqjGMVRAS3+W1p0VM49+L/93isJ//8XnmfzCfx75/7GynI2ZW+1pJ2f+PfnwVCnaO5FVR02h2WtoDhk53ouKonUC10cTuzPI+D2HqiqRoA0aTlcySWl478hqNlkZ+M+U38NndUHQMfvIaDHNd+qEz9mSVkRwTRKDeda+bw0P8CNBpOteP7KwaoxdVGnaH/daytzBLM2ZpbrPfqFKx5vgbPb6eQufszLDlp521kAht530MZrp11EIIvRBijxDikBDimBDij+4wbCDxTXoJZqvs1xCmzkiOtUWVP+Rk8V76e1w6+lJGp26wTZWb/xiMXej0e3ZHo9nCgZxKl5TltUalEiRFdyFLtuBx28Jia/pYpTE1YipNFsdzVQotDUr6wwnszCglMSqQsD6M/u2MCMPQUCPvSUTdCMyXUk4BpgJLhRCzXWrVAGNLahHD/H3aVGY4i4TwAHw0Kj7JfBOL1cLtw86BzX+ACZfB+fc7/X494dCZKhrNVpfmp+0kxRhILajBYnXgKCdfaxvaFDQcELaPl/2tz1UanYovmM22+muFPmM0WdiXXeHU/DSAXqsm2E876CW5unXU0oZ96Vvb/E8JL5oxW6xsTy/h4vER/Row0xlatYqEqCZONWxhxYhFDP/iAQgdAyv+6dbFw9bszrS9wro6ogZbnrrBZCG7rM7xAZOvhXuPwupK28d+lNI56nTUq3WsqqqFox/3+boKsD+ngkaz1an5aTuRgYNf6aVHOWohhFoIcRAoBr6SUu52tiHrD+Qx75ltjHokhXnPbGP9gbzuT/IC9p2uoKrBxKJ+DmHqChGyBQnclr4TLCa4/h3QBbrsft2xJ7ucxKhAgv18XH6vpO4WFJ3I8tHLWzod7awYcyXLo+fBsU/AanW5DYOVnafKUKuES8YNRAyBNvIeOWoppUVKORWIA2YKISa2P0YIcZsQYp8QYl9JSSczEzph/YE8fvfxEfIqG5BAXmUDv/v4yIBw1luOF+GjVnH+2HCXXD+nOodc87fMrPIjuuCIbTBR2FiX3KsnmCxWfjxd4Za0B8DYiEC0atHlFEFnsnz0cjZfvZlDPz/E8MDhpFekw8SroCYfcve4xYbByM6MUibFumbx2aadOMRTH62RUlYC24GlDvb9W0o5Q0o5Izy8d07r2U3pNJjazrNoMFl4dlN6r67jbqSUbEktYk5CKP461/QO/evQv/BB8Ex1KplJd0LiJS65T085kldFfZOFWaOd/wrrCB+NinGRgZ1XfrgIlVBx/fjrOVB8gNTw0aDWKemPPlJjNHEot4p5/Rxr2hmRBj0ltY2O1zEGCT2p+ggXQgQ3/98XWAikOdOI/ErH08062+4tZJTUkV1W75JqD4DMyky+yPyClVVVHG6awuchP3fJfXrDnmZ5sP6Oce0N9imC7hazXTF2Bb4aX97N/BTGLoLj68E6+PX5nM3e7HIsVun0hUQ7kQYdFqukrHbwRtU9iaijga+FEIeBvdhy1F8404iYYN9ebfcWtjZ3I853Qf00wD/2PY+v1cotIpgXDQ/YpKk8zO7MMhLC/QkPdF6JVXckRRsoq2ty++utwcfAZaMv48vML6lIXAK1RZCzy602DAZ2nCrDR6NyyfAuaN30MoQdtZTysJRympRyspRyopTS6VPUH1wyHt92sxx8tWoeXDLe2bdyKltSi0iKNhDrggdKeskRNuV9y421DYRc9y4jY2PcsqDWFRarZF92hdvSHnaSYzsZeeoGViaupMnaxEeyFjS+SvqjD+zMKOOcESFOF5ewMxS6E72iM3HFtFievmoSBr0tzxsdpOfpqyaxYlqshy3rnPK6Jn48XeGytMffv7qbQIuVn1/4J4iYQJJd8aTe5JL79YTUgmpqGs1uW0i0M6FZ7NYTD6oxIWOYFTWL9zM+wTx2MaR+BhZz9ycqAFBW20hqQbXL8tMwNOZ9eIWjBpuz/uu1UwF4+afTvdpJA3ydVoxVwkInDmGyc/S7Z9huKuWm4IkETb4eODtIvzuxW1fyQ3P9tDMVx3uCXezW3QuKdlZOWElhXSHbYxNtKjCnv/eIHQORHzJtaxpzXJSfBggL8EEloFhx1O5hTEQAABnF3j9acmtaERGBOiY2O1CncWYvfz/2GsGouHH5f1o226WpPPH6b2d3VjkjQ/2ICtJ3f7CTSYo2cKzAM1/7hXEXEu0fzbs16aD1V9IfvWBHRikBOg1T4vr5d9KFOrxGrSIsQDe0c9TuJC7EFx+1ilMl3u2oG80WvkkvYcGESFTO7EasKWL/xz9nh6+OWybfhr/O0LKrU8UTN2G1SvZml7s97WEnKcbAmfJWYrduRKPScN3469hbvJ+TYy+E1M9tjUcK3bIro4yZo4ah6Y+GaA/U4SMNeiX14S40ahXxYX5kFHfSLuwl7M4sp67J4tz8tLkJ+eHPeUlvJdQniOsn3dLhkE4VT9zAieIaKutNzHRz2sOOfeRpqoe+/p+M/Qk6tY61gf7QUA5Z3/TovIHacesM8isbyCqt63/b+NYnSPERLI6LYXL8cBbHxZDiI9qow0caBncbuVc5arANIcr08oh6S2oReq2KeWOcmHfb/Ci7iw+wT+/DrVN/g6+mYyVJUoyBk8W1GE3ur+Xd3Zxr9FREfVbs1rWOujPHGqwP5pJRl/BF2SGq9EFw9JMeXWugdtw6A/tY0/7WT6eYy1kdNowCrQYpBAVaDavDhpFiLm85JtKgo7hm8KY+vE6Ka0xEAJuPF9FktuKj8brnCFJKtqYWc96Y8P6XGx3+oEWfTwIvjRpHpF8QV4+72uHhdsWTk0W1TOpvzq+X7MkqJzbYl+HDnCOh1FvCA3WEB+pcmqO3O1Z7l6zdsQIsnxzNJSOv5pNTn/BuTDK/OvYZG0Y8QGWjoKreRFVDx38ni2qwtOvRsXfcevtiuTPYeaqUYf4+JEb1by7NmtBhGNulGI0qFWtCh7G8+fNIg57yuiYazRZ0GteUAXoSr3PUCeEBWKySnPI6xkR4bvBQZ6QV1pBX2cDdC8b070L2vFuz4vV3vnoOY+Tx0AXo1I6bSc4qnlS51VFLKdmdVeayeSY9JbmdLJmz6WyUwb3vH+Se9w8C4Dsynv9qSrjNVM0n697ha+s0APx81AT5agny1WLw1TJ8mB9pnTQoeXvHbX9ZfyCPZzelkVdpxFer4rND+f16MBWqHa8Dtd5uFxAorm70WDDhSrzSUQOcKq71Ske95bitG/HixH7mp7c+0eKkJfD3kGDiTCZWHPoC5j/j8BS74om764kzSuoorW3yWNrDTnKMge87E7t1Ap05UAncv2gcQX5achpv4P3TT/GNYRh/i86k8bIHMei1Dt/+5j2zjTwH17QPux+MtH8raTBZW95K+uqsw1BTQsd0X1SrKYf272lxjdEjjjolM4U1+9dQWFdIlH8Uq6avYvno5d2f2EO8LrcwOtym/ZdR4l0Livbc5fNfnUCrFuxslr3vE011zSvYNrb5+ZKq8+H2ymq0nekAclbxxN0Liruzmuun3dyR2J7kmCDMzakfVxAd7LjsMDbYl7sWjOXnc+J5+IKrifCL4P3I4QRmf0WYTnaaonPUcQtglZKKOsdqMgMdpw9Yyz/I8IaOvkCv1rNq+qqWz1uaXqrcn6dOyUxh9c7VZ7U26wpYvXM1KZkpTruH1zlqf52G6CC9V9VSt14UAjBZZN8XhTK/gX/OBSDF34/FcTHcExGGWkqElN3q/dkUT6rdOilsT1Y54YE64kM9+0rp6lryacODO2xrP8pAq9Jy7bhr2WmuIMtaDxlbO72eveM2NtgXgc3h371gDFUNZn75xl6PLAq7GmcPWDvzzVMc0uuYFzWLKL8okOAnBavnrm4TsXpSO3HN/jUYLW3va7QYWbN/jdPu4XWOGmwLihleVPnhlCjBWAWfr4I3LwehImX6NawOC6VAqwEhsAjBk2GhpEy7ssvLJMUYqG+ycLozxRMnYnuL2MqnB/OpNZr49GC+y+/ZFSOGuS71czy/ms3Hi5g6PKiNY3U0yuAn436CVqXlvZCwbptfVkyLZccj88l6Zjk7HpnPfYvGs+a6qRw4U8ndaw8MutGcTh2wVpLOf0v3ohZqnjj/Kb665isu8o0hyGLikqg5bQ4N9tPio1ZRVON+R11YV9ir7X3BKx11QngAGSV1bh9r2Rn9jhJObIKXZ8P+N2HuXXD7DtY05ThYyRasKe1aPCfZTYonZ98ibL/49lyjJ0vLXJX6aTJbue+DgwT5+vDazTPbOFZHedUw3zCWxC/hU39f6k5saFlr6CnLJkXz+KVJbD5exB8/P+Y1v+fO4N6FY2m/9NfXAWuF3zzNp4H+XJlwORF+tjWh+fGLKdBoOH50bZtjhRA2pZcq9zvqTrU2O9neF7zUUftT22j2mpbQPkcJ9eXw8W3w7rWgD4JfboHF/wc+fn1+Cp9VPHGto/ZWMQdXpH7+tvUkaYU1PHPVJIb590xe7IbEG6jDwqc6ASc39/qev5g3itsuGM2bu07zr28ye32+t2KySiQQ6u/T5VtJt5Rn8VrBNyBU/HLqb1s2XzTx56ikZKuD/K+nlF5WTV+FVrRdi9ALbZscen/xUkfdPPPDS9IffRrDemw9vDwTjn4EFz4Mv/4G4s5p2W1o1R7emu6ewj4aFWMjAl2+oOitYg721E+nYre95EBOBf/Yfoprzonr1YCtSeGTmBQ6kbXBwViPfNSnez+yNJHLpsTw541pg6IJpsls5e/bTjF1eDD7HlvY5VtJd5R892c+CvTnipFLiQ44W90R4h/ODHUg2+pyoN2biKe6E5fX1pHQaETVvM4UbTKzurSM5bXOS096paNuGc7kJY56xbRYfr88seXzLqOEmiJ4/0b48CYwxMJt2+Hi34PmbEnWsbJj1DbWomr37W+/kt0ZtnriKpe+MnurmIM99eOMeuqGJgv3f3CIKIOeP1yW1OvzV064gWyNih9yvrZV8vQSlUrw3DWTmTM6lAfXHeL7k6W9voY38dH+XPIqG7hn4ViE6McMnOp8Xs/ZjEWo+OX0uzrsnh85iwyNICvzqzbbPSVyW7XtCU5pNfy8qobD2WfYnJvP8urKNi3u/cUrHXV4oI5AnYZTXlT5MTrM9vB4+5ezHEcJUsLBtbYo+sRmWLgafrUVoia1OazCWMF9X99HuH84v5/9e6L9oxEIov2jO6xkd0ZSjIHS2iaXtsze48RcozNxZurn2U3pZJbW8ew1UzD0QXR1SfwShmkDeTdAByc29skGnUbNv352DqPDArj97R89Oh2xP7SOpi8c17/GqPLvn+PDAF8uibuI4YbhHfYvmHwzAFuPt81TRxn01DVZqG1077zwreYKzEKwtK6+zXbZRaltb/FKRy2EYLSXVX7Yu8wSox004VSegXeuhvW3Q3gi/GYHnHcvqNv2E1msFh769iFKG0p58aIXuW78dWy+ejOHbzrM5qs397hAvmU2tQvz1GV1TUhss377lWt0Mnax2/46tF0ZZby2I4ub5ozs88wWH7UP1ySu5FtfPWeOvNdnW4J8tbx+y7kE6jX84n97ya2o7/4kL8Np0XRtCW9lfIpRpeJXM+5xeEhU1FSSLSq2lR1us/1sLbV7o+pPA4IYbjKR1NS2Nr4I580C8kpHDbYFRW+aopdeWE1YgA9hAa26yqxW2Ptf+MdsOL0Llv0FfrEBwsY6vMbfD/6dHwp+4NHZj5IcltxnWyY0PyxcFX1V1DXx8tenuHh8OPseW9SvXKMrSIrun9htbaOZBz48RHyoHw8vS+z+hC64NvE61ELF+6X7obHvmpbRQb68/ouZNJgs3Py/vVTWD5yGGGdG01U7X2RtgJ4l0fMYHTS60+MWBidyRDRRWJHRsq2lO9GN6Y+yhjIO6DUsqjW2eQOtlz483XSN0+7jxY46gMJqo9tfYzojIvszNsjfnh1cvvMleOMySLkP4mbAb3fBrF+DyvG3dOvprbx65FWuHnc1V429ql+2BOq1jAz1c9mC4stfn6Ku0cwjyya45Pr9JTnGJnbb19TPn1KOU1DVwPPXTsHPp39TFCL8IlgYfg4f++upP/5pv641PiqQf/9sBjll9dz25o8DpiHGadF0QyXvnPiAOpWKW8+9r8tD54//CQDbDr3Wsq2l6cWNtdRbTm9BCsk5tRqapAarFORaw3jE9Cv2GRY57T5e66jtC4reMPLUeugD7qj5G+GWYloGl29+DHJ/hMtfgp+th5CRnZ6fWZXJozseZVLYJH4383dOsSk5xuCSEr0z5fW8ues0P5kex/h+Tj1zFf0Ru/06rZi1e85w2wUJnDPSObNLVk6/gxq1ipRjb/X7WnMSQnn+2insyS7nvg8OYvXyhhhnRtO1P7zM234+zI+YwbiQcV0eO3r8CkaZLGzL+65lmyfayDdkbyBSG8355hKeNV/L6MZ3OK/pb3ylvtCp6zle66hbD2fyNJYtq/EVDl5F/UJg+s+hiyiizlTHPV/fg06t468X/RUfdc/qdLsjKdrA6bJ6qo3OVRp5fnM6QsB9i7v+Q/EkE+yt5Hm9e1BV1jfx8EeHGRcZwL2LHKen+sK0yHNI1Bh4t+E0sr6i39e7bEoMjy2fwJdHCnky5bhXN8TYo+lV/Y2mG2tZe+wNatQqbpv5QPfHa3xYoItin6mCygbb99xfpyFQp3Fb5UdRXRH7i/bzk4BoBLDBOtNl6zle66hHhvqhUQmvWFDU1HTSOl1T0OV5Ukr+sOMPnK4+zbMXPOvUTiX7gmJaQd/zou05mlfF+oP53HLeKKKDPFuG1xU2sVu/Xr9R/L/PjlFe18Rfr53q1JnFQghuGPMTTvlo2bf/X0655q/OH80t80bxvx3ZvPpdllOu6Wzs0fSU4cFc1M9oun7vv3nTT8P5oZNJDu3Z+s2CkQuwCPgm7awkV4RBR7GbUh+bT29GIllalElhwATyieD4E0tdsp7jtY5aq1YxItQ7ZLlqdJ00QnQzQOl/x/7HV6e/4r5z7mNm9Eyn2pQU4/wBRc9sSCPYT8vtFyY47ZquIjkmqFc5+i+PFPDpwXzuXjCWibHOn+W9bPpvCLbCu5mfOe2ajy2fwPJJ0fzpy1Q+O+TZOSuOcFpu2mTkw0P/oVKt5tezHurxackTf0qk2czWjC9atkUa9G6r+tiYvZEJQQmMyjvMDp/ziA/zx9fHNaIFXuuoAcaEe0eJ3odBt9BEuzpbrS8seLzTc34o+IE1+9ewJH4JP0/6udNtigjUERbg47QSvW9PlPD9qVLumj+WIN/e1xS7m6QYAznlPUv9lNQ08ugnR5gcF8RvLnLNQ0iv9eWqwLFss9ZQUJrqlGuqVILnr53CzFHDeOCDQ+zK6MdoXSfjzGjauP91/qcXzA4ez5TwKT0+TwyLZ4FVx86609SbbCWN7mojz63J5XDJYZZobOscH9Sfw4Qox93GzsCrHXVCRADZZXWYLVaP2vF2/SxydGMBYfsXNBwu+xtMvtbh8QW1BTz0zUOMMoziiblP9C/a6AQhBEkxQU5ZULRaJU9vSCMuxJcbZ49wgnWup6cdilLaRtLWNVl4/popaPujht0N1029HYAP9jzvtGvqtWr+87MZjAz147a39pFW6Blx3/Y4LZq2mPjox5cp06i5bebDvT59QeS5NCLZkbMNgAiDnuIao8vz+puyNwGwtCATS9QUdlcaWspmXYF3O+rwAEwWSU655xoAGposlJaVMNKUATNugdWVcO/RTp10o6WRe7ffi8lq4sWLX8RP67oZzknRBk4W19Bk7t+DbP3BPFILqnlwyfgBozeX1MMpgh/tz2NLahEPLRnP2EjXVrHEJCziYpPgo5K9NFqcF9UF+Wl5/ZaZ+Pmoufm1vR6ft9JktvLy186JppsOvstrOjPTA0dxbvS5vT5/+oSrCbZY2Jr2IWATuTVZJOUuFmbYlL2JySGJxOYdoCB2CXB2kdsVeLmj9rzay4miGpaq9qC1NsKUld0e/9TupzhWdow/nfcn4oPiXWpbcowBk0VysrjvC4pGk4XnN59gUmwQl02OcaJ1riUiUE94oK7LiDq/soE/fnaMmfHD+MW8Ua43SghuiLmACqxsaHYcziI22Jf/3TyT2kYzV/1jB3Oe3tpBKd1dfLw/l9wKJ0TTVguf7n2BYo2GX896pE+X0MRfyIVGE9+WHcZkMbUSEHBd+iOrKovU8lSWqkMA2Od3AQCJnnTUQojhQoivhRCpQohjQgjnze7rhgQvGM6UXljDT9Tf0RQ82tbY0gXrTqzj45Mfc+ukW5k/Yr7LbetpVNkVb+7KtqltL0tEpXJ+isaV2GrJHS+mSil5aN1hLFLy3DVTULvpazt32m2MaWri3aP/c/rrd1KMgZ/PGUlhdSMFVUYkZ5XS3eWsm8xW/u6kaNp07GP+q2lksn8cc2LmdH+CI7R6FgYmUCPN7C3cS4Qbml42Zm9EIFhccAqiJrG3OhiDXkNMkGMpN2fQk4jaDNwvpZwAzAbuEEL0ftRYHzDotUQE6jwqy1Vw+gSzValopt3QZb30kZIjPLX7KebFzOOOqXe4xbZRof74+aj7vKBYWd/E37ed4sJx4czt47wLT5IcY+BUcS2N5o4dfG//cJrvT5Xy6PIJjHCjhJiInsxKiy+pxmIOlRxy+vUdqey4c06406JpKUnZ9RfytBpum/lQv641J+FSfK1Wtpz8mMjmNnJXCQhIKdmYtZFzQicSmfsjJK0gtaCaCdEGl6xF2enWUUspC6SU+5v/XwOkAm4b+pAQHsApD0bUMTm2tmDVlOs6PaasoYx7t99LhF8Ef77gz6hV7snzqlSCCc1zL/rCP7ZnUNNo5pF+zrvwFEnRjsVus0vreOrLNM4fG8YNM928OCoEl465gkCrlXePvNb98b3Ek3PCnRlNW9K/5FVVDRN8I7lg+EX9upZu/DLOazDyde53hAXYGspclfo4WXmSzKpMlqqDAbBOuIK0whqX5qehlzlqIUQ8MA3ooBclhLhNCLFPCLGvpKTESeZBQoQ/GcW1nunOkpKZ1ZvJ8J8GwY7/4M1WMw9++yCVjZW8cNELBOmcX6PbFXZpqt62GudW1PP6zmyumhbn8l8yV5HsoJbcYpU88OEhNGrBX66e7NIopzP8Jl3Hippavsr9huL6Yqdeu7N54J0pqDsTZ0bTG3c8xWmtltvOfaD/P6Nho1mAP6WWetIqjhLq7+Oy1MfGrI2ohZqF+ScgciJnVDHUN1lcWvEBvXDUQogA4CPgHillhxBOSvlvKeUMKeWM8PD+PW1bkxAeQLXRTGmt+6eJVZzYyUgKyBu5otNj1uxfw97CvTw+53EmhLp/iFFyjIHaRjNnejka86+bTwDe3SreHY7Ebl/9LpN9pyt44opkz3VXRiaxUhOBRVpZd2KdUy/tSG0IYESIn0vngjgzmrZmfs1/rKWM0YUyP35x/40TgguGX4RGSrZmbybCoHdJ6kNKycbsjcwKn0romb0taQ9wbcUH9NBRCyG02Jz0O1LKrmWXnYwn1V6M+96mQfrgM2mFw/0bszfy+rHXuX789VyecLl7jWumLwuKx/Or+eRgHr+YF0+shxVb+oMt9RPY8rWfKKrh+c0nWJIcyYqpnh3JOjzpJ5zf0MCHae9jsjhvHsuKabE8fdWkNkrpS5Ii+CGrnMc/O+qyN8+WaHpBP6NpYOt3T5Lh48Nt59yHSjin8Cxw3CXMajCyJWsDkQbXRNTHy45zpuYMS1XNTjnpCo4X1KASMM7FpZ89qfoQwH+BVCnlX11qjQM8NpzJ3MiwrM/ZaD2XscOjO+w+VXGKx3c8ztTwqTx0bs/bXp3NuMhA1CrRq1byZzamYdBr+e2FY1xomXtIjgkitaCaRrOF+z44SKBew5+unOSRlEdbw67ihuoaShvL2Xy69+K3XbFiWmwbpfR//WwGv75gNG//kMOTX6Q63Vm3RNNxQVw0vn/RtDz9A/9uyideG8TihJ4JZfSI+PNY0NBErrEUP/8Sl+SoN2ZvRKPSMD8vHSKSIHwcaQXVjArzR+/gLceZ9ORxNg/4GTBfCHGw+d8lLrWqFVEGPX4+avdH1Cc2ojPXsEU7n9DWYgFATVMN92y/Bz+NH89f9DxatedarvVaNWMjAnq8oPj9yVK+PVHCXfPHEOTn/a3i3WE0W6hvsjD+sY0czavm8qnRbcUdPEX4OOYYEoiXat5Ne9eltxJC8MiyRH4xL57XdmTx543pTnXWZ3PT4/r9APzm29Wk6Xz41fS7nbvorgvg4tBJCAk1moOU1jZicmJHs1Va2Zi9kXkRMwjK2Q1JKwBILax2yxpPT6o+vpdSCinlZCnl1OZ/X7rcsmZUKsHocH/3N70ceo8yEUp19FwAUjJTWLxuMZPfmMzFH1xMTnUOz1/0PBF+Ee61ywFJ0T2bTW1rFU8lNtiXn83pfH72QGH9gTw+2d+2fvi9Pbleo+itSr6SSbVVHC45zOQ3JrN43WJSMlNcci8hBI9fmsRPZ43gX99k8OKWk065rlOj6fyDvNKQSawmgEvGXekU+1oTNnYpUxuNFBp3IiWU1jovqj5UcojCukKWCgMgIekKaowmzpQ3eIej9gYSwgPcW0tdV4o8uZmPLXMZFx1MSmYKq3eupqCuAImk0dKIWqWmsK7QfTZ1QVKMgeKaRkq6UTz57FA+x/KreWDJuAHTKt4Vz25Kp7Fd+7w7a4q7IyUolK/8bTXcEklBXQGrd652qbN+8oqJXDsjjjVbT/Ly16f6fU1nRtO7tv8/jup0/Grqb9GqXPA2N2YhC+oaKDSdQWjLnJr+2Ji1EZ1ax8V5x226qBGJpDfrqLq64gMGiKMeEx5AXmUD9U1ukuU6sg5hNfOh6TwSowJZs38NRkvbxQmz1cya/WvcY0832BcUuxr72Wi2ObCkaANXTPEO7cP+4sma4p6w5uR7GNtJsxktxv7/3hz+wCYHZ5eFO3x2HrNKJXj6qslcOS2WZzel8+p3mX2+jVOj6eI0XqlJJVLtyxWJ1/frWp0SnsiC5oU+TeAxp407tVgtbMrexAWRM/E/vets2sNNFR8wQBx1Qossl5vSH4fWUhWcxAk5nMQoQ6eRs7dE1MnR3UtTvbXrtK3d+JKB1yreGZ3VFHe23d109vtRUFfA0dKjfbvo4Q/g87ttcnB2WbjP727jrNUqwbNXT2b5pGj+LyWVN3Zm9+lWzoym921fzX69jlsm3eq6NR0hiBu9gPFNZrSBR50mILCvaB9lxjKWigDsaQ+A4wU1BPlqW+aLuJKB4ajD3ViiV5wKBQfZH7IUlYCxkQGdKrM4U7GlPwT5aYkL8e10QbGq3sRL205x/tgwzh/rvBp3T+OopthXq3aqVl1/6Oz3QyBYmbKS67+4nvWn1mM098ChWMxQXw5fPQ6mdm8MpgbY+kSbTRq1ihevn8qipEj+32fHeHd3Tq9sd2Y0TXkWr5T/SJhKx1XJP+vftbpjzCIW1NWi9s0hu8I5gdTG7I34afw4P/cohI2DCFu/RFphNROiA91SYdQ/CWY3ER/mh0q4aYreofdApSHFOpf4UC16rZpV01fxhx1/wGQ9Ww+rV+tZNd1t86m6JamLVvJ/fHOKaqNpwLaKd4Zd7ujZTenkVzYQE+zLg0vGO10Gqa+sCpvF6ppPMLZ6g9FbrTwy7Fwag2J5v3AHf9jxB57b9SRX6odzrSaM4SYTNFaDsbrtR1M3DU1VuR02adUq/n7DNG5/60ceXX8EH42Kq8/pWpXIjj2afvKKif12RAe3/5HdvnoeSLoJvcbF0efoC5lf38g/QuBY5S5s44n6jslq4qvTX3FR9Bx8v3sDzn8AhMBqlaQX1nDducOdY3c3DAhHrdOoGTHMz/URtdUCh9+HMYvYl6chKca2SLB89HLePv42x8uOI5FE+Uexavoqlo92Yh1oP0mOCeKr1CLqGs34687+WPMqG/jfjmxWTI1t0VkcTKyYFus1jrk9yw98AuYy1oQEU6hRE2W2sKqikuWnPwFgJbBXr+M9QyBvWRp5g1PMM6tZSQDzfMJQB8WCzgD6oOaPBvjmL9BQ3vFmfsNAyg6Dw3QaNf+88Rx+9cY+Hlp3CK1acEU3zUAmixOj6ep8XineSYhfANdM/mX/rtUT9EGMi5pOlCmPM9Y9/b7cD/k/UNVYxTL8QVpb0h6ny+ttreMuVHVpzYBw1OCmyo+sb6CmgMaFT3H6SD1XTrNFH3WmOk5VnuKa8dfw2OzHXGtDH0mOMSCl7XXsnJHDWrb/dfMJkHD/AG4VH7BU5bIcyfK69tGwgNu2I/QGZuqCmKkLpKixgo9OfsS6E+u4o6GE2IBArh1/EVeOuZIQfcjZU/1CbTnpNukPAfVl8N9FsPCPED+vzd30WjX/+fkMbv7fHu774BA+ahXLJnVs4rLjzGj62DdP8r2vjlWJK10qotEaMXYRS/f/jdeDUqlpqiHQp+9VGRuzNxLoE8jcM4chdAxE2oR309y4kAgDJEcNtgXFzNI6LC6cZ8Ch90AfRJphLlLC+CjbD3j7me0YLUaWjVrmunv3E0et5KkF1Xx8IJeb5o4kLsR9oz4VmulM/DgoDmKmwrDR4B8KGh8i/SP57dTfsunqTTx74bNE+0fzwo8vsPDDhTz6/aMcKTlia2KZfC0p825l8YjhTI4fzuIRw0m54E64/CWoyoPXL4F3roWiY21u6euj5rWbz2Xq8GDuWnuALceLHJpmslh5adspJjsjmq4r5ZXcrRjQcP3U3/TvWr1hzEIW1NeDsPBt7rd9vkyjpZFtOdtYEDMPn+zvbdUezQ+u1ILqljUsdzBwHHW4P01mK3kVLiq9aqyB1M8h+SrSSm0DoBKbHfXGrI1E+kUyLWKaa+7tBKKD9IT4advkqf+8MY1AnYY7Lh74reIDkgWP20SQW9ONKLJWpWVp/FL+t/R/fHz5x1w59kq2nN7CDV/ewPUp17N652pW526kQC2QQlCgFqzO30RKcCjc9SMsXA1nfoB/zoNPbofKs4uI/joN//vFuSTHGPjtO/vZnt5xsp8zJuSlbP8Di1+byOQPL+JrPx0z/WMJ8HGPQwMgchJjrIEEmFVszt7S58vsyNtBramWZdKvTdoDbBUfo8MDXN46bmcAOWoXV34c/8y2YDNlJWmFNfhqbXnxqsYqvs//nqXxS502QMYV2MRuz3Yo7jxVyvb0Eu64eAzBfj4etm6IMvlamwhy0HB6IorcnrEhY3ls9mNsvWYrj856lEZzIx+d/KhDTX9LbbaPH5x3L9x9EObeBUc/hpfOgU2P2ipGsIlxvHnLLMZGBvDrt35kx6nSluu0jqYvHt+3jtuU7X9gddYnLQ8SgO9qsknZ/oc+Xa9PqFSURp7Hgvp6dubv6FlVjQM2Zm0kRBfCzNMHbW8/UZNa9qW5qXXcjvd6nna4fDjTobW2H8bwmaQX1jAuKhCVSrA1Zytmq9mr0x52kmOCSC+yid0+vSGN2GBfbpob72mzhjaTr7WJIXcjitwVAT4BXJ94PZ9c8Umnx7Sp2fYbBoufhLv32+73wz9gzRT49jloqiPIT8tbv5zFqDB/fvXGPp7fnM68Z7Yx9tEN5FY0MGvUsN5H06YGyP6eNRkft6lyAWhUCdZkdm67KzCOvJhL6qsxWhrYlb+r1+fXm+rZnrudRbHno8n+rk3ao9poIreiwS0diXa8xlG3nqXhaCZCiL8Pof4+romoK3Mg+zuYshIJpBXWkNg8tnBD1gaGBw4nKdQt6mP9oqHJTJPZyrjHNnAkr4oLx4W57dVMwfUIIYj2d7wI2GbB0U5QHFzxMvxmJ8SfB9uehL9Nh33/Y5hexVu/nEWgXs1L206R16qb8+0fTnc/L6WuDNK+hM2PwasL4enh8PpyCtWOHXyhmz2NdtwCzqlvxA8NW3O29vr8b/O+pcHcwFLpC9LSJu2RVtDcOu6mig/wEkfdfpZGZzMREsIDXOOo7V1dk6+lpLaR8romxkcFUtpQyp7CPSyNX+r5sZndsP5AHh/sa1tL+8mBPK8ZUKTgHFZNX4Ve3bYWWSAoN5Zz3/b7KKgt6HhSxARYuRZu2QQh8fDFPfCPWYSf2YgKweWq7/ne524ydTfwvc/dLLJ823ZeipRQkQ0H18Jnd8PfZ8Kzo+G9lbD7FRBqjLNv57ULb+/U7ijnDbLrEWER0RyXCcxpFHyT+w1ma+/GT2zM2ki4bzjTs3+0fc+ip7TsSyt0b8UHeEl5nqNZGva8W+ta5YQIfzYdc7xa3WektFV7jJwHIfGkn7TJiCVGBbI5ezNWaeWSUW6b6tpnHA8osvLspnSvrTNW6D32v4c1+9dQWFdIlH8Uv536W4rqinj1yKt8l/sdt06+lZuTb8ZH3W5tYsRsuGUjpG+ArX+ED37Oe9YIorUV6IStmStOlPKM9lX+XlMCu09Czi7I+QFqmh8AuiAYMQumXA8j5mCJmkxK7lZeOvAShXWFJPpGkVVfSGObJh/JqtHOn5bXFQa9hu+ZyiWVX7JVF8b+ov3MjJ7Zo3Nrm2r5Lvc7rk24HPWWv8GcO9vUp6cWVBPip20R0nUHXuGoezpLIyE8gPK6M5TXNTHM30kLZHk/QtlJmHc3QMtErPFRgfzru42MCR7DmBDvr5rw9gFFCs5j+ejlDputLku4jGf3PstLB17i01Of8vDMh7kg7oK2BwkBiZfAuCVwaC3DP70LNW0f8H6iiYe0H8KGD8EQZ0ubjJgNI+ZA+ARoHjS1M28nf910E+kV6SSHJvOneX9iZvRMUrb/gTWZn1CoskXSq0ZfyfKLnnTZ98MRQgiO+83klw0foRO29EdPHfXXZ76mydrEUqserOY2aQ+wVXwkRrlWdbw9XuGoo/yjKKjr+MrWflZCQitZrmH+wzoc3ycOrQWNvmUiVlphDeGBOpoo50DxAe6adpdz7uNiYoJ92+QZW29XGBrEBMTwwsUvsDNvJ0/veZo7tt7BRXEX8dDMhxge2K7VWaWGaTei+vROh9eSgLjnKAR3bJFOLUvlhR9fYFfBLmIDYvnLBX9hSfySlqqo5Rc96XbH7Ijy4GTMRYHMVRvYmrOVR2Y+0iPnuiFrAzH+MUzO2mMTtY45W5ZrsUpOFNaw0s3q9l6Ro141fRV60Xaill5oO8zSGGMv0XNW5Ye5EY5+BImX2tpzseWfEqMC2ZS9CYCl8Uudcy8X4+0DihTcx9zYuXx8+cfce8697C7czYr1K3j54Ms0mDs+yEUnTTkiaHgHJ51fm8/vv/s9131xHcfLj/PQuQ/x2YrPWDZqmVeWrkYG+bNHNYUF5YUU1RdxrOxYt+dUGivZlb+LJXEXIrK+aVPtAXC6rI4Gk+tVx9vjFd/d5bV1rC4tI9pkBikRUvJIWRnLa9sOYYoJ9kWnUTlvQfHkZmiogCkrAdvT8mRRLeMjA9mQvYHk0GRGGNz75OwrjkRPn75qkpKfHqJo1VpumXgLn6/4nAUjF/CvQ/9ixfoVbM3Z2lamqwdNOVWNVfx131+57JPL2Hx6M7+Y+Au+vOpLfpb0s455cC8i0qDjK9MkLqooRi1UPar+2JqzFbM0s8yqA6up5U3bTqq94sONC4ngJakPtj7B8upKlldXcljnw09jomi0mm2jG1vVnapVglFhTpTlOvQeBETC6IsAyC6ro9FsJXxYDcdPHueBGQ845z5uwpsHFCl4hkj/SP5ywV+4Ztw1PLX7Ke75+h7mxczjkZmPEB8Uf/bva+sTtgl8QXE2Jz35WposTaxNW8u/D/+bmqYaLk+4nDun3ek14327I9Kg55OmSfxZZWWGLpKtOVu7nXi5IXsDIw0jScz8AYJGQOz0NvtTC6pRqwRjItzYaYmXRNStRzRObmxisrGRd4ICsToY3TgmIsA5TS91ZXBiE0y6BtS255V9IbHI8gMAS+KX9P8+CgpewLlR5/LBZR/w8LkPc6jkEFd+diUv/vgi9aZ6UgL8WTw8hsmjRrB4eAyf+/vyReYXXL7+cp7b9xyTwifx4WUf8n/n/d+AcdJgc9SlBGEMm8SC+gayqrLIrOxc8aa0oZS9hXtZGncRIvNrSLq8wzTCtMJqEsJdrzreHu+IqIPimhUrbPysuoYHI8L4LiyOC9sdmhAeQMqRAowmS/++WUc/sr3aTL2hZVNaQTUqIdlXuo3pEdMH1C+lgkJ3aFVabky6kaWjlvLCjy/w36P/5YP0DzBajC2z1gvqCnj0+0eRSCYMm8D/W/T/mBMzx8OW943IZuWV4sjzmJ/+Kk8Nj2ZrzlZGB492ePxXp7/CKq0ss/g4THuALfVxzkgHzUUuxjsi6nZ5sgV19USazbwV3rELKyEiACltaYp+cWitrXe/eWwh2Co+4iKryarOHBAt4woKfSHMN4w/nfcn3lr2Fg2WhjaCGGAT4g3WBfPepe8NWCcNtNQ5ZxhmE2k2Mck/rss89cYsWzluQuZOW1li3Iw2+6vqTeRVukd1vD3e4ajbDa/RBsawssHKbmMhJ05taHNoQrg/ABnF/XDUJemQv79lEdFOelENAcOOoBZqFo1c1PfrKygMAKZGTMVitTjcV9VY5ZWVHL3BHlGnacaDzsB8q5ZjZccc9m0U1hWyv3g/y4ZfDBlbbbXTDtIe4B7V8fZ4z0+i9fCa+1O5+pp16KXknS33QelZ2fvRYQEI0c8peofWglDDxKtbNtU1mskpr6NStZdZ0bMI9Q3txxejoDAw8HY90P6g16oJ8tVSUGuB0RexMN/mRxxF1S3luGYtWJo6NLmAe1XH2+M9jrodQZGTuGzEEr7Qqyl/63KoOA3YBqDHBvv2fUHRarHN9hizAAIjWzafKKpB6HKpsRQNmNppBYX+4mh2iLfpgfaHSIOOwiojjFlIfGUeCQFxbMvZ1uG4jVkbSQpNYkTGdxAYA3HndjgmtaCGYf4+RAS6r3Xcjtc6aoCfTv8tTUKwTmuCNy+H6nygn8OZsr+D6ryOaY/CGrSGQ2hUWhaMXNBf0xUUBgTLRy9n9dzVRPtHI7BN51s9d7VX6YH2h0iDnqKaRhizEID52jD2Fe2jwljRcsyZmjMcLTvKsriL4dQWWzSt6uga3ak63h6vdtQJwQnMi5nHe6FRmOrK4M0roK6UhPAAMkvqsPZFluvQe7bBMuPbLhamFlShDTrM+THnYfBx/6uNgoKnWD56OZuv3szhmw6z+erNg8ZJg81RF1cbISgWIpJYUFGMVVrZfmZ7yzH2tMcSswYsjQ7THharJL3INuPDE3i1owa4MelGSpoq2bTgPqg8A2+uYEKIhQaThYLqXio3NNbalFySV3ToxjpQvB+hqWbZaKXaQ0FhsBBp0FFc02jTWh2zgKSc/UT7RbVJf2zI2sDU8KlEn/oaAqNh+KwO18kqrcNosnokPw09cNRCiNeEEMVCiKPuMKg9c2PmMipoFG+X7EFe9zaUprPs4J3409D7PHXq52Cq65D2kFKS3bgDNToujGtfua2goDBQiTLosVglZXW29IewNLHAMIad+TupN9WTWZnJiYoTLB0+35b2mHB5p2kP8EzFB/Qson4d8NjqmkqouHHCjRwrO8bBoFC45nX8y47wX5/nyC4o6d3FDq21DQEfMbvN5oKqOiy+hxkTMNNtkvYKCgquJ8Le9FLdaBvTqvVnfoORJmsT3+V9x8bsjQgEi02A2egw7QG2ig+NB1rH7XTrqKWU3wLlbrClUy4dfSkGHwNvHX8LEpfDla8wU5XGnH2rbBPwekJVLmR9a4um2y0GfHbiW1SaOuYPX+wC6xUUFDyFvZa6sMoIGh2MuoDp2fsI0YWwNWcrG7I2cG7UuYSf3Gab+9MuiLOTWlBDQngAOo1npO2clqMWQtwmhNgnhNhXUtLLSLcb/LR+XD3uarbmbCW/Nh8x+Rr+EbiKcTV7YN0tYDF1f5HDHwASJl/XYde2M5uQFj1XJS50qt0KCgqeJarZURfVNK9njVmAujKHBL9oNmRtILs6m/TyNFLyvm1Oezh2xGkF1R5Le4ATHbWU8t9SyhlSyhnh4eHOumwLKxNXIhCsTVsLQM7Iq3hW9UtI+wLW/8ZWH925cbZqjxFzYNioNrsaLY2crN2FxjiZKINnXmsUFBRcQ1iAD0JAUXXzm/eYhaT4+3Go8qwmZFVTNatDAkgJdSwcXFnfRH6VkUQPLSTCAKj6sBPlH8XCkQv56ORH1JvqSQgP4OX6BRgv+AMc+dAm2Ck7KdfL3w+l6Tadt3Z8n/c9ZhoYqZ/r2i9AQUHB7WjUKsICdBRVNUfUw0axJiwMk2wb2BlVKtbkbnJ4DU/NoG6Nd0zP6yE3TriRTdmb+CzjMxLCLwbg+JhfMl02wHfPgdYflj7dIQfNofdArXM4DevLzA1Isz8zIjt2Iik4B5PJRG5uLkZjL8spBzF6vZ64uDi0Wm33Byv0iyiD/mzqAyjsJDztTLvV0xUf0ANHLYRYC1wEhAkhcoH/J6X8r6sNc8SU8ClMCpvEO6nv8OJ5tqL8jOJaps9/DEz18MM/wMcfFvzh7EnmJjiyzrYI6Rvc5nr1pnq+yf0GU81UkqY6SYNRoQO5ubkEBgYSHx/vka4ub0NKSVlZGbm5uYwaNar7ExT6RaRBR17lWUcdpQuhoLGiw3GdzTdJLagm1N+H8AD3t47b6UnVx0opZbSUUiuljPOUkwabsvCNE24kuzqbnIb9+KhVNrUXIWDJUzD9Jltk/d3zZ0869RU0lHeonQbYfmY7jRYj5qrJjI/y3NNysGM0GgkNDVWcdDNCCEJDQ5U3DDcRYdBT1Ko5btX0e9G362ruar5JakENE6LdqzrengGTo7azKH4REb4RrE17h/gwv7NNL0LApS/ApGttskIf3wYvTIT3bgChgvqOFYYbsjfgpxoGjaM8Vh85VFCcdFuU74f7iDLoKa9rotFsy0svH3clq5v0RJvNCCmJNltYHbfUYeu82WLlRFGNR9MeMMBy1GBTqVg5YSVr9q9hWtilZBa1ejKq1LDin1B6Ag6/f3a7tELKPbaOo2aNuKrGKr7P+54w63wCQgPcLq2joKDgHuwCAiU1jcSF+MHhD1helM3y1oIJhf+BYZPaaLTCWR1VT834sDPgImqAq8dejU6to1b3NafL62kyW8/uVGugvrTjSaYGW6TdzLacbZitZurKJnm07EahI+sP5DHvmW2MeiSFec9sY/2BPE+bpDCAsXcntqQ/tj5hk9pqTTv/YOe4F1R8wAB11MH6YC5LuIws47dYRS055e3UXqo6+cNuJZa7IWsDsf5x5BeHkRip5Ke9hfUH8vjdx0fIq2xAAnmVDfzu4yOKs1boMy1NL/Zaagei2Z1tT/Nw67idAZf6sHPjhBtZd2Id2uA9nCo+nzERrZxtO7HcNtuBsoYydhfuZvmIn5KGUBYS3cgfPz/G8fzqTvcfyKmkyWJts63BZOGhdYdZuyfH4TlJMQb+32XJDvfZOXLkCLfffjs7duwAYP/+/TzwwANs29ZxiLzC4KJNGzl06x9ak1pQzZiIAHw0no1pB2REDbZZ1bOi5qAN2cWJ4qq2O9uJ5QK2zxc8DpxVG45U2fr6PZ1/UjhLeyfd3faekpycTEZGBhaLbUHp/vvv57nnnuvXNRUGBiF+WnzUqrO11N34h9bYKz48zYCNqAF+nnwjuwvvYFfh19xN4tkd9gWBrU/YXmeC4mw/hObtG7I2kBCUQHnFMPx96okL8XVwdQVX0F3kO++ZbeRVNnTYHhvsy/u/7rsitkqlIjk5mWPHjnHy5ElGjBjB9OnT+3w9hYGDEIIIg842QQ+69Q92KuqaKKw2erziAwa4oz4v9jx0MpKTxhSkvL1tydPkazt84+Gs2vCdU+/k6901jIsKRKVSSqW8hQeXjOd3Hx+hwXS2xddXq+bBJeP7fe3Zs2ezY8cO/vGPf7Bx48Z+X09h4BBp0J9NfUCn/qE1qc0did7wxj1gUx9gm1Wd6HcJjerTHCo51KNzWmR34peQVlhDopKf9ipWTIvl6asmERvsi8AWST991SRWTIvt97Vnz57NY489xpVXXklsbP+vpzBwaN9G3hO8YcaHnQEdUQNcHLuMg+lrefXwG/x94dRuj9+QtYGk0CR8RRSV9ccYr1R8eB0rpsU6xTG3JzExEZ1Ox8MPP+z0ayt4NxEGHd+c6OHs+mbSCqoJC9AR7gHV8fYM6IgaYEJkOKbKmXyXt42C2oIuj82pzuFY2TGWxS8jrdD2tBzvBa81Cu5hzZo1PP300/j7+3vaFAU3E2nQU9toprbR3ONzUgs9O4O6NQPeUSdEBNBUPhcJrE1f2+WxG7Nteckl8UtIb8k/eccPQsF1ZGRkkJiYSENDAzfddJOnzVHwAFHtm166wdY6XusVaQ8YBI46IlCHvzqMGO25rDuxjnpTfafHbsjawPSI6UQHRJNWUEOkQUeIv48brVXwBAkJCaSlpfHf/3psnpiCh4lobiPvqaPOKq2jyWxVImpnIYQgISIAfcNF1DTV8EXmFw6PO1lxklOVp1g6yqbTm1ZYo6Q9FBSGCJG9jKiPF3hPxQcMAkcNkBDuT0FRFBNDJ/J26ttYZcfmiA1ZG1AJFYtGLsJssXKqpFZJeygoDBEi27eRd0NqQQ1atSAh3Dumag4SRx1AUXUjV4+9gayqLHbm72yzX0rJxuyNzIyaSZhvGNllttcaxVErKAwNAnQaAnSaHkfUaYXVjIkI9HjruB3vsKKf2J96o3xnE+4bztvH326z/1jZMc7UnOGSUZcAtKr4UBy1gsJQIcKg67GjTi2oZoIX+YdB4ajHRNjKrU6XNXJ94vXsyN9BRmVGy/4NWRvQqDTMHzEfgLSCGtReMBFLQUHBfUQZ9D1KfZTXNVFU3eg1FR8wSBz1yFB/NCrBqeJarh5nm1X9Tuo7AFillY3ZGzkv5jyCdEGALaIeFeaPTqOIBSgoDBUi20lydUZagV3MVnHUTkWrVjEi1I+M4jqG6Ydx6ehL+TzjcyqNlRwoPkBxfXFLtQdAelG1kp9WUBhi2AczSSm7PK6l4sNLSvNgkDhqsOWpM0ps+ok/nfBTjBYj606uY0PWBvRqPRcPvxiA2kYzZ8obFEc9RDly5AhRUVEcPXrU06YouJkog54mi5WKelOXx6UW1BAeqCPMg6rj7RlUjjq7rA6zxcrYkLEkBCXw0oGXeD/dpp349ZmvAUhXWseHNE899RQ7d+7kqaee8rQpCm6mp7XUaYXVXpX2gEEwlMnOmIgATBZJTnk9qTXfkFOT01JPbbQYWb1zNQBVJZMApXV8qLJ2rW3MwLvvvuthSxTcjV3k1jZj2rEjNlmsnCyq5bwxYe40rVsGUURtq/zIKKljzf41mNqJVxotRtbsX0N6YTUBOg2xwYpYgILCUMIeURd3EVFnltTRZLF6XUQ9aBz16OZa6oySWgrrCh0eU1hXSFphDeMiAxSxAG/m8AfwwkRYHWz7ePiDfl/yyJEjzJs3r+Xz/fv3M3/+/H5fV2HgYB9X2lWJXlqh91V8wCBy1EG+WsIDdWQU1xLlH+XwmCj/KGXGh7dz+AP4/O5m8VFp+/j53f121opmooJOo2aYvw+FXUTUxwuq8VGrGB3uXaNwB02OGmzpj1MltaxauorVO1djtJz9gejVem5K/A2P7jMp+WlPsuERKDzS+f7cvWBpF/GYGuDTO+HHNxyfEzUJlj3T5W0VzUQFsKU/ukp9pBbUMCYiAK3au2LYHjlqIcRSYA2gBl6VUnb9V+EhxkQE8NnB/JZW8TX711BYV0iUfxSrpq/C33QusFdpHfdm2jvp7rb3AkUzUSHSoOs69VFQzfljw91oUc/o1lELIdTAy8AiIBfYK4T4TEp53NXG9ZaE8ACqjWZKa5tYPno5y0cvb7P/lW9sbeVKRO1Buol8eWFic9qjHUHD4Rcp/br17Nmzufnmm7njjjsUzcQhSmSgnmP51Q73ldU2UlzT6DUzqFvTk/h+JnBKSpkppWwC3gOucK1ZfSOh1YKiI9IKa4gy6An2U8QCvJYFj4O2XUWO1te2vZ8omokKkUF6SmsbMVs6jkL2JjHb9vTEUccCrUOc3OZtXkdCRPeOWkl7eDmTr4XL/maLoBG2j5f9zba9nyiaiQqRBh1SQmltU4d9aV4sz9eTHLWjOrYOzfJCiNuA2wBGjBjRT7P6RrRBj5+PmlPFHR21yWIlo7iWC8Z6VyG7ggMmX+sUx2wnIyOD5cuXM2/ePEUzcYgTGWirpS6sNhIVpG+z73hBNRGBOkK9qHXcTk8cdS4wvNXncUB++4OklP8G/g0wY8aMrqeeuAiVSjA63J+MkroO+7JLbYXs3jRoRcE92DUTFRTsztlRG3lqQY1Xpj2gZ6mPvcBYIcQoIYQPcD3wmWvN6jsJ4QFkOIioU+0zPiK98wehoKDgeuwit+1L9EwWK6eKB7CjllKagTuBTUAq8IGU8pirDesrCeEB5FU20NBkabM9vbAatUqQEKHkJxUUhiqh/jrUKtGh6SWjpBaTRXplxQf0sI5aSvkl8KWLbXEKrSs/JsYGtWxPL6xhtCIWoKAwpFGrBBGBHWupU71QLKA13tV+4wTGdFL5kVZYQ6KX/hAUFBTcR4QDpZe0ghpb63iYd75xDzpHPTLUD5WgzYJijdFEboUiFqCgoABRDkRujxdUMzYyAI2XtY7b8U6r+oFeq2b4ML82EfWJIvtCouKoFRSGOpEORG69ueIDBqGjho6VH2ktqi6Ko1ZQGOpEGvRUNZgwmmwFByU1jZTWepfqeHsGpaMeExFAZmkdFqutnDu9sIYAnYa4EEUsQEFhqNNekqtlBrUXB3KDasypnYRwf5rMVvIqGhgR6tfSOi6EIhag4By2bdvGa6+9RlVVFQ0NDURFRXHVVVdx1VVXedo0hW6wS3IVVTcyMtS/peLDm4sNBqmjPlv5MXyYL2kF1Vw6JcbDVil4mgcffJDIyEgeeOABAKSUhISEkJ+fj5+fX6+Onz9/PvPnz+fo0aMUFhaycOFCt34tCn3HHlHba6nTCmqINOgY5u+9w9oGZeqjtaMurDZSbTQrFR8KHD16lMmTJ7d8npWVRXh4uEMn3ZfjFQYG7bUTjxd4n+p4ewalow7x92GYvw+nimtbFhITFfmtAUNKZgqL1y1m8huTWbxuMSmZ/ZtDbefIkSNtHO/hw4fbfN7f4+3nKNqM3o1Br0GvVVFUbaTJbCWjpNbrHfWgTH0AjAkPIKOklvRCpTRvIJGSmdJGRq2groDVO1cDdBCC6A0VFRU0NTURFXVWT7Mrx9vb4+201mZUq9Xcf//9PP/88322W8H5CCGINOgprG5saR339jfuQeuoEyL82XSsiNjgaqKD9AT5aT1tkgLw5z1/Jq2880l2h0sO02RtOyvYaDHy+I7HWXdincNzEocl8vDMrsUA2kfHYHO8P/3pT/t1fGJiImPHjm35XNFmHBhENncn2hcSk7w8oh6UqQ+w5anL65rYnVWu1E8PINo76e6295T09HQSEhJaPrdarezYsYPzzz+f9evXc+utt3LFFVewefPmbo9vjUajQadrO7/Yrs24evVqnnrqqX7ZreAa7CK3qQXV+GhUjPLS1nE7gzeibl5QLKgycvlUpeLDW+gu8l28bjEFdQUdtkf7R/O/pf/r831HjhzJP//5T4xGI3q9nqeeeooLLriAsLAwVqxYwYoVK6ioqOCBBx5g8eLFXR6/fv16UlJSKC4u5o477mDx4sVt7qVoM3o/kYE6vqo2klpQwzgvbh23493W9YPM0rOdiR/sPcP6A3ketEahp6yavgq9uq3yhl6tZ9X0Vf267uLFi7n44otJTExk/PjxnDx5kn/9619tjvm///s/7rjjjm6PX7FiBf/5z394/fXXef/99zvcS9Fm9H6igvQYTVYO5FQwYQAUGgzKiHr9gTye25Te8nlFvYnffXwEgBXTlAjHm7EvGK7Zv4bCukKi/KNYNX1VvxYS7Tz//PMOF/aklDzyyCMsW7asTT65s+PttHbsrVG0Gb2fiOYSvbomi9dXfMAgddTPbkqnwdRWZbjBZOHZTemKox4ALB+93CmOuae89NJLbNmyhaqqKk6dOsXtt9/e5fGdOXZFm3HgEBl4dl1hIMjzDUpHnV/Z0KvtCkObu+++m7vvvrvHx3fm2BVtxoHDkbyqlv/f/8EhHl6a6NVB3KB01DHBvuQ5cMoxwcpQJoX+01vHruBdrD+Qx3Obz6ZGC6qMXp8aHZSLiQ8uGY+vtq3klq9WzYNLxnvIIgUFBW/h2U3pGDtJjXorgzKitj8Vn92UTn5lAzHBvjy4ZLzXPi0VFBTcx0BMjQ5KRw02Z604ZgUFhfYMxNTooEx9KCgoKHTGQEyNDtqIWsG7kFIqwg2tkFJ62oQhy0BMjSqOWsHl6PV6ysrKCA0NVZw1NiddVlaGXq/v/mAFlzDQUqOKo1ZwOXFxceTm5lJSUuJpU7wGvV5PXFycp81QGCAojlrB5Wi1WkaNGuVpMxQUBizKYqKCgoKCl6M4agUFBQUvR3HUCgoKCl6OcEWZkBCiBDjdx9PDgFInmuNsvN0+UGx0Bt5uH3i/jd5uH3iXjSOllOGOdrjEUfcHIcQ+KeUMT9vRGd5uHyg2OgNvtw+830Zvtw8Gho2gpD4UFBQUvB7FUSsoKCh4Od7oqP/taQO6wdvtA8VGZ+Dt9oH32+jt9sHAsNH7ctQKCgoKCm3xxohaQUFBQaEVXuOohRBLhRDpQohTQohHPG1Pe4QQw4UQXwshUoUQx4QQqzxtkyOEEGohxAEhxBeetsURQohgIcQ6IURa8/dyjqdtao8Q4t7mn/FRIcRaIYRHpycJIV4TQhQLIY622jZMCPGVEOJk88cQL7Tx2eaf82EhxCdCiGAPmujQxlb7HhBCSCFEmCds6w6vcNRCCDXwMrAMSAJWCiGSPGtVB8zA/VLKCcBs4A4vtBFgFZDqaSO6YA2wUUqZCEzBy2wVQsQCdwMzpJQTATVwvWet4nVgabttjwBbpZRjga3Nn3uS1+lo41fARCnlZOAE8Dt3G9WO1+loI0KI4cAiIMfdBvUUr3DUwEzglJQyU0rZBLwHXOFhm9ogpSyQUu5v/n8NNgfjVXMShRBxwHLgVU/b4gghhAG4APgvgJSySUpZ6VGjHKMBfIUQGsAPyPekMVLKb4HydpuvAN5o/v8bwAp32tQeRzZKKTdLKc3Nn/4AeHRcYCffR4AXgIcAr12w8xZHHQucafV5Ll7mBFsjhIgHpgG7PWxKe17E9gtn7eY4TzEaKAH+15yeeVUI4e9po1ojpcwDnsMWXRUAVVLKzZ61yiGRUsoCsAURQISH7emOW4ANnjaiPUKIy4E8KeUhT9vSFd7iqB1Nk/fKp5sQIgD4CLhHSlntaXvsCCEuBYqllD962pYu0ADTgX9KKacBdXj+lb0NzbneK4BRQAzgL4S40bNWDWyEEI9iSx2+42lbWiOE8AMeBR73tC3d4S2OOhcY3urzODz8uukIIYQWm5N+R0r5saftacc84HIhRDa21NF8IcTbnjWpA7lArpTS/iayDpvj9iYWAllSyhIppQn4GJjrYZscUSSEiAZo/ljsYXscIoS4CbgU+Kn0vlrgBGwP5EPNfzdxwH4hRJRHrXKAtzjqvcBYIcQoIYQPtsWbzzxsUxuETUPqv0CqlPKvnranPVLK30kp46SU8di+f9uklF4VCUopC4EzQgi7iugC4LgHTXJEDjBbCOHX/DNfgJcteDbzGXBT8/9vAj71oC0OEUIsBR4GLpdS1nvanvZIKY9IKSOklPHNfze5wPTm31OvwiscdfOCw53AJmx/FB9IKY951qoOzAN+hi1SPdj87xJPGzUAuQt4RwhxGJgKPOVZc9rSHO2vA/YDR7D9jXi0e00IsRbYBYwXQuQKIX4JPAMsEkKcxFax8IwX2vh3IBD4qvnv5V9eaOOAQOlMVFBQUPByvCKiVlBQUFDoHMVRKygoKHg5iqNWUFBQ8HIUR62goKDg5SiOWkFBQcHLURy1goKCgpejOGoFBQUFL0dx1AoKCgpezv8HoZnDTsU2/wcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(y, 'o-', label='$y$')\n",
    "ax.plot(yhat, 'o-', label='$\\hat y$')\n",
    "ax.plot(yhat_2, 'o-', label='$U_2 U_2^\\\\top y$')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 410,
   "id": "78f372f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.91923137428117"
      ]
     },
     "execution_count": 410,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(y - yhat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 411,
   "id": "79b2a7fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.007402518166287"
      ]
     },
     "execution_count": 411,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(y - yhat_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "610e6214",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}