lab9/lab9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metody Obliczeniowe w Nauce i Technice\n",
    "## Laboratorium 9 - Faktoryzacja QR\n",
    "### Albert Gierlach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Użyta wersja pythona:\n",
    "\n",
    "* Python 3.8.3 64bit\n",
    "\n",
    "Użyte biblioteki:\n",
    "\n",
    "* numpy 1.18.2\n",
    "* matplotlib 3.2.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false; \n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false; \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "import matplotlib as mtplt\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "np.set_printoptions(precision=3)\n",
    "mtplt.rcParams['figure.figsize'] = [18, 9] # plots size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Faktoryzacja QR metodą Grama-Schmidta\n",
    "Na podstawie podanego pseudokodu zaimplementowałem algorytm ortogonalizacji Grama-Schmidta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_qr(A_orig):\n",
    "    A = A_orig.copy()\n",
    "    n = A.shape[0]\n",
    "    Q = np.zeros(shape=A.shape, dtype=np.float64)\n",
    "    R = np.zeros(shape=A.shape, dtype=np.float64)\n",
    "\n",
    "    def norm_col(col):\n",
    "        return col / np.linalg.norm(col)\n",
    "\n",
    "    def back_sum(U, A, n):\n",
    "        return sum((U[:, i].dot(A[:, n])) * U[:, i] for i in range(n))\n",
    "\n",
    "    for k in range(0, n):\n",
    "        Q[:, k] = norm_col(A[:, k] - back_sum(Q, A, k))\n",
    "\n",
    "    for i in range(n):\n",
    "        for j in range(i, n):\n",
    "            R[i][j] = np.dot(Q[:, i], A[:, j])\n",
    "\n",
    "    return Q, R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Później zaimplementowałem funkcję sprawdzającą zaimplementowany przeze mnie z algorytmem bibliotecznym."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|  n  | my_qr == linalg |\n",
      "| 291 |      equal      |\n",
      "| 101 |      equal      |\n",
      "| 129 |      equal      |\n",
      "| 108 |      equal      |\n",
      "| 93  |      equal      |\n"
     ]
    }
   ],
   "source": [
    "def get_random_square_matrix(n):\n",
    "    return np.random.rand(n, n)\n",
    "\n",
    "\n",
    "def compare_matrices(A, B, eps=1e-6):\n",
    "    return np.allclose(np.abs(A), np.abs(B), eps)\n",
    "\n",
    "\n",
    "format_str = \"|{:^5}|{:^17}|\"\n",
    "print(format_str.format(\"n\", \"my_qr == linalg\"))\n",
    "for i in range(5):\n",
    "    n = random.randint(50, 300)\n",
    "    A = get_random_square_matrix(n)\n",
    "    qq, rr = my_qr(A)\n",
    "    q, r = np.linalg.qr(A)\n",
    "\n",
    "    equal = compare_matrices(qq, q) and compare_matrices(rr, r)\n",
    "    print(format_str.format(n, \"equal\" if equal else \"not equal\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Następnie wygenerowałem 80 macierzy 8x8 z różnym wskaźnikiem uwarunkowania. Do tego celu użyłem faktoryzacji QR do znalezienia dwóch macierzy ortogonalnych oraz rozkładu SVD w celu znalezienia macierzy diagonalnej z malejącymi wartościami osobliwymi. Uwarunkowanie macierzy S jest realizowane poprzez zapewnienie, że pierwsza wartość na diagonali jest wielokrotnością ostatniej wartości na diagonali."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrices_num = 80\n",
    "n = 8\n",
    "\n",
    "matrices = []\n",
    "conds = []\n",
    "for i in range(matrices_num):\n",
    "    cond = random.randint(500, 2000)\n",
    "    m1, m2, A = tuple(get_random_square_matrix(n) for _ in range(3))\n",
    "\n",
    "    q1, _ = np.linalg.qr(m1)\n",
    "    q2, _ = np.linalg.qr(m2)\n",
    "    _, S, _ = np.linalg.svd(A)\n",
    "    S[0] = cond * S[-1]\n",
    "\n",
    "    matrices.append(q1 @ np.diag(S) @ q2)\n",
    "    conds.append(cond)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "W ostatnim kroku rozłożyłem wygenerowane macierze za pomocą zaimplementowanego algorytmu QR. Następnie przemnożyłem macierz \\\\(Q^T\\\\) przez \\\\(Q\\\\) i wynik odjąłem od macierzy identycznościowej. Z tak powstałej macierzy wyliczyłem normę, która oznacza błąd jaki powstał przy różnych wartościach wskaźnika uwarunkowania macierzy. Na wykresie przedstawiłem zależność współczynnika uwarunkowania od otrzymanego błędu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABBYAAAImCAYAAADniAk9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzde5hlZ10n+u+PTrBaugVMIBNJIOgRnDQokRYRmLYbGGUYBI6oRBAConm6R9HOIwqoA4ozwOg59MhxJsccZKABCRcFGUQBgaYa5RpMlCR4o3NpQIEg2g1dgYT3/LF2JUVTXVW9unbt2+fzPPXsvda+rLfqV7tqre9633dVay0AAAAAfdxh1A0AAAAAJpdgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAMIaq6teq6tWjbsfxquq6qnrk4P4vV9XLRt2mcVBVr6iq/zK4/++q6m9WeO49q+poVW3auBYCwPAIFgCYeYOD5WODg73Fr98ZdbvGXWvtha21nxp1O8ZNa+1ga+2+i8tLw5jB4ze01ra01m4dTQsBYH2dNuoGAMCY+KHW2p+t9qSqOq21dstx6zadzEHiyT4fAGCc6bEAACuoqqdV1Z9X1b6q+nySXxt0e7+0qt5WVV9Msquq/m1VHaiqL1TV1VX12CXv8XXPX2Y7966q91bVkap6Z5Izj3v8wVX1F4P3v6qqdi557EBVvaiqPlRV/1JVf1RV33wSr/2Nwfd4pKreUVVnLnn8KVV1fVXdVFW/clybbhuuUVXnVVWrqouq6oaq+tzS51fV5qp6ZVX9c1VdW1W/VFWHV/i5P6SqPjz4fj5cVQ9Za5uXea/HVdWVVfWvVfUPVfWowfpvqaq3VNXnq+rvq+qnj/veXl9V+wfbuLqqti95/IKq+ujgsdclmVvy2M7F762qXpXknkn+96AnzC8t+VmddqrtAIBxIFgAgNV9b5JPJLl7kv86WPekwf2tST6Y5H8necfgOc9M8pqquu+S91j6/Pcts43fT3JFukDhN5JctPhAVd0jyR8n+S9JvjnJs5L8QVXdbcnrn5rkJ5N8S5Jbkrz0JF77pCRPH7T9joPnpKrOT3JpkqcM3veMJOes/KPKw5LcN8kjkjyvqv7tYP3zk5yX5FuT/PskP3GiNxiEIn88+B7OSPKSJH9cVWes1uZl3utBSfYn+cUkd0myI8l1g4dfm+Tw4Hv7kSQvrKpHLHn5Y5NcPnjdW5L8zuA975jkzUlele5n+oYkT1hu+621pyS5IV2PmC2ttd9c5mm92gEA42LsgoWqenlVfaaqPrZO7/engzM0bz1u/c8Ozgq0lc5yADAz3jz4f7H49dNLHvtUa+3/aa3d0lo7Nlj3R621P2+tfTXJA5JsSfLi1tqXW2vvTvLWJD++5D1ue35rbWHphqvqnkm+J8l/bq3d3FqbTxdULPqJJG9rrb1t8Pp3JvlIkkcvec6rWmsfa619Mcl/TvJj1U0OuJbX/q/W2t8OvrfXD76fpDvIfWtrbb61dvPgfb+6ys/x11trx1prVyW5Ksl3Ddb/WJIXttb+ubV2OIPg4wT+Y5K/a629avAzf22Sjyf5oTW0+XjPSPLy1to7B9//J1trH6+qc9OFIM9urS201q5M8rJ0Icqi9w1+bremCxEWv5cHJzk9yX9vrX2ltfbGJB9e5eeyrFNsBwCMhbELFpK8Ismj1vH9fitf+8950Z8neWSS69dxWwBMrse31u6y5Ov/W/LYjcs8f+m6b0ly4yBkWHR9knus8h5LX//Pg1Bg6esX3SvJjy4NPtIdjJ59gve/Pt2B75lrfO0/Lrn/pXQhyW3f1+IDg/bdtML3seb3yuo/j+P/Px//8zzRdo53bpJ/OME2Pt9aO3IS25gbDF/4liSfbK21417bx6m0AwDGwtgFC4OzNJ9fuq6qvm3Q8+CKqjpYVd9xEu/3riRHlln/l6216065wQDMgrbKuk8lObeqlv5fvWeST67yHos+neSuVXWn416/6MZ0PRKWBh93aq29eMlzzj3utV9J8rk1vnaldt32vlX1jemGJvTx6XztMIpzT/TEdD/Pex237vif51rdmOTbTrCNb66qrT228ekk96iqOu61J7JS7U+lHQAwFsYuWDiBy5I8s7X2wHRjKP/niNsDAEt9MMkXk/xSVZ1e3eSIP5RuXPyqWmvXpxue8OtVdceqeli+ttv/q5P8UFX9YFVtqqq5wQSBSw/Uf6Kqzh8c/L8gyRsHXefX8toTeWOSx1TVwwbzCrwg/fcdXp/kuVV118G8Dz+7wnPfluQ+VfWkqjqtqp6Y5Px0w0tO1u8leXpVPaKq7lBV96iq72it3ZjkL5K8aPAz+c50wyZes4b3fH+6eSx+btC+H07yoBWe/0/p5pb4OqfYDgAYC2MfLFTVliQPSfKGqroyye9m0H2zqn64qj62zNfbR9lmACbS4qz9i19vWusLW2tfTjfB3n9I10vgfyZ5amvt4yex/SelmyTy8+kmOty/5P1vTPK4JL+c5LPpzsL/Yr72//ir0g0n/Md0Vyj4uZN47Ym+r6uT/Ey6iSU/neSf000y2McLBq89lOTP0oUWN59guzcleUySX0g39OKXkjymtfa5k91oa+1D6SZ53JfkX5K8N7f3hvjxdBNKfirJm5I8fzAHxWrv+eUkP5zkael+Jk9M8ocrvORFSX51MBRluUkme7UDAMZFfe3wwPFQVeelmyzqflX1TUn+prV29sqvWvH9diZ5VmvtMcs8dl2S7X12VgBgHFTVgSSvbq29bNRtWauq2pPkwtba94+6LQDAqRn7HguttX9NcqiqfjRJqmM2ZACYIFV1dlU9dDAc4b7peiOsuVcIADC+xi5YqKrXphu7eN+qOlxVz0jy5CTPqKqrklydrkvnWt/vYLrrSz9i8H4/OFj/c1V1ON1EUn9VVRNzlgcAJtAd0w1nPJLk3Un+KOZMAoCpMJZDIQAAAIDJMHY9FgAAAIDJIVgAAAAAejtt1A1Y6swzz2znnXfeqJvBGnzxi1/Mne50p1E3gw2g1rNDrWeDOs8OtZ4daj0b1Hl2jGutr7jiis+11u623GNjFSycd955+chHPjLqZrAGBw4cyM6dO0fdDDaAWs8OtZ4N6jw71Hp2qPVsUOfZMa61rqrrT/SYoRAAAABAb4IFAAAAoDfBAgAAANDbWM2xAAAAAOPmK1/5Sg4fPpyFhYWhb+vOd75zrr322qFv50Tm5uZyzjnn5PTTT1/zawQLAAAAsILDhw9n69atOe+881JVQ93WkSNHsnXr1qFu40Raa7npppty+PDh3Pve917z6wyFAAAAgBUsLCzkjDPOGHqoMGpVlTPOOOOke2YIFgAAAGAV0x4qLOrzfRoKAQAAAGNu06ZNuf/973/b8oUXXpjnPOc5I2zR7YYeLFTVdUmOJLk1yS2tte3D3iYAAABMk82bN+fKK69c8Tm33nprNm3adMLltb7uZG3UUIhdrbUHCBUAAABg/Zx33nl5wQtekIc97GF5wxve8HXLr33ta3P/+98/97vf/fLsZz/7ttdt2bIlz3ve8/K93/u9ef/7339KbTDHAgAAAKynhYVk9+5k27Zkz55u+RQdO3YsD3jAA277et3rXnfbY3Nzc3nf+96XCy+88GuWd+zYkWc/+9l597vfnSuvvDIf/vCH8+Y3vzlJ8sUvfjH3u9/98sEPfjAPe9jDTqltGzHHQkvyjqpqSX63tXbZBmwTAAAARmPv3mT//uTYseTQoW7dpZee0luuNBTiiU984rLLH/7wh7Nz587c7W53S5I8+clPzvz8fB7/+Mdn06ZNecITnnBKbVq0EcHCQ1trn6qquyd5Z1V9vLU2v/hgVV2c5OIkOeuss3LgwIENaBKn6ujRo2o1I9R6dqj1bFDn2aHWs0OtZ4M6j9ad73znHDlyZM3P/8YDB7Lp2LFu4dix3Pqe9+RLa3z9rbfeesJtLbe+tZbW2m2PLV3+0pe+lK985Su3PbawsJAvf/nLOXLkSObm5vKlL31p2e0sLCyc1O/b0IOF1tqnBrefqao3JXlQkvklj1+W5LIk2b59e9u5c+ewm8Q6OHDgQNRqNqj17FDr2aDOs0OtZ4dazwZ1Hq1rr702W7duXfsLdu1Kbrih67GweXM27dq15tcfOXLkhM9dbn1VZcuWLbc9tnR5586dec5znpObb745d73rXfOmN70pz3zmM2977om2Mzc3lwsuuGBN7U2GHCxU1Z2S3KG1dmRw/weSvGCY2wQAAICR2revu52fT3bsuH35FCzOsbDoUY96VF784hev+Jqzzz47L3rRi7Jr16601vLoRz86j3vc4065Lccbdo+Fs5K8qaoWt/X7rbU/HfI2AQAAYHTm5k55ToXj3Xrrrcuuv+6661ZcftKTnpQnPelJX/e6o0ePrlfThhsstNY+keS7hrkNAAAAYHRcbhIAAADoTbAAAAAwKxYWkt27k23bkj17umU4RRtxuUkAAADGwd69yf793dUKDh3q1q3zXADTqrWWwfyBU621dtKv0WMBAABgVhw82IUKSXc7Pz/a9kyIubm53HTTTb0OuidJay033XRT5ubmTup1eiwAAADMih07up4Kx44lmzd3y6zqnHPOyeHDh/PZz3526NtaWFg46QP79TQ3N5dzzjnnpF4jWAAAAJgV+/Z1t/PzXaiwuMyKTj/99Nz73vfekG0dOHAgF1xwwYZsa70IFgAAAGbF3Jw5FVh35lgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDehh4sVNWmqvrLqnrrsLcFAAAAbKyN6LHw80mu3YDtAAAAABtsqMFCVZ2T5D8medkwtwMAAACMRrXWhvfmVW9M8qIkW5M8q7X2mGWec3GSi5PkrLPOeuDll18+tPawfo4ePZotW7aMuhlsALWeHWo9G9R5dqj17FDr2aDOs2Nca71r164rWmvbl3vstGFttKoek+QzrbUrqmrniZ7XWrssyWVJsn379rZz5wmfyhg5cOBA1Go2qPXsUOvZoM6zQ61nh1rPBnWeHZNY62EOhXhoksdW1XVJLk/y8Kp69RC3BwAAAGywoQULrbXnttbOaa2dl+TCJO9urf3EsLYHAAAAbLyNuCoEAAAAMKWGNsfCUq21A0kObMS2AAAAgI2jxwIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAASBYWkt27k23bkj17umVYg9NG3QAAAADGwN69yf79ybFjyaFD3bpLLx1tm5gIeiwAAACQHDzYhQpJdzs/P9r2MDEECwAAACQ7diSbN3f3N2/ulmENDIUAAAAg2bevu52f70KFxWVYhWABAACAZG7OnAr0YigEAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvQwsWqmquqj5UVVdV1dVV9evD2hYAAAAwGqcN8b1vTvLw1trRqjo9yfuq6k9aax8Y4jYBAACADTS0YKG11pIcHSyePvhqw9oeAAAAsPGqO/4f0ptXbUpyRZL/I8n/aK09e5nnXJzk4iQ566yzHnj55ZcPrT2sn6NHj2bLli2jbgYbQK1nh1rPBnWeHWo9O9R6Nqjz7BjXWu/ateuK1tr25R4barBw20aq7pLkTUme2Vr72Imet3379vaRj3xk6O3h1B04cCA7d+4cdTPYAGo9O9R6Nqjz7FDr2aHWs0GdZ8e41rqqThgsbMhVIVprX0hyIMmjNmJ7AAAAwMYY5lUh7jboqZCq2pzkkUk+PqztAQAAABtvmFeFODvJKwfzLNwhyetba28d4vYAAACADTbMq0L8VZILhvX+AAAAwOhtyBwLAAAAwHQSLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6O22lB6tqxxrf57rW2g3r0B4AAABggqwYLCR5+hrf501JBAsAAAAwY1YMFlpraw0WAAAAgBlkjgUAAACgt9XmWDiUpK30lMHj/7219tL1bBgAAAAw/lYbCnHvjWoIAAAAMHkMhQAAAAB6EywAAAAAvQkWAAAAgN5WDBaq6oUb1RAAAABg8qzWY+FRG9IKAAAAYCKteFWIJJuq6q7pLiv5dVprn1//JgEAAACTYrVg4TuSXJHlg4WW5FvXvUUAAADAxFgtWLimtXbBhrQEAAAAmDiuCgEAAAD0tlqw8Nsb0goAAABgIq0YLLTWXnGix6rq4nVvDQAAADBRTmUoxLJXigAAAABmR+9gobX2u+vZEAAAAGDyrBosVNUFVfWaqvro4Ouyqvr2wWOrXVUCAAAAmGIrBgtV9YQkb0jyriRPS/L0JO9P8oaq+r4kbx92AwEAAIDxtVqPg+cneWRr7bol666qqvck+XiSlwyrYQAAAMD4W20oxGnHhQpJksG661trvzyMRgEAAACTYbVg4StVdc/jV1bVvZLcPJwmAQAAAJNiLUMh/qyqXpjkiiQtyfckeU6SZw+5bQAAAMCYWzFYaK29uaoOJfmFJM9MUkk+luTHWmtXbUD7AAAAgDG26uUiBwHCUzegLQAAAMCEWW2OBQAAAIATEiwAAAAAvQkWAAAAgN5WnWMhSarqgenmWfi2JP+U5NWttfcMs2EAAADA+Fs1WKiqlyQ5I8n/neRvk9w9ya9W1dlJ3tJaOzrcJgIAAADjasVgoaqekuQOrbWLquoXkjx+8NA/JflPSeaq6l9aa38w5HYCAAAAY2i1ORaekeTXBvfvmq7nwoeSbEny7iRvSHLxsBoHAAAAjLfVhkLctbX2hcH9R7bWHjy4/6dV9YHW2vOq6swhtg8AAAAYY6v1WPhMVZ07uP+Jqvq5qrpPVT0zyXVV9c1J/nW4TQQAAADG1WrBwv+b5IWD+88Y3O5NUkl+MslvJHnZcJoGAAAAjLsVg4XBpIyHq+o1Se7RWntpa+0/JXlrklcnuaW19poNaCcAAAAwhla93GRr7blV9f1Jfr2q7pHkq+muCnFpa+2dw24gAAAAML5WDRaSpLX23iTvHXJbAAAAgAmz2hwLAAAAACckWAAAAAB6EywAAAAAva0aLFTVBVX1mqr66ODrsqr69sFja5qjAQAAAJhOKwYLVfWEJG9I8q4kT0vy9CTvT/KGqvq+JG8fdgMBAACA8bVaj4PnJ3lka+26Jeuuqqr3JPl4kpcMq2EAAADA+FttKMRpx4UKSZLBuutba788jEYBAAAAk2G1YOErVXXP41dW1b2S3DycJgEAAACTYi1DIf6sql6Y5IokLcn3JHlOkmcPuW0AAADAmFsxWGitvbmqDiX5hSTPTFJJPpbkx1prV21A+wAAAIAxturlIgcBwlM3oC0AAADAhFltjgUAAACAExIsAAAAAL0JFgAAAIDeTjpYqKqPDqMhAAAAwOTp02Oh1r0VAAAAwETqEyz88bq3AgAAAJhIJx0stNZ+dRgNAQAAACaPyRsBAACA3gQLAONmYSHZvTvZti3Zs6dbBgCAMXXaqBsAwHH27k3270+OHUsOHerWXXrpaNsEAAAnsGKwUFV/naQt91CS1lr7zqG0CmCWHTzYhQpJdzs/P9r2AADAClbrsfCYDWkFALfbsaPrqXDsWLJ5c7cMAABjasVgobV2/UY1BICBffu62/n5LlRYXAYAgDFkjgWAcTM3Z04FAAAmhqtCAAAAAL0JFgAAAIDeVrsqxPOz/FUhjnegtWbacgAAAJgxq82xcN0a3+cLp9gOAAAAYAKtdlWIV25UQwAAAIDJs9pQiHuu8X2+0Fr713VoDwAAADBBVhsK8cp0cyzUCs9pSV6RZP86tQkAAACYEKsNhdjV942r6tx0YcO/SfLVJJe11n677/sBAAAA42e1Hgun4pYkv9Ba+2hVbU1yRVW9s7V2zRC3CQAAAGygOwzrjVtrn26tfXRw/0iSa5PcY1jbAwAAADbe0IKFparqvCQXJPngRmwPAAAA2BjVWhvuBqq2JHlvkv/aWvvDZR6/OMnFSXLWWWc98PLLLx9qe1gfR48ezZYtW0bdDDaAWs8OtZ4N6jw71Hp2qPVsUOfZMa613rVr1xWtte3LPTbUYKGqTk/y1iRvb629ZLXnb9++vX3kIx8ZWntYPwcOHMjOnTtH3Qw2gFrPDrWeDeo8O9R6dqj1bFDn2TGuta6qEwYLQxsKUVWV5PeSXLuWUAEAAACYPMOcY+GhSZ6S5OFVdeXg69FD3B4AAACwwYZ2ucnW2vuS1LDeHwAAYEMsLCR79yYHDyY7diT79iVzc6NuFYyNoQULAAAAU2Hv3mT//uTYseTQoW7dpZeOtk30IyQaCsECAADASg4e7EKFpLudnx9te+hPSDQUw5xjAQAAYPLt2JFs3tzd37y5W2YyCYmGQo8FAACAlezb193Oz9/efZ7JtGNH11Ph2DEh0ToSLAAATBpjhGFjzc3pLj8thERDIVgAAJg0xggD9CMkGgpzLAAATBpjhAEYI4IFAIBJYyI5AMaIoRAAAJPGGGEAxohgAQBg0hgjDMAYMRQCAGBYFhaS3buTbduSPXu6ZQCYMnosAAAMi6s3ADAD9FgAABgWV28AYAYIFgAAhsXVGwCYAYZCAAAMi6s3ADADBAsAAMPi6g0AzABDIQAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAAAAoDfBAgAAANCbYAEAAADoTbAAAAAA9CZYAAAAAHoTLAAAAAC9CRYAAJh8CwvJ7t3Jtm3Jnj3dMgAb4rRRNwAAAE7Z3r3J/v3JsWPJoUPduksvHW2bAGaEHgsAAEy+gwe7UCHpbufnR9segBkiWAAAYPLt2JFs3tzd37y5WwZgQxgKAQDA5Nu3r7udn+9ChcVlAIZOsAAAwOSbmzOnAsCIGAoBAAAA9CZYAAAAAHoTLAAAwKlYWEh27062bUv27OmWAWaIORYAABh/CwvJ3r3dZSUXJ2ecmxt1qzp79yb793eXuTx0qFtnvgdghggWAAAYf+N88H7wYNeupLudnx9tewA2mKEQAACMv3E+eN+xI9m8ubu/eXO3DDBD9FgAAGD87djR9VQ4dmz8Dt737etu5+dvH6YBMEMECwAAjL9xPnifmxufYRkAIyBYAABg/Dl4Bxhb5lgAAAAAehMsAAAAAL0JFgAAAIDeBAsAAABAb4IFAGB1CwvJ7t3Jtm3Jnj3dMgBAXBUCAFiLvXuT/fuTY8eSQ4e6dWboBwCixwIAsBYHD3ahQtLdzs+Ptj0AwNgQLAAAq9uxI9m8ubu/eXO3DAAQQyEAgLXYt6+7nZ/vQoXFZQBg7RYWuuGFBw/e/v90bm7UrTplggUAYHVzc+ZUAIBTNaVzFhkKAQAAABthSucsEiwAAADARpjSOYsMhQAAAICNMKVzFgkWAAAAYCNM6ZxFhkIAAAAAvQkWAAAAgN4ECwAAAEBvggUAAIBJsbCQ7N6dbNuW7NnTLcOImbwRAABgUuzdm+zfnxw7lhw61K2bwskAmSxD7bFQVS+vqs9U1ceGuR0AAICZcPBgFyok3e38/GjbAxn+UIhXJHnUkLcBAAAwG3bsSDZv7u5v3twtw4gNdShEa22+qs4b5jYAAABmxr593e38fBcqLC7DCFVrbbgb6IKFt7bW7neCx+V6iYcAABSqSURBVC9OcnGSnHXWWQ+8/PLLh9oe1sfRo0ezZcuWUTeDDaDWs0OtZ4M6z46xr3VryY03JkeOJFu3Jueem1SNulUTaexrzbpQ59kxrrXetWvXFa217cs9NvLJG1trlyW5LEm2b9/edu7cOdoGsSYHDhyIWs0GtV7BwkI3gdLBg7efMZibG3WrelPr2aDOs2Psa7179+0T0G3enFx0kQnoehr7WrMu1Hl2TGKtRx4sAEwsszID9GcCOoCpMezJGwGml51imF6uEz98JqADmBrDvtzka5O8P8l9q+pwVT1jmNsD2FB2inHwOb0WeyRdc03yylcml1wy6hZNn337uuEP55/f3ZqADjaO/1+ss2FfFeLHh/n+ACNlVmYMh5leeiQN39ycz8ukmrI5hmaS/1+sM3MsAPRlpxgHn9Nrx45uZ3txYkE9kuB2Dkonn/9frDNzLABAX4bDTC/d9FlP09bt3EHp5PP/i3WmxwIA9GU4zPTSI4n1NG1n+PXomXz+f7HOBAsA0JeDT2Atpu0Mv4PSyef/F+tMsAAAAMM0bWf4HZQCxxEsAADAMDnDD0w5wQIAAAyTM/zAlHNVCAAAAE7etF3xhN70WAAAAODkTdsVT+hNjwUAAABO3rRd8YTeBAsAAACcvB07uiudJNNxxRN6MxQCAACAk+eKJwwIFgAAADh5rnjCgKEQAAAAQG+CBQAAgD5cbhGSGAoBAADQj8stQhI9FgAAAPpxuUVIIlgAAADox+UWIYmhEAAAAP243CIkESwAAAD043KLkMRQCAAAAOAUCBYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAE5sYSHZvTvZti3Zs6dbBljitFE3AABO2cJCsndvcvBgsmNHsm9fMjc36lYBTIe9e5P9+5Njx5JDh7p1l1462jYBY0WPBQAm3+JO7zXXJK98ZXLJJaNuEXCqnCUfHwcPdqFC0t3Oz4+2PdPE7zlTQo8FACafnV6YPs6Sj48dO7oaHDuWbN7cLbM+/J4zJQQLAEw+O70wfQSG42Pfvu52fv724WasD7/nTAnBAgCTz04vTB+B4fiYm3MWfVj8njMlBAuMnknXgFNlpxemj8CQWeD3nCkhWGD0jC0DAI4nMGQW+D1nSrgqBKNnbBkAAMDEEiwwejt2dGPKEmPLAAAAJoyhEIyesWUAAAATS7DA6BlbBgAAMLEMhQAAAAB6EywAAAAAvQkWAAAAYCULC8nu3cm2bcmePd0ytzHHAgAAAKxk795k//7k2LHk0KFunXnibqPHAgAAAKzk4MEuVEi62/n50bZnzAgWgPGm2xmszGcEAIZvx45k8+bu/ubN3TK3MRQCGG+6ncHKfEYAYPj27etu5+e7UGFxmSSCBWDc6XYGK/MZAYDhm5sT3K/AUAhgeNaji7ZuZ+NBd/vx5TMCAIyYHgvA2iwsdF2uDx68vfvX3NzKr1mPLtq6nY0H3e1PXp/PTB8+I5xIa10gOOzfQQBmnmABWJs+B5br0UVbt7PxoLv9yduoMMZnhBO58UaBIAAbwlAIYG36HFjqoj091PLkCWMYtSNH/A4Ck8fwy4kkWADWps+B5b59yUUXJeef393qoj0+Tvaf9qzW8lR2boQxjNrWrdPxO+ggA2bLYo+/a65JXvnK5JJLRt0i1sBQCGBt+ozj1kV7fC3XTf+JTzzx82e1lqcynMHcB4zaued2QeCk/w6a4wVmix5/E0mwAKzNrB5YTqvl/mmvFCzMqlPZufGZYdSqpuN30EEGfW3UJLqsrx07uhDx2LHJ7m01YwyFgGmk2yir0U1/bdbz5+RzCf34e0Vfy3Wp97d4/M3q8MsJp8cCTKPluo3u29cvtXe5sum0XDf9D3xgtG0aR+s5nEF3buhn6efwIQ9JbrmlOyj0P4nVLNfbxd/i8afH30QSLMA0Ws9/pC5XNp38016b9fw56c4N/Sz9HO7e7X/SuBun4QfLdamfn/e3GIbAUAiYRst1G+17UONyZbA+dOeGUyegG3/jNKP/cl3qx+FvseEYTCE9FmAaLdd9+5JL+k2Es3i5MhPowKmZxKtEjNOZR0hM6jYJxin8Wa7X2Tj8LTYcgykkWIBptJ7/SKflcmUwapM4/MTOL+NmHA4KWdm4hz/j8Ld4nMIXWCeCBZgVff+RTsvlymCjTcPZfju/jJtxOChkZcKf1Y17+AI9CBYAYBim4Wy/nV/gZAl/Vid8YQoJFjbCNJy1Am7nM81aTMPZfju/AOtP+MIUEixshGk4awXczmeatZiGs/12fgGANRAsbIRpOGsF3M5nmrVwtp9EDycAZoJgYSNMw1krVmbHcbb4TLMWzvaT6OEEwEwQLGwEZ62mnx3H2eIzDayVHk5MIydUgOMIFjbC3Fz3B3fv3m6H4pJL/AGeNnYcZ4sz0cBa6eHENHJCBTiOYGGj+AM83ew4ArAcPZzGizPt68MJFeA4goWN4g/wdLPjCMBy9HAaL070rA8nVIDjCBY2ij/A082OI8DkcNZ6djnRsz6cUAGOI1jYKP4Arx87hACcCmetZ5cTPevDCRXgOIKFYXIAPBx2CGFy+bvIOHDWenY50QMwFIKFYVrrAbAd7ZNjh7Dj94ZJNAvBoM/m+HPWenY50w4wFIKFYVrrAbAA4uTYIezMwgEa02cWgkGfzfHnrDUArCvBwqlY7UB/rQfA6x1ATDs7hJ1ZOEBj+sxCMOizOf6ctQaAdSVYOBWrHeiv9QC4bwBx4ECye/fs9WCwQ9iZhQM0ps8sBINr+WzqgQYATBHBwqlY7azUWg+A+wYQVXowzLLF35MDB7rfhQMHkj17puMAxUHX9JqFYHAtf9P1QAMApohg4VSs1xnjvgHEgQO6286yxd+b3btvP0C5/vrusUk/QHHQxSRby990wyUAgClyh1E3YKLt25dcdFFy/vnd7bC79C7urF59dXe7c2cXaCS6wk+zhYUuPNi2reuRsLDwtY9P4wHKNH5PsNSOHf5+M51W+58FwFTSY+FUjLpL7yyMVR5XG9FVf3Ebr3tdcuRIcuuty5+9n8a5Fqbxe4Kl/P1mWulxBjCTBAuTbNTBxsmYtjHzG7HjtHQbi5Y7ez+NByjT+D3BUpP09xtOhh5nwzFJ+1GT1FZg3Qw1WKiqRyX57SSbkrystfbiYW6PMTZtZzA2Ysdp6TYWLXf2fhoPUKbte7KTNd3UF26nx9lwTNJ+1CS1FVg3QwsWqmpTkv+R5N8nOZzkw1X1ltbaNcPaJhtkYSG54YZu/ORad6Kn7QzGRuw4Ld3Gpk3J1q3JhRc6ez+J7GRNN/WF2+lxNhyTtB81SW0F1s0weyw8KMnft9Y+kSRVdXmSxyURLEy6vXuT+943ueaate9ET9sZjI3YcVpuG86CTiY7WdNNfeF209bjbFxM0n7UJLUVWDfVWhvOG1f9SJJHtdZ+arD8lCTf21r72eOed3GSi5PkrLPOeuDll18+lPawjq6+OkfPPDNbDh/ulufmut4LK2ktufHGbhLCrVuTc89NqobfVk7Z0aNHs2XLllE3Y7LdcENy003JV7+a3OEOyRlnJPe856hb9XXUuqcJqe8idZ4daj1FVtmPGqta2+cbmrGqM0M1rrXetWvXFa217cs9NsweC8v9Bfm6FKO1dlmSy5Jk+/btbefOnUNsEuvida/LgfvcJzuf9awuib7oouRnfmbUrWJIDhw4EJ/LU7SwkFxyydj3PlHrniakvovUeXao9exQ69mgzrNjEms9zGDhcJJzlyyfk+RTQ9weG2XfvuT1r0/OP9/4SVgLXYOnm/oCADNumMHCh5N8e1XdO8knk1yY5ElD3B4bZW6u6+Z79dWjbgkAAAAjNrRgobV2S1X9bJK3p7vc5Mtba45EAQAAYIoMs8dCWmtvS/K2YW4DAAAAGJ07jLoBAAAAwOQSLAAAAAC9CRYAAACA3gQLAAAAQG+CBQAAAKA3wQIAAADQm2ABAAAA6E2wAAAAAPQmWAAAAAB6EywAAAAAvQkWAAAAgN4ECwAAAEBvggUAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0Vq21UbfhNlX12STXj7odrMmZST436kawIdR6dqj1bFDn2aHWs0OtZ4M6z45xrfW9Wmt3W+6BsQoWmBxV9ZHW2vZRt4PhU+vZodazQZ1nh1rPDrWeDeo8Oyax1oZCAAAAAL0JFgAAAIDeBAv0ddmoG8CGUevZodazQZ1nh1rPDrWeDeo8Oyau1uZYAAAAAHrTYwEAAADoTbDAsqrqLlX1xqr6eFVdW1XfV1XfXFXvrKq/G9zedcnzn1tVf19Vf1NVPzjKtnNyquqSqrq6qj5WVa+tqjm1ng5V9fKq+kxVfWzJupOubVU9sKr+evDYS6uqNvp7YWUnqPVvDf6G/1VVvamq7rLkMbWeQMvVecljz6qqVlVnLlmnzhPqRLWuqmcO6nl1Vf3mkvVqPaFO8Pf7AVX1gaq6sqo+UlUPWvKYWk+gqjq3qt4zOK66uqp+frB+evbLWmu+fH3dV5JXJvmpwf07JrlLkt9M8pzBuuck+W+D++cnuSrJNyS5d5J/SLJp1N+DrzXV+R5JDiXZPFh+fZKnqfV0fCXZkeS7k3xsybqTrm2SDyX5viSV5E+S/IdRf2++1lTrH0hy2uD+f1Pryf9ars6D9ecmeXuS65Ocqc6T/3WCz/SuJH+W5BsGy3dX68n/OkGt37FYqySPTnJArSf7K8nZSb57cH9rkr8d1HNq9sv0WODrVNU3pfsj93tJ0lr7cmvtC0kely5wyOD28YP7j0tyeWvt5tbaoSR/n+RBYVKclmRzVZ2W5BuTfCpqPRVaa/NJPn/c6pOqbVWdneSbWmvvb91/s/1LXsOYWK7WrbV3tNZuGSx+IMk5g/tqPaFO8JlOkn1JfinJ0omz1HmCnaDWe5K8uLV28+A5nxmsV+sJdoJatyTfNLh/53T7ZolaT6zW2qdbax8d3D+S5Np0J/imZr9MsMByvjXJZ5P8r6r6y6p6WVXdKclZrbVPJ92HI8ndB8+/R5Ibl7z+8GAdY6619skk/1eSG5J8Osm/tNbeEbWeZidb23sM7h+/nsnyk+nOaiRqPVWq6rFJPtlau+q4h9R5+twnyb+rqg9W1Xur6nsG69V6+uxN8ltVdWO6/bTnDtar9RSoqvOSXJDkg5mi/TLBAss5LV2XrEtbaxck+WK6rjknsty4HpcbmQCDcVyPS9fF6luS3KmqfmKllyyzTq2nw4lqq+YTrqp+JcktSV6zuGqZp6n1BKqqb0zyK0met9zDy6xT58l2WpK7Jnlwkl9M8vrB2Gq1nj57klzSWjs3ySUZ9CKOWk+8qtqS5A+S7G2t/etKT11m3VjXWrDAcg4nOdxa++Bg+Y3pgoZ/GnS/yeD2M0uef+6S15+T27tsMd4emeRQa+2zrbWvJPnDJA+JWk+zk63t4dzehX7peiZAVV2U5DFJnjzoMpmo9TT5tnTB8FVVdV26mn20qv5N1HkaHU7yh63zoSRfTXJm1HoaXZRunyxJ3pDbh52q9QSrqtPThQqvaa0t1ndq9ssEC3yd1to/Jrmxqu47WPWIJNckeUu6P3QZ3P7R4P5bklxYVd9QVfdO8u3pJhVh/N2Q5MFV9Y2Dsx6PSDfmS62n10nVdtAt70hVPXjwO/LUJa9hjFXVo5I8O8ljW2tfWvKQWk+J1tpft9bu3lo7r7V2Xrodzu8e/B9X5+nz5iQPT5Kquk+6ybU/F7WeRp9K8v2D+w9P8neD+2o9oQZ1+b0k17bWXrLkoanZLztt1A1gbD0zyWuq6o5JPpHk6emCqNdX1TPSHZD+aJK01q6uqtenCx9uSfIzrbVbR9NsTkZr7YNV9cYkH01Xu79MclmSLVHriVdVr02yM8mZVXU4yfOTvDgnX9s9SV6RZHO6cfp/EsbKCWr93HSzSb9zcCWqD7TWdqv15Fquzq2131vuueo82U7wmX55kpdXd1nCLye5aNATSa0n2Alq/dNJfnswsfZCkosTn+sJ99AkT0ny11V15WDdL2eK9svq9p6RAAAAACfHUAgAAACgN8ECAAAA0JtgAQAAAOhNsAAAAAD0JlgAAAAAehMsAAAAAL0JFgCAoaqqV1TVjyxZfmNVfeuS5QuqqlXVDy5Zd8eqmh9cxx0AGGOCBQBgw1TVtiSbWmufWLL6x5O8b3CbJGmtfTnJu5I8cWNbCACcLMECAPA1quqpVfVXVXVVVb2qqu5VVe8arHtXVd1z8LxXVNVLq+ovquoTi70SqvM7VXVNVf1xkrsvefsnJ/mjJduqJD+S5GlJfqCq5pY8982D5wMAY0ywAADcZtCj4FeSPLy19l1Jfj7J7yTZ31r7ziSvSfLSJS85O8nDkjwmyYsH6/7PJPdNcv8kP53kIUue/9AkVxy3fKi19g9JDiR59JLHPpbke9blGwMAhkawAAAs9fAkb2ytfS5JWmufT/J9SX5/8Pir0gUJi97cWvtqa+2aJGcN1u1I8trW2q2ttU8lefeS55+d5LNLln88yeWD+5fna4dD3Jrky1W1dV2+MwBgKEyIBAAsVUnaKs9Z+vjNx712uecsdSzJXJJU1aYkT0jy2Kr6lcHrz6iqra21I4Pnf0OShTW2HQAYAT0WAICl3pXkx6rqjCT5/9u5Y9SqoiCOw7+DINhqJUG0cBc2BlyAZYxYuAF7KytbW1eRWGkrsUgfJYi4BBdgYXEs3kVe9YRLfFh8XzUwA/ee9s8wY4yb1Xl1tPSftjm0uMun6miMcW2Mcbs63Op9re4v9aPqYs55Z855b855tzqpHi/fvlX9mHP+uoJ3AQD/iGABAPhjznlZva7OxhgX1ZvqRfV8jPG5etbm7sIu76rv1ZfqbXW21XtfPVzqJ8vstpPqeKkPqw+rHgIA7M2Y82/bjgAAV2OMcaP6WD1Ybijsmj2tXs45v+3l5wCAVWwsAAB7M+f8Wb2qDnbNjTGutzkMKVQAgP+cjQUAAABgNRsLAAAAwGqCBQAAAGA1wQIAAACwmmABAAAAWE2wAAAAAKz2G/g7HxczLJbOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1296x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "norms = []\n",
    "for matrix in matrices:\n",
    "    Q, _ = my_qr(matrix)\n",
    "    norms.append(np.linalg.norm(np.identity(n) - (Q.T @ Q)))\n",
    "\n",
    "plt.scatter(conds, norms, s=13, c=\"red\")\n",
    "plt.xlabel(\"cond(A)\")\n",
    "plt.ylabel(\"|| I - Q @ Q.T ||\")\n",
    "plt.legend([\"Error\"])\n",
    "plt.title(\"Error depending on condition\")\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wykres przedstawia tendencję wzrostową. Im większy współczynnik uwarunkowania tym większy błąd obliczeń, czyli macierz Q traci swoją ortogonalność. Wynika to z tego, iż algorytm ortogonalizacji Grama-Schmidta staje się niestabilny numerycznie, czyli mały błąd obliczeń powoduje coraz to większe zaburzenia w dalszych krokach obliczeń."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Rozwiązywanie układów równań metodą QR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "W tym zadaniu mamy układ nadokreślony (mamy 3 niewiadome i 11 równań), więc musimy znaleźć taki zbiór parametrów, który minimalizuje błąd naszego rozwiązania. Rozwiązując to zadanie metodą QR minimalizujemy błąd średniokwadratowy i otrzymujemy dokładniejsze rozwiązanie.\n",
    "\n",
    "Mając równanie postaci \\\\(Ax = y\\\\), możemy je rozwiązać w następujący sposób:\n",
    "* rozkładamy macierz A metodą QR i otrzymujemy: \\\\(QRx = y\\\\)\n",
    "* wiedząc, że \\\\(Q\\\\) jest ortogonalne mnożymy układ lewostronnie przez \\\\(Q^T\\\\) i otrzymujemy: \\\\(Q^TQRx = Q^Ty\\\\)\n",
    "* ostatecznie otrzymujemy: \\\\(Rx = Q^Ty\\\\)\n",
    "\n",
    "\\\\(A\\\\) jest macierzą obliczonych kolejnych potęg x w danym punkcie \\\n",
    "\\\\(y\\\\) to wartości funkcji w poszczególnych punktach \\\n",
    "\\\\(x\\\\) to wektor współczynników funkcji \\\\( f(x) = \\alpha_0 + \\alpha_1 x + \\alpha_2 x^2 \\\\)\n",
    "\n",
    "Poniżej zaimplementowałem obliczanie aproksymacji funkcji."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[13.972  0.182 -0.434]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABBwAAAImCAYAAADwqpn0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdebyV4/rH8c/dnAYdlSPSYKpooiJHqRBROjIUIhlCOegYM0fmKfohJJJTmeWQg04J4VASiqS0JUMndTSnXT2/P+4tuzS393728Hm/Xvu12ms9az3XXrFb67vu67pDkiRIkiRJkiTlpGJpFyBJkiRJkgofAwdJkiRJkpTjDBwkSZIkSVKOM3CQJEmSJEk5zsBBkiRJkiTlOAMHSZIkSZKU4wwcJEkSACGEh0MI16Vdx4aEEIaEEG7exO03hxB+DiH8lMd15dvnTJKktIUkSdKuQZKkIiGEMA5oBOySJMmvKZeTp0IICbB3kiQztvH+Q4A5SZJcu4HbdgemAzWTJPnvdhW66Rq6A+ckSdIit84hSVJh4goHSZLyQAihFtASSICOuXieErn12PlYTWB+boYNkiRp6xk4SJKUN7oB/wGGAGdkvyGrXeDhEMLoEMLiEMLbIYSa2W5PQggXhRC+yWobuCuEUCzrtu4hhPdCCP1DCAuAviGEHUMIQ0MI80II34YQrg0hFAsh7BRCmBNCODbrvuVDCDNCCN2y1XFz1p9bZx17RQjhvyGEH0MIx4UQjgkhTA8hLAghXJ2txgNDCB+EEH7JOvaBEEKprNveyTrs0xDCkhBCl6zrO4QQJmfd5/0QQsNsj7d/CGFS1vPxDFBmQ09qCOEIYDSwa9ZjD/mt9vWOy8g6lhBC3xDCs1nP0eIQwtQQQtNsx+4eQngx6/mbn/Wz1AMeBg7OOs8v6z9nWd/3yHpOF4QQ/hlC2HW9v8fzQwhfhxD+F0J4MIQQNvpfjCRJBZyBgyRJeaMbMCzr66gQwp/Xu70r0A+oAkzOOi67TkBT4ADgr8BZ2W47CPgG2Bm4Bfg/YEdgD6BV1rnPTJJkQdb9BoUQdgb6A5OTJBm6kZp3Ib7R3w24HhgEnAY0Ia7WuD6EsEfWsauBv2fVfzBwONALIEmSQ7OOaZQkSfkkSZ4JIRwAPA6cB1QGHgH+GUIonRVUjASeAnYCngNO2FCBSZL8Gzga+CHrsbtv5GdZX0fgaaAS8E/gAYAQQnHgVeBboFbWz/50kiRfAucDH2Sdp9L6DxhCOAy4DegMVMt6jKfXO6wD0IzYWtMZOGoL65UkqcAxcJAkKZeFEFoQl/0/myTJx8BM4NT1DhuVJMk7WbMdriF+kr57ttvvSJJkQZIks4H7gFOy3fZDkiT/lyTJKmAl0AW4KkmSxUmSZAD3AKcDJEnyJvEN/BigPfEN/8ZkArckSZJJfONcBbg/63GnAlOBhlmP+3GSJP9JkmRV1jkfIYYdG9MDeCRJkg+TJFmdJMmTwK9A86yvksB9SZJkJknyPDBhE4+1LcYnSfJakiSricFGo6zrDwR2BS5PkmRpkiQrkiQZv4WP2RV4PEmSSVl/j1cR/x5rZTvm9iRJfsn6e3wLaJwTP4wkSfmRgYMkSbnvDODNJEl+zvp+OOu1VQDf/faHJEmWAAuIb3z/cDvxk/ON3VYFKJV1TPbjd8v2/aNAfeCJJEnmb6Lu+VlvyAGWZ13OzXb7cqA8QAhhnxDCqyGEn0IIi4Bbs2rZmJrApVntFL9ktSjsnvVz7Qp8n6w72frbDT3Idsi+m8UyoEzW/IvdgW+zwputtSvZ6sz6e5zPus/9+uctvw3nkSSpQDBwkCQpF4UQyhKXzrfKejP+E7H1oFEIoVG2Q3fPdp/yxFaCHzZ0O1BjvduyvzH/mbgyoeZ6x3+f9djFiasPhgI9Qwh7beOPtr6BwDTiThQVgauBTc0n+I64eqJStq8dkiQZAfwI7LbefIMaW1HLUmCH377J+pmrbuF9vwNqbGT45ua29vqBbM97CKEcsV3k+y08tyRJhYqBgyRJues44nyDfYnL5xsD9YB3ibMVfnNMCKFF1vyCfsCHSZJkX7lweQjhT1ltFhcDz2zoZFkrEp4FbgkhVMgaPnkJ8I+sQ34b9HgWcDcwNOsN+faqACwCloQQ6gI917t9LnGmxG8GAeeHEA4KUbkQQvsQQgXgA2AVcFEIoUQI4Xhiq8OWmk5csdA+hFASuBYovYX3/YgYeNyeVVOZEMIh2X6G6r8Nw9yA4cCZIYTGIYTSxFUeH2a1mEiSVOQYOEiSlLvOILYuzE6S5KffvohDCrtm+yR9OHADsZWiCXEeQHYvAx8TB0qOAgZv4pwXEj/l/wYYn/XYj4cQmhDDh25ZwcQdxE/t+2z/j8llxLkUi4lhwvqBSF/gyaz2ic5JkkwkznF4APgfMAPoDpAkyUrg+Kzv/0ecSfHilhaSJMlC4sDKx4irC5YCczZ5p9/vuxo4FtgLmJ11vy5ZN48lzq34KYTw8wbuOwa4DniBGFrsCZy8pXVLklTYhHXbIyVJUl4LIQwB5iRJcu1Gbk+IrQoz8rQwSZKk7eAKB0mSJEmSlOMMHCRJkiRJUo6zpUKSJEmSJOU4VzhIkiRJkqQcZ+AgSZIkSZJyXInNH5K+KlWqJLVq1Uq7DEmSJEmSlM3HH3/8c5IkVTd0W4EIHGrVqsXEiRPTLkOSJEmSJGUTQvh2Y7fZUiFJkiRJknKcgYMkSZIkScpxBg6SJEmSJCnHFYgZDpIkSZKkLZeZmcmcOXNYsWJF2qWokChTpgzVq1enZMmSW3wfAwdJkiRJKmTmzJlDhQoVqFWrFiGEtMtRAZckCfPnz2fOnDnUrl17i+9nS4UkSZIkFTIrVqygcuXKhg3KESEEKleuvNUrZgwcJEmSJKkQMmxQTtqW/54MHCRJkiRJueKll14ihMC0adNSq+GHH37gxBNPzJHHGjlyJF988cVW3698+fIbvH7AgAHUq1ePrl27bm9pa2VkZDB8+PC130+cOJGLLrooxx5/axg4SJIkSZJyxYgRI2jRogVPP/10jj3mqlWrtur4XXfdleeffz5Hzr2tgcPGPPTQQ7z22msMGzYsxx5z/cChadOmDBgwIMcef2sYOEiSJEmSctySJUt47733GDx48DqBw7hx4zj00EPp1KkT++67L+effz5r1qwB4kqASy+9lAMOOIDDDz+cefPmAdC6dWuuvvpqWrVqxf3338+YMWPYf//9adCgAWeddRa//vorEyZMoGHDhqxYsYKlS5ey3377MWXKFDIyMqhfvz4AQ4YM4bjjjuPYY4+ldu3aPPDAA9x7773sv//+NG/enAULFgAwaNAgmjVrRqNGjTjhhBNYtmwZ77//Pv/85z+5/PLLady4MTNnzmTmzJm0a9eOJk2a0LJly7UrOWbNmsXBBx9Ms2bNuO666zb4/Jx//vl88803dOzYkf79+9O3b1/uvvvutbfXr1+fjIwMMjIyqFevHj169GC//fbjyCOPZPny5QDMmDGDI444gkaNGnHAAQcwc+ZM+vTpw7vvvkvjxo3p378/48aNo0OHDgAsWLCA4447joYNG9K8eXM+++wzAPr27ctZZ51F69at2WOPPXIsoHCXCkmSJEkqxHr3hsmTc/YxGzeG++7b9DEjR46kXbt27LPPPuy0005MmjSJAw44AICPPvqIL774gpo1a9KuXTtefPFFTjzxRJYuXcoBBxzAPffcw0033cSNN97IAw88AMAvv/zC22+/zYoVK9h7770ZM2YM++yzD926dWPgwIH07t2bjh07cu2117J8+XJOO+20tW/as5syZQqffPIJK1asYK+99uKOO+7gk08+4e9//ztDhw6ld+/eHH/88fTo0QOAa6+9lsGDB3PhhRfSsWNHOnTosLZF4/DDD+fhhx9m77335sMPP6RXr16MHTuWiy++mJ49e9KtWzcefPDBDT4/Dz/8MK+//jpvvfUWVapUoW/fvht9Lr/++mtGjBjBoEGD6Ny5My+88AKnnXYaXbt2pU+fPnTq1IkVK1awZs0abr/9du6++25effVVIAY8v7nhhhvYf//9GTlyJGPHjqVbt25MzvqPY9q0abz11lssXryYOnXq0LNnz63aAnNDXOEgSZIkScpxI0aM4OSTTwbg5JNPZsSIEWtvO/DAA9ljjz0oXrw4p5xyCuPHjwegWLFidOnSBYDTTjtt7fXA2uu/+uorateuzT777APAGWecwTvvvAPA9ddfz+jRo5k4cSJXXHHFButq06YNFSpUoGrVquy4444ce+yxADRo0GBtODFlyhRatmxJgwYNGDZsGFOnTv3D4yxZsoT333+fk046icaNG3Peeefx448/AvDee+9xyimnAHD66advw7O3rtq1a9O4cWMAmjRpQkZGBosXL+b777+nU6dOAJQpU4Yddthhk48zfvz4tfUcdthhzJ8/n4ULFwLQvn17SpcuTZUqVdh5552ZO3fudtftCgdJkiRJKsQ2txIhN8yfP5+xY8cyZcoUQgisXr2aEAJ33nkn8McdDza2A0L268uVKwdAkiQbPe+CBQtYsmQJmZmZrFixYu19sitduvTaPxcrVmzt98WKFVs7H6J79+6MHDmSRo0aMWTIkHVWCfxmzZo1VKpUae0KgU3VviVKlCixtrUEWGcLyuw1Fy9enOXLl2/yediYDd3ntzrXP8fWzsrYEFc4SJIkSZJy1PPPP0+3bt349ttvycjI4LvvvqN27dprVyx89NFHzJo1izVr1vDMM8/QokULIL6J/23A4/Dhw9den13dunXJyMhgxowZADz11FO0atUKgHPPPZd+/frRtWtXrrzyym2uf/HixVSrVo3MzMx1BjpWqFCBxYsXA1CxYkVq167Nc889B8Q3859++ikAhxxyyNq5FVs6ELJWrVpMmjQJgEmTJjFr1qxNHl+xYkWqV6/OyJEjAfj1119ZtmzZOjWu79BDD11bz7hx46hSpQoVK1bcovq2hYGDJEmSJClHjRgxYu1S/9+ccMIJa3dPOPjgg+nTpw/169endu3aa48tV64cU6dOpUmTJowdO5brr7/+D49dpkwZnnjiCU466SQaNGhAsWLFOP/88xk6dCglSpTg1FNPpU+fPkyYMIGxY8duU/39+vXjoIMOom3bttStW3ft9SeffDJ33XUX+++/PzNnzmTYsGEMHjyYRo0asd9++/Hyyy8DcP/99/Pggw/SrFmztS0Lm3PCCSewYMECGjduzMCBA9e2jGzKU089xYABA2jYsCF/+ctf+Omnn2jYsCElSpSgUaNG9O/ff53j+/bty8SJE2nYsCF9+vThySef3IpnZeuFbVmGkdeaNm2aTJw4Me0yJEmSJKlA+PLLL6lXr17aZWzQuHHj1hlqmF358uVZsmRJClVpS2zov6sQwsdJkjTd0PG5tsIhhPB4COG/IYQpG7jtshBCEkKoklvnlyRJkiRJ6cnNloohQLv1rwwh7A60BWbn4rklSVJ+NG8eTJgQL5W7fK4l5VOtW7fe4OoGwNUNhUyuBQ5JkrwDLNjATf2BK4D838shSZJyzogRULMmtG0bL7Ntj6Yc5nMtScoH8nRoZAihI/B9kiSf5uV5JUlSyubNg7PPhuXLYeHCeHn22X76nht8riVJ+USJvDpRCGEH4BrgyC08/lzgXIAaNWrkYmWSJCk3LV4M741YwLg1tzCOg/mRalTjR3bNnMeuF6xh18aw666w227xctddoVIl2Mrty/WbjAwoVSoGDb8pWTJeX7VqWlVJkoqgPAscgD2B2sCnIb6CqA5MCiEcmCTJT+sfnCTJo8CjEHepyMM6JUnSdli8GN57D8aNi18TJ8Lq1XUowR4cxIe04m1+Yhemr96DcaOr8r/n/vgYZcr8Hj6sH0Zk/ypfPq9/ugKgVi1YuXLd6zIz4/WSJOWhPAsckiT5HNj5t+9DCBlA0yRJfs6rGiRJUs7bcMAQP1Q/8EDo0wdat4aDv3uJchd0jzdkZsLgwXDKfixfDj/+CD/8AN9/Hy+zf02eDKNGwdKlfzx3xYobDiKyhxTVqkHp0nn7nKSqatX43J599rrPtasbJOWxOXPmcMEFF/DFF1+wZs0aOnTowF133UWpUqX+cOwPP/zARRddxPPPP7/JxzzmmGMYPnw4lSpV2up6+vbtS/ny5bnsssu2+r7ZTZw4kaFDhzJgwICNHvPLL78wfPhwevXqtV3nKuhCkuTO4oEQwgigNVAFmAvckCTJ4Gy3Z7CFgUPTpk2TiRMn5kqdkiRp62wuYGjdOitgOBjKlVvvzvPmxaX9tWpt9RvgRYv+GEZsKKRY/8N9gMqVN79a4s9/hhJ5ufYzt23Hcy2p4Pvyyy+pV69eaudPkoSDDjqInj17cuaZZ7J69WrOPfdcdtppJ+666651jl21ahUl8uAXcE4FDlsiIyODDh06MGXKlFw/V17a0H9XIYSPkyRpuqHjc+1vNUmSUzZze63cOrckSco5ixfD+PG/Bwwff7yRFQwbChjWV7XqNr/5rVgxftWtu/FjkgQWLNj4aokffoApU+Cnn+LPkF2xYjF02NyKicqV47H53nY815KKqBwMKseOHUuZMmU488wzAShevDj9+/endu3a3HjjjTz77LOMGjWKFStWsHTpUh5//PG1b9CXLVtG9+7dmTZtGvXq1SMjI4MHH3yQpk2bUqtWLSZOnMiSJUs4+uijadGiBe+//z677bYbL7/8MmXLlmXQoEE8+uijrFy5kr322ounnnqKHXbYYaO1du/enTJlyjB16lTmzp3LvffeS4cOHVixYgU9e/Zk4sSJlChRgnvvvZc2bdowbtw47r77bl599VX69u3L7Nmz+eabb5g9eza9e/fmoosuok+fPsycOZPGjRvTtm1bLrnkErp06cKiRYtYtWoVAwcOpGXLltv1HBcEhSnHlyRJOSBHA4Y8FkIMBCpXhgYNNn7c6tXw3/9ufLXE7Nnwn/9seGOHkiVjm8bmVkzsuKODLyUVICNGxFasUqXiUrHBg+GUTX6GvElTp06lSZMm61xXsWJFatSowYwZMwD44IMP+Oyzz9hpp53IyMhYe9xDDz3En/70Jz777DOmTJlC48aNN3iOr7/+mhEjRjBo0CA6d+7MCy+8wGmnncbxxx9Pjx49ALj22msZPHgwF1544SbrzcjI4O2332bmzJm0adOGGTNm8OCDDwLw+eefM23aNI488kimT5/+h/tOmzaNt956i8WLF1OnTh169uzJ7bffzpQpU5g8eTIA99xzD0cddRTXXHMNq1evZtmyZVv2RBZwBg6SJBVxmwoYDjoIrrrq94BhEx8QFSjFi8fQoFo1WO/18DpWroyrITa2WmLaNBgzJu4+ub4ddtj8aolddy08z6mkAiz7drq/7XBz9tlwxBHbvNIhSRLCBlLX7Ne3bduWnXba6Q/HjB8/nosvvhiA+vXr07Bhww2eo3bt2mvDiCZNmqwNLaZMmcK1117LL7/8wpIlSzjqqKM2W2/nzp0pVqwYe++9N3vssQfTpk1j/Pjxa4OKunXrUrNmzQ0GDu3bt6d06dKULl2anXfemblz5/7hmGbNmnHWWWeRmZnJcccdt9EQpbAxcJAkqYhZtGjdGQxFIWDYVqVKQY0a8WtTli79ffDlhmZLTJwYv8++U+Vvdtxx86slqlWLtUhSrsiF7XT3228/XnjhhXWuW7RoEd999x177rknH3/8MeU2skxuS+cMls42Ebh48eIsz6q/e/fujBw5kkaNGjFkyBDGjRu32cdaPxwJIWxzHatWrfrDMYceeijvvPMOo0aN4vTTT+fyyy+nW7duW/T4BZmBgyRJhZwBQ+4rVw722it+bUySxL+Lja2W+P77+Pfz449xY4n1Va26+dUSO+8cV29I0lbJhe10Dz/8cPr06cPQoUPp1q0bq1ev5tJLL6V79+6bnKcA0KJFC5599lnatGnDF198weeff75V5168eDHVqlUjMzOTYcOGsdtuu232Ps899xxnnHEGs2bN4ptvvqFOnToceuihDBs2jMMOO4zp06cze/Zs6tSpwwcffLDZx6tQoQKLFy9e+/23337LbrvtRo8ePVi6dCmTJk0ycJAkSQXPokV/bJFYs8aAIW0hxNUMO+4I++678ePWrIH58ze9E8fkyTB3bjw2u2LFYJddfg8gateG7t2hiKzclbStcmE73RACL730Er169aJfv36sWbOGY445hltvvXWz9+3VqxdnnHEGDRs2ZP/996dhw4bsuOOOW3zufv36cdBBB1GzZk0aNGiwzhv/jalTpw6tWrVi7ty5PPzww5QpU4ZevXpx/vnn06BBA0qUKMGQIUPWWc2wKZUrV+aQQw6hfv36HH300dSvX5+77rqLkiVLUr58eYYOHbrFP09BlmvbYuYkt8WUJGnjNhcwZN+m0oCh8Fi1KoYOG1ot8VtIMWNGXCF9xBFw2WVw5JEOspSKim3aFjOfbKe7evVqMjMzKVOmDDNnzuTwww9n+vTplMql3rLu3bvToUMHTjzxxFx5/MIk32yLKUmScsemAobmzeGaa2LA0Ly5AUNhVqJEbKfY1ErhX36BRx+F+++Hdu3izh2XXhoHzzsTQtIf5JPtdJctW0abNm3IzMwkSRIGDhyYa2GDcpcrHCRJyuc2FzD8toLBgEEbs3IlPP003H03fP55bLe4+GI491yoVCnt6iTlhm1a4SBthiscJEkq4FzBoJxWqhR06wannw5vvhmDhyuvhJtvhh49YviwuZ04JEnaWgYOkiSlbOHCdQOGSZMMGJQ7QoCjjopfkyfH4OH+++NXly5xzsP++6ddpaSckiTJH7Z7lLbVtnRH2FIhSVIe21jAUKrUukMeDRiUF2bPhgED4qyHxYvh8MNj8HDUUQ6YlAqyWbNmUaFCBSpXrmzooO2WJAnz589n8eLF1K5de53bNtVSYeAgSVIu21TAsP4MhrJl061VRdfChb8PmPz+e6hfPwYPDpiUCqbMzEzmzJnDihUr0i5FhUSZMmWoXr06JUuWXOd6AwdJkvKQAYMKspUr4Zln4K67fh8wedFFcN55DpiUJP2RgYMkSbnIgEGFUZLA6NFxzsPo0VC+PJxzDvTuDTVrpl2dJCm/MHCQJCkHLVwI7777e8DwyScGDCrcPv00Bg9PPx2DiM6dY7vFAQekXZkkKW0GDpIkbQcDBin67rs4YPKRR+KAycMOi8FDu3YOmJSkosrAQZKkrfDLL+u2SGQPGA4++PeA4aCDDBhUNC1cCIMGwX33xQGT++0Hl14Kp54KpUunXZ0kKS8ZOEiStAmLFsE77xgwSFtr5Up49tnYbvHpp1Ct2u8DJv/0p7SrkyTlBQMHSZI2YNUqeOghuO66GDoYMEjbJkng3/+OwcObb0K5cr8PmKxVK+3qJEm5aVOBQ7G8LkaSpPzgvfegSRO4+OI4e2HMmNhKMW4c9O0bAwfDBmnLhABt28Ibb8DkyXD88fDgg7DXXnDKKfDxx2lXKElKg4GDJKlImTsXuneHFi1gwQJ4/nl4/fU4/M6AQdp+jRrB0KEwaxZccgm89ho0bRr/H3vttdiuJEkqGgwcJElFwqpV8MADUKcODB8OffrAtGlwwglO15dyQ/XqcOedcWeLu++Gr7+G9u2hQQN4/HH49de0K5Qk5TYDB0lSoffBB9CsGVx4Ybz87DO47bbYZy4pd1WsGHew+OYb+Mc/oGRJOPvsONvhttvgf/9Lu0JJUm4xcJAkFVrz5sU3Nn/5S/zzs8/GgXZ166ZdmVT0lCwJXbvGXWBGj46tF1dfDbvvHodLZmSkXaEkKacZOEiSCp3Vq2HgQNhnn9hLfsUVsX3ipJNsn5DSFgIccUScnfLpp7Gt6cEHYc894eSTwY3JJKnwMHCQJBUqH34IBx4IvXrB/vvH9ok77oDy5dOuTNL6GjaEJ5+MAyYvuwz+9a/Y9tS6NYwa5YBJSSroDBwkSYXCzz9Djx5xi8uffoKnn45bXdarl3ZlkjanevUYDH73HdxzT5z30KED1K8Pgwc7YFKSCioDB0lSgbZ6NTzySGyfGDIkfko6bRp06WL7hFTQVKwYt9KcOROGDYPSpeGcc6BmTbj11riVrSSp4DBwkCQVWBMmxBUN558fl2ZPngx33QUVKqRdmaTtUbIknHoqTJoE//53bI+65po4YPKii2ILhiQp/zNwkCQVOPPnw3nnwUEHwZw58ZPQt96C/fbbxgecNy+mF/Pm5WidkrZPCHD44XG2w2efQefO8PDDsNdecRXThAlpV5iNv0ck6Q8MHCRJBcaaNTBoUGyfGDw4bqX31Vfxk9Btbp8YMSKu127bNl6OGJGjNUvKGQ0awBNPxNUNl18Ob7wRB8S2agWvvprygEl/j0jSBoUkSdKuYbOaNm2aTHSPJEkq0iZOhAsugI8+gpYt4zZ6DRps54POmxffHCxf/vt1ZcvCt99C1arb+eCSctPixfDYY3DffTB7NtStC5deCqedBmXK5GEh/h6RVMSFED5OkqTphm5zhYMkKV9bsAB69oyfZH77LTz1FLz9dg6EDQAZGVCq1LrXlSwZr5eUr1WoAH//O8yYAcOHx/f4PXpArVpwyy2x9SpP+HtEkjbKwEGSlC+tWRPbJurUgUcfjYPivvoqfnqZY7tP1KoFK1eue11mZrxeUoFQsiSccgp8/HHcCveAA+Daa6FGDbjwwrjFZq7y94gkbZSBgyQp35k0CQ45JG6HV6dO/P6++2DHHXP4RFWrxlSjbNm4H1/ZsvF7l0FLBU4IcNhh8Npr8PnncajkI4/A3nvHYZMffZRLJ/b3iCRtlDMcJEn5xv/+B9ddBwMHQpUqcYvL00/PwRUNGzNvXlz+XKuWbxKkQuSHH+D//i/+Tlm4MM5/ufxyaN8eiuX0x27+HpFURG1qhoOBgyQpdWvWwNChcMUVse/6ggvgppugUqW0K5NUGCxeHBcd9O8fB0zWqRMHTJ5+eh4PmJSkQsihkZKkfGvy5Pip45lnwl57xT7sAQMMGyTlnAoV4ja6M2fGHSvLlYNzz42bS/Trl4cDJiWpiDFwkCSl4pdf4iDIJk3g66/hiSdg/Hho3Ao0Gx8AACAASURBVDjtyiQVViVKwMknx212x46Fpk3h+uth993hb3+LgYQkKecYOEiS8lSSxPaJOnXgwQfjlpdffQXdu+dCT7UkbUAI0KYNjBoFU6bEXS4GDYJ99oGTToIPP0y7QkkqHHxpJ0nKM599BoceCmecAbVrw4QJ8MAD8Kc/pV2ZpKJqv/3ifIeMDLjySvj3v6F589jq9fLLccaMJGnbGDhIknLdwoWxf/qAA+DLL+Gxx+D99+P3kpQfVKsGt94K330H998Pc+bAccdBvXrw6KOwfHnaFUpSwWPgIEnKNUkC//hHbJ8YMAB69IDp0+Hss22fkJQ/lS8f58t8/TU8/XQcOHneeb8PmPz557QrlKSCw5d7kqRcMWUKtG4dt52rUQM++ggGDoSddkq7MknavBIloEuX2Pr11ltw4IFxwGSNGnHr3hkz0q5QkvI/AwdJUo5atAguuSTuNjFlSlyK/J//xGnwklTQhBDD01dfhalT4dRTY1vYPvvACSfE32+SpA0zcJAk5Ygkifvb160L990X2yamT49tFLZPSCoM9t03hg3ffgtXXRVXPhx8MLRoASNHOmBSktbnS0BJ0nabOhUOOyx+8rfbbvETv0cegcqV065MknLeLrvALbfA7NlxPs3330OnTtCuHSxYkHZ1kpR/GDhIkrbZ4sVw+eWxfeLTT+Hhh2PYcOCBaVcmSbmvfHm48MI4YHLgQHj7bWjWDD7/PO3KJCl/MHCQJG21JIFnnontE3ffDd27x/aJ886D4sXTrk6S8laJEnD++TFwWL48tlm88ELaVUlS+gwcJElb5csv4Ygj4OST47LiDz6AQYOgSpW0K5OkdDVvDh9/DA0awIknwrXXOtdBUtFm4CBJ2iJLlsCVV0LDhjBpEjz0UNzqsnnztCuTpPyjWjUYNy4Ozr3lFujYERYuTLsqSUqHgYMkaZOSBJ57DurVgzvvhG7dYvtEz562T0jShpQuHVd+PfggvPEGHHQQTJuWdlWSlPcMHCRJG/XVV3DkkdC5c2yZeO89GDwYqlZNuzJJyt9CgF69YMyYuHPFQQfBq6+mXZUk5S0DB0nSHyxdGveYb9AAJkyA//u/ePmXv6RdmSQVLIceGuc67LVXbK+4+WbnOkgqOgwcJElrJUmcrF6vHtx+O5x6alzl8Le/xSnskqStt/vuMH48dO0K110HJ50U5+JIUmFn4CBJAuI+8kcfHSer/+lP8O67MGQI/PnPaVcmSQVf2bIwdCjccw+MHBm3zpw5M+2qJCl3GThIUhG3bFncuq1+/bjF5f33x+W/LVqkXZkkFS4hwCWXxEGSP/wAzZrB6NFpVyVJucfAQZKKqCSJn7LVqxe3buvSJbZPXHSR7ROSlJuOOCLOxaleHdq1g7vvjr+TJamwMXCQpCJoxgxo3x46dYKKFeHtt+NS3112SbsySSoa9tgjrio74QS4/HI47bS44kySChMDB0kqQpYtg+uvh/32iwPM+veHSZPiFHVJUt4qVw6eeQZuvRVGjIitbN9+m3ZVkpRzDBwkqYj45z9j0NCvX5yQ/tVX0Ls3lCyZdmWSVHSFELchfuWVOESyaVMYNy7tqiQpZxg4SFIhN3MmdOgAf/1r/DRt3Dj4xz+gWrW0K5Mk/aZ9e/joI6hSJc54eOAB5zpIKvgMHCSpkFq+HPr2jasa3n47bsX2ySfQqlXalUmSNqROHfjwQzjmGLjwQjjnHFixIu2qJGnbGThIUiH06qsxaLjxRjj++Ng+cckltk9IUn5XsWLcQej66+Hxx6F167iFpiQVRAYOklSIzJoFHTvCscdCmTIwZgwMHw677pp2ZZKkLVWsWAyMX3gBpkyBJk3ijhaSVNAYOEhSIbBiRRwGue++MHYs3HknTJ4Mhx2WdmWSpG11/PHwn//E+TutWsFjj6VdkSRtnVwLHEIIj4cQ/htCmJLturtCCNNCCJ+FEF4KIVTKrfNLUoE3bx5MmBAvN+Ff/4L69ePy244dYdq0uKd7qVJ5VKckKdfUrx//KTjsMOjRAy64AFauTLsqSdoyubnCYQjQbr3rRgP1kyRpCEwHrsrF80tSwTViBNSsCW3bxssRI/5wSEYGdOoUh4uVKAGjR8f93KtXz/tyJUm5509/glGj4Ior4KGH4i4Wc+emXZUkbV6uBQ5JkrwDLFjvujeTJFmV9e1/AF8WS9L65s2Ds8+O20wsXBgvzz577UqHX3+FW26J7RNvvgm33w6ffRZfgEqSCqfixeGOO2L+PHEiNG0aLyUpP0tzhsNZwL82dmMI4dwQwsQQwsR5m1lOLEmFSkbGH/shSpaEjAzeeCMur7322rhn+7RpcOWVtk9IUlFx8snw3ntxsGTLlvDUU2lXJEkbl0rgEEK4BlgFDNvYMUmSPJokSdMkSZpWrVo174qTpLTVqvWHBt3Zv/6ZE25sSLt28UXmG2/Ac8/B7runU6IkKT377x9XNzRvDt26xW2PV63a/P0kKa/leeAQQjgD6AB0TZIkyevzS1K+V7UqDB4MZcvya4Uq3FbiOurxBf8aW5pbb43tE0cemXaRkqQ0Va0a2+ouugj694d27WD+/LSrkqR1lcjLk4UQ2gFXAq2SJFmWl+eWpALllFMYXeJo/nZZaaYvLsvxHeHee+P8SEmSIHbb3X9/XPFw/vlxrsPLL0PDhmlXJklRbm6LOQL4AKgTQpgTQjgbeACoAIwOIUwOITycW+eXpILsjjvgyM6VWFOqLP/6F7zwgmGDJGnDuneHd96J3XgHHxxb7iQpP8i1FQ5JkpyygasH59b5JKkwSBK4/nq4+WY45RR4/HEoUybtqiRJ+d2BB8LHH8MJJ0DnznDVVdCvX9zdQpLSkuYuFZKkbJIkDv66+WY455w4edywQZK0pXbZBd56C849F267DTp2hF9+SbsqSUWZgYMk5QOrV8N558F998HFF8Ojj/qplCRp65UqBY88Ag8/HIdKHnggfPll2lVJKqoMHCQpZatWxW3NBg2Ca66J08ZDSLsqSVJBdt55cbXDokVw0EFxmKQk5TUDB0lK0a+/xl7b4cPh1ltjO4VhgyQpJ7RoARMnQp06cNxxcNNNsGZN2lVJKkoMHCQpJcuWwV//Ci+9BAMGxAFfkiTlpOrV4w4W3brBDTfEoZKLF6ddlaSiwsBBklKweDEcfXTsrx08GC68MO2KJEmFVdmyMGRInBP0yivQvDnMmJF2VZKKAgMHScpjCxbAEUfAe+/FVoqzzkq7IklSYRdCHEr8xhswdy40awavv552VZIKOwMHScpD//0vtGkDkyfDCy/AySenXZEkqSg5/HCYMAFq1oT27eHOO+O2zJKUGwwcJCmPfP89HHoofP01vPpqnN8gSVJeq107rrI76SS48ko49dQ4V0iScpqBgyTlgVmzoGVL+OGHuJy1bdu0K5IkFWXlysGIEXD77fDMM3DIIZCRkXZVkgobAwdJymXTpsWwYeFCGDs2/lmSpLSFEFc4vPZaDMabNo3/TklSTjFwkKRc9OmnsY1i1SoYNy6+mJMkKT9p1y7Oddh5ZzjySLj/fuc6SMoZBg6SlEs+/BBat4bSpeMe6A0apF2RJEkbtvfe8d+tDh2gd28480xYsSLtqiQVdAYOkpQL3n47bn25007w7ruwzz5pVyRJ0qZVqAAvvgh9+8KTT8YVenPmpF2VpILMwEGSctjrr8flqbvvHsOGWrXSrkiSpC1TrBjccAOMHAlffhlbAd97L+2qJBVUBg6SlINeegk6doS6deMqh113TbsiSZK23l//GlssKlSANm3gkUfSrkhSQWTgIEk5ZNiwuKd5kybw1ltQtWraFUmStO323TcOkzziCDj//Pi1cmXaVUkqSAwcJCkHDBoEp58e+11Hj4ZKldKuSJKk7VepErzyCvTpE1c5HHYY/PRT2lVJKigMHCRpO913H5x7Lhx9NIwaBeXLp12RJEk5p3hxuO02eOYZ+OSTONdhwoS0q5JUEBg4SNI2ShK4+Wb4+9/hhBPi/IayZdOuSpKk3NG5M7z/PpQsCS1bxp0sJGlTDBwkaRskCVx9NVx3XWylePppKFUq7aokScpdjRrF1Q2HHALdu8PFF0NmZtpVScqvDBwkaSutWQMXXQS33x4HaA0ZAiVKpF2VJEl5o0oVeOMN6N0bBgyAo46Cn39OuypJ+ZGBgyRthdWr4Zxz4IEH4NJL4aGH4p7lkiQVJSVKQP/+sa3i/ffjXIfJk9OuSlJ+48tkSdpCmZnQtSs88QT07Qt33QUhpF2VJEnp6dYNxo+Pgfxf/hJbDCXpNwYOkrQFVqyIgyGfeSYGDTfcYNggSRLE1Q0TJ0KTJnDKKXDllTGAkCQDB0najKVL4dhj4z7kDz0El12WdkWSJOUvf/4zjBkDPXvCnXdC+/bwv/+lXZWktBk4SNImLFwYh2GNHRuHQ/bsmXZFkiTlT6VKxWD+0Ufjv5sHHghTp6ZdlaQ0GThI0kbMnw+HHw4ffhhbKc44I+2KJEnK/3r0gHHjYMkSaN4cRo5MuyJJaTFwkKQN+OknaN0apkyJL5ROPDHtiiRJKjj+8pc412HffaFTpzj7aM2atKuSlNcMHCRpPbNnQ8uWMGsWvPZa7EOVJElbZ7fd4O23oXt3uOmmGDwsWpR2VZLykoGDJGUzY0YMG+bNgzffhMMOS7siSZIKrjJl4PHHYcAAGDUqtlhMn552VZLyioGDJGWZOjWGDUuXxmFXf/lL2hVJklTwhQAXXgj//ncM9A88MK4glFT4GThIEjBpErRqFV8UvfMOHHBA2hVJklS4tG4d5zrUrg0dOsDtt0OSpF2VpNxk4CCpyHv//dg6Ua5cDBv23TftiiRJKpxq1oT33oOTT4arroIuXeLKQkmFk4GDpCJt7Fg48kjYeWd4913Ya6+0K5IkqXDbYQcYNgzuugteeCG2MM6alXZVknKDgYOkImvUKDjmmLi08513oEaNtCuSJKloCAEuuyzOcpg9G5o2hTFj0q5KUk4zcJBUJD33HBx3HNSvD+PGwS67pF2RJElFz1FHwYQJUK1a/PN99znXQSpMDBwkFTlDh8be0YMOip+mVK6cdkWSJBVde+0FH3wAHTvC3/8OZ5wBy5enXZWknGDgIKlIGTgwvpA57DB44w3Ycce0K5IkSRUqwPPPQ79+8NRTcZvq775LuypJ28vAQVKRcffd0KsXHHssvPJK3JVCkiTlD8WKwbXXwj//CdOnx7kO776bdlWStoeBg6RCL0mgb1+4/PK4/dYLL0CZMmlXJUmSNuTYY+Gjj6BSpbgiceBA5zpIBZWBg6RCLUli0HDjjXDmmXEbrpIl065KkiRtSt268OGHcevqXr3gvPPg11/TrkrS1jJwkFRorVkTX6Tccw/87W/w2GNQvHjaVUmSpC1RqVJsr7jmGhg0CNq0gR9/TLsqSVvDwEFSobRqFXTvDg8/DH36wIABsTdUkiQVHMWLw803x+2sP/ssznX48MO0q5K0pXz5LanQWbkybnv51FPxRcptt0EIaVclSZK21Yknxq0zS5eGQw+N/8ZLyv8MHCQVKsuXw3HHxcGQ/fvHZZiSJKnga9AAJkyAFi3iKsZXXkm7IkmbY+AgqdBYvBiOOQZefx0efRR69067IkmSlJMqV45zHQ44IK5m/PjjtCuStCkGDpIKhf/9L06yfvfduMyyR4+0K5IkSbmhXLm4uqFqVejQAWbPTrsiSRtj4CCpwJs3L+7TPWlSHCrVtWvaFUmSpNy0yy4walRspWzfHhYuTLsiSRti4CCpQPvhB2jVCqZNi0ssO3VKuyJJkpQX9tsvzmyaNg1OOgkyM9OuSNL6DBwkFVgZGdCyJXz3XZzbcNRRaVckSZLy0uGHx7lNo0dDr16QJGlXJCm7EmkXIEnbYvr0+CJjyRIYMwYOPDDtiiRJUhrOPBNmzoRbboE994Q+fdKuSNJvDBwkFTiffw5t28KaNTBuHDRqlHZFkiQpTf36wTffwFVXQe3a0KVL2hVJAlsqJBUwEybEmQ0lSsA77xg2SJIkCAGeeAJatIAzzoD33ku7Iklg4CCpABk/PrZRVKoUt7+sWzftiiRJUn5RujSMHAk1asBf/wozZqRdkSQDB0kFwujRcOSRsOuucWVD7dppVyRJkvKbypXhtdfin485BubPT7ceqagzcJCU7/3zn9ChA+y9dwwbqldPuyJJkpRf7bUXvPwyzJ4dt8v+9de0K5KKLgMHSfna00/D8cdD48bw1luw885pVyRJkvK7Qw6BoUNjC+ZZZ7ldppQWd6mQlG89/jiccw60bAmvvgoVKqRdkSRJKig6d/5954o99og7WUjKWwYOkvKlAQPg4ovhqKPgxRdhhx3SrkiSJBU0V14JM2fCzTfH0OHMM9OuSCpabKmQlO/cdlsMGzp1ij2Yhg2SJGlbhAAPPRQHT597LowZk3ZFUtFi4CAp30gSuOYauPpq6NoVnn02bnElSZK0rUqWjK8p6taNc6GmTk27IqnoMHCQlC8kCfz973DrrdCjBzz5JJSw6UuSJOWAHXeEUaPiqsljjoGffkq7IqloMHCQlLrVq+Myx/vvh9694ZFHoHjxtKuSJEmFSY0acQj1zz/DscfC0qVpVyQVfgYOklKVmQmnnw6PPQbXXQf33hv7LSVJknJakyZxy+1Jk2L75urVaVckFW65FjiEEB4PIfw3hDAl23U7hRBGhxC+zrr8U26dX1L+9+uvcNJJMGIE3H473HSTYYMkSdqMefNgwoR4uQ2OPRbuuy8Opr7sshyuTdI6cnOFwxCg3XrX9QHGJEmyNzAm63tJRdCyZdCxY/zH/v/+L25bJUmStEkjRkDNmtC2bbwcMWKbHubCC+OOWPfdBw88kMM1SlorJEmSew8eQi3g1SRJ6md9/xXQOkmSH0MI1YBxSZLU2dzjNG3aNJk4cWKu1Skpby1aBB06wHvvxVYK98SWJEmbNW9eDBmWL//9urJl4dtvoWrVrX641avhhBPglVfiByAdOuRgrVIREkL4OEmSphu6La9nOPw5SZIfAbIud97YgSGEc0MIE0MIE+dt43IpSfnPggVwxBHwwQcwfLhhgyRJ2kIZGVCq1LrXlSwZr98GxYvDsGFwwAHQpUuc6yApZ+XboZFJkjyaJEnTJEmaVt2GxFJS/jN3LrRuDZ99Bi++GP9xlyRJ2iK1asHKletel5kZr99G5crFFQ5VqsQVDt99t10VSlpPXgcOc7NaKci6/G8en19SSubMgUMPhZkz45ZUxx6bdkWSJKlAqVoVBg+ObRQVK8bLwYO3qZ0iu112gddei9tktm8fWz8l5Yy8Dhz+CZyR9eczgJfz+PySUvDNN9CyJfz0E7zxRmypkCRJ2mqnnBJnNvz73/HylFNy5GH32w9eeAG+/DLuoJWZmSMPKxV5ubkt5gjgA6BOCGFOCOFs4HagbQjha6Bt1veSCrEvv4xhw6JFMHYstGiRdkWSJKlAq1oVmjXb7pUN6zviCHjkEXjzTbjgAsjF2fpSkVEitx44SZKNxY2H59Y5JeUvkyfDkUdCsWLw9ttQv37aFUmSJG3cWWfF9s9bb4U993Tbbml75VrgIKlo+89/4OijoUIFGDMG9t477YokSZI2r1+/2A7apw/Urg2dO6ddkVRw5dtdKiQVXOPGQdu2ULkyvPuuYYMkSSo4ihWDJ56AQw6Bbt3g/ffTrkgquAwcJOWof/0rrmyoUSOGDTVrpl2RJEnS1ilTBkaOhN13h7/+NbZZSNp6Bg6ScsyLL8Z/lOvVizMbqlVLuyJJkqRtU6VK3C4zSeCYY2D+/LQrkgoeAwdJOeIf/4g9jk2bxt0oqlRJuyJJkqTts/fe8PLLcQfOTp3g11/TrkgqWAwcJG23Rx+NPY6tWsWtpCpVSrsiSZKknHHIIfDkk7FV9Kyz3C5T2hruUiFpu/TvD5dcAu3bw3PPQdmyaVckSZKUs7p0iTtXXH113C7zppvSrkgqGAwcJG2TJIGbb4brr4eTTootFaVKpV2VJElS7ujTJw6P7NcP9tgDundPuyIp/zNwkLTVkiT+o3vnnXDGGfDYY1DC3yaSJKkQCwEGDoTZs6FHj7gj12GHpV2VlL85w0HSVlmzBi68MIYNPXvC448bNkiSpKKhZMnYQlq3Lhx/PHzxRdoVSfmbgYOkLbZqFZx9Njz4IFx+ebws5m8RSZJUhOy4I4waFedWHXMM/PRT2hVJ+ZdvFSRtsYsvhiFD4MYb4Y474tJCSZKkoqZGDXj1VZg3Dzp2hGXL0q5Iyp8MHCRtkWHD4KGH4LLL4qBIwwZJklSUNWkCI0bAxInQtSusXp12RVL+Y+AgabO++ALOPRcOPRRuuy3taiRJkvKHjh3hvvtg5MjYbippXY56k7RJS5bACSdAhQrw9NMOiJQkScruoovidpn9+8Oee8IFF6RdkZR/+NZB0kYlSdz2afp0+Pe/oVq1tCuSJEnKf+69FzIyYvhQsyZ06JB2RVL+YEuFpI0aODCuarj5ZmjTJu1qJEmS8qfixWH4cNh/fzj5ZJg0Ke2KpPzBwEHSBn30EfTuHRP6K69MuxpJkqT8rVw5eOUVqFw5vn767ru0K5LSZ+Ag6Q/mz4eTToLddoMnn4Ri/qaQJEnarGrVYNQoWLoU2reHRYvSrkhKl28jJK1jzRro1g1++gmeew522intiiRJkgqO+vXh+efhyy+hc2fIzEy7Iik9Bg6S1nHbbfDaa3D//dC0adrVSJIkFTxt28LDD8Mbb8Df/hYHcUtFkbtUSFprzBi4/nro2hXOOy/taiRJkgqus8+O22XedlvcLvOKK9KuSMp7Bg6SAPj+ezj1VKhbNybyIaRdkSRJUsF2883wzTdxAHft2nFGllSUGDhIIjMzbuG0dCmMGwfly6ddkSRJUsFXrBgMGQJz5sDpp0P16nDwwWlXJeUdZzhI4uqrYfx4eOwxqFcv7WokSZIKjzJlYORI2H136NgxtllIRYWBg1TEvfQS3H03XHBBXOUgSZKknFWlShzKnSRxu8wFC9KuSMobBg5SETZzJnTvDs2awT33pF2NJElS4bX33nGlw6xZ0KkT/Ppr2hVJuc/AQSqili+HE0+E4sXhueegdOm0K5IkSSrcWrSAJ5+Ed96Bc85xu0wVfg6NlIqoiy6CyZNh1CioWXML7zRvHmRkQK1aULVqLlYnSZJUOJ18cty54pprYI894MYb065Iyj2ucJCKoCFD4oDIa66BY47ZwjuNGBGTibZt4+WIEblZoiRJUqF11VVw9tlw001xxYNUWIWkAKzjadq0aTJx4sS0y5AKhc8+g+bN45ZMb74ZWyo2a968GDIsX/77dWXLwrffutJBkiRpG2Rmxg9+3n4b3ngD2rRJuyJp24QQPk6SpOmGbnOFg1SELFoU5zZUqgTDh29h2ACxjaJUqXWvK1kyXi9JkqStVrIkPP887LNPHCL55ZdpVyTlPAMHqYhIkrh075tv4Jln4M9/3oo716oFK1eue11mZrxekiRJ22THHeM8rbJl42qHuXPTrkjKWQYOUhExYEBM0W+/HVq23Mo7V60KgwfHfw0rVoyXgwfbTiFJkrSdataEV16B//4XOnaEZcvSrkjKOc5wkIqA99+HVq2gQwd48UUIYRsfyF0qJEmScsXLL8fWik6d4pblxfxoWAWEMxykImzePOjcOabnTzyxHWEDxJChWTPDBkmSpBz2179C//7xw6Errki7GilnlEi7AEm5Z/Vq6NoVfv4Z/vOfOCxSkiRJ+dPFF8PMmXDPPbDHHtCrV9oVSdvHwEEqxPr1g9Gj4bHHoHHjtKuRJEnS5vTvHztYL7wwrlBt3z7tiqRtZ0uFVEi98QbcdBN07w5nnZV2NZIkSdoSxYvH7csbN4YuXeCTT9KuSNp2Bg5SIfTdd7GVon59ePDB7ZzbIEmSpDxVvjy8+irstFMc+j1nTtoVSdvGwEEqZFaujEMiV66M22DusEPaFUmSJGlrVasGo0bB4sWxrWLRorQrkraegYNUyFxxRRwQ+fjjsM8+aVcjSZKkbdWgQfwAaerU2F6xalXaFUlbx8BBKkSeew7uvx9694YTT0y7GkmSJG2vI4+Ehx+G11+Hv/0NkiTtiqQt5y4VUiHx1VdxOOTBB8Mdd6RdjSRJknLKOefE7TJvvx323BMuvzztiqQtY+AgFQLLlsUVDWXKwDPPQKlSaVckSZKknHTLLTBrVmyfrV3b1awqGAwcpAIuSaBnz9jb9/rrsPvuaVckSZKknFasGAwZEnesOP10qF4dmjdPuypp05zhIBVwgwfD0KFwww2xx0+SJEmFU5kyMHJkDBs6doRvvkm7ImnTDBykAuyTT+LwoCOPhGuvTbsaSfr/9u482u7p/v/4c2cUIsZbJEUUNVfpVaWouabyraEaX0MJaTOQiCkSJIYQQkQSEZFrrqjEt1TR4tv40dZ0QyqIWmjElLqKmDJn//7Y8S2aiOSec/YZno+1sk7u5xz3vFZ7ctc9r7P3e0uSim3tteG++2DhQjjgAHjvvdyJpKWzcJAq1AcfpL17dXVw663QsmXuRJIkSSqFTTeFu+9OMx0OPRTmzs2dSFoyCwepAsUIxx0HM2bAHXek0kGSJEm1Y5dd0kyH//f/4KSTPC5T5cmhkVIFuvxy+N3vYPjwdAymJEmSak+XLmmOwznnwLe+BYMG5U4kfZGFg1RhHnkEzj4bjjgCTjkldxpJkiTl1L9/Kh3OPz+VDscemzuR9G8WDlIFmTkTjjwSNt4Yxo2DEHInkiRJUk4hwJgx8NprcOKJsMEGsPvuuVNJiTMcpAqxYAEcdRTMmgUTJ0KHDrkTSZIkqRy0bp1+P9x0U/jpT2HatNyJpMTCQaoQAwfCpEmpwd5mm9xpJEmSVE5WXz0dl9m2LRx4ILzzTu5EkoWDVBHuvRcuvjhNIHZfniRJkpZkLTtd7gAAIABJREFUww3hnnvSNtyDD4bZs3MnUq2zcJDK3PTpcMwxsN12MGJE7jSSJEkqZzvsAOPHw5NPpt8hFy3KnUi1zMJBKmNz56bTKBYtggkTYKWVcieSJElSuTvkEBg2DO68E846K3ca1TJPqZDKWN++0NgId92VTqaQJEmSvo7eveGVV+Dyy9Pvkb/6Ve5EqkUWDlKZuu02GD0azjgjtdSSJEnS1xUCDB+etuf27JnmO+y/f+5UqjVuqZDK0AsvQLdusOuuMHhw7jSSJEmqRC1bpnkO224LP/sZTJmSO5FqjYWDVGY+/hgOPxxWWQVuvz2dqyxJkiStiPbt4fe/hzXWSMdlvvFG7kSqJRYOUhmJMa1s+PvfU9nQsWPuRJIkSap0HTumY9Y/+ggOOijdSqVg4SCVkWuuScveLrwQ9tgjdxpJkiRVi222gYkT4bnn0vaKBQtyJ1ItsHCQysSTT0KfPnDAAdCvX+40kiRJqjb77ps+4PrDH+Dkk9PqWqmYshQOIYRTQwjPhxCeCyGMDyGslCOHVC7+9S844oi03O2WW6CFVaAkSZKK4KST0odbY8bAFVcs48FNTfDUU+lWWgFLfVsTQrgvhNC50E8YQugEnALUxxi3BloCPy/080iVYtEiOPZYePttmDAB1lwzdyJJkiRVs8GD4cgj0/Hrd965lAeNH5/O0txnn3Q7fnxJM6o6fNXnqDcCD4QQBoQQCj0nvxXQLoTQClgZeKvA31+qGEOGwH33pXOSd9ghdxpJkiRVuxYt4MYbYeed4eij4fHHv/SApibo2hVmz4ZZs9Jt166udNByW2rhEGO8A9gO6AA0hhBODyH0/ezPij5hjPFN4HJgBvA2MCvG+MCXHxdC6BZCaAwhNDb5wlaV+tOf4NxzoUsX6N49dxpJkiTVipVWgrvvhk6d4OCD4dVXP3fn9OnQps0X/4PWrdN1aTksa6f4fOAToC2w6pf+rJAQwhrAIcBGQEdglRDC0V9+XIxxbIyxPsZYX1dXt6JPJ5Wtt95KRcNmm8HYsRBC7kSSJEmqJWuvnVbaLlwIBx4I77+/+I7OnWHevC8+eP78dF1aDl81w2E/YAppy8P2McaBMcbzP/vTjOfcG/hHjLEpxjgf+B9g52Z8P6nizJ+f9s198knaN9e+fe5EkiRJqkXf/jbcdVda4XDooYt7hro6aGiAdu2gQ4d029CQrkvL4atWOAwAjogx9osxflrA55wB/CCEsHIIIQB7AdMK+P2lsjdgAPz5z2llwxZb5E4jSZKkWrbrrnDDDfDww3DiiYuPy+zSBV57DR56KN126ZI7pipQq6XdEWPctRhPGGN8IoQwEXgaWAA8A4wtxnNJ5eiuu2DoUOjRA446KncaSZIkKf1e+uqrab7YxhvDwIGkFQ2ualAzhBhj7gzLVF9fHxsbG3PHkJrtlVfge9+DTTdNKxzats2dSJIkSUpihBNOSCdY3HwzHHNM7kSqBCGEyTHG+iXdt9QVDpIKa/ZsOPzwdAzRhAmWDZIkSSovIcC118KMGekUzPXXh913z51KlWxZp1RIKpDevWHKFLjlFgf8SpIkqTy1aZOGmm+6Kfz0p/Dii7kTqZJZOEglcNNNcN110L9/OnJIkiRJKlerrw733ptW5B5wwOeOy5SWk4WDVGRTp0L37rDHHnB+cw6UlSRJkkqkc2e4++60vaJv39xpVKksHKQi+vBDOOyw1BLfdhu0cmqKJEmSKsSOO0K/fmmI5P33506jSmThIBVJjGnYzquvwu23w7rr5k4kSZIkLZ9zz4WttoKTToJZs3KnUaWxcJCKZMQImDgRLrkEdtstdxpJkiRp+bVtCzfcAG+/DaedljuNKo2Fg1QEjz0Gp58OhxySbiVJkqRKtcMOcOaZ0NAADzyQO40qSYgx5s6wTPX19bGxsTF3DOlraWqC7beH1q3h6afT/AZJkiSpks2Zk37H/fhjeO456NAhdyKVixDC5Bhj/ZLuc4WDVEALF8LRR6fSYeJEywZJkiRVh5VWguuvhzffTKsdpK/DwkEqoIsuSsvMRo5MDbAkSZJULX7wg3RE5rXXwv/+b+40qgRuqZAK5IEHYL/94Jhj0tFBIeROJEmSJBXW7Nnw3e/C3LkwdSqsumruRMrNLRVSkb3+Ohx1VDoyaPRoywZJkiRVp3bt0qkVM2ZAv36506jcWThIzTRvHhx5ZLqdOBFWWSV3IkmSJKl4dt4Z+vRJH7Q9/HDuNCpnFg5SM511VjoGs6EBNtssdxpJkiSp+C66CDbZBLp2hU8+yZ1G5crCQWqGiRNh+HDo3RuOOCJ3GkmSJKk0Vl45nVrxj39A//6506hcWThIK+ill+CEE9K03ssuy51GkiRJKq1dd4VevWDECHj00dxpVI4sHKQV8OmncPjh0KYN3HFHupUkSZJqzSWXwLe+lT6I+/TT3GlUbiwcpOUUI3TvDs89B7/+Nay/fu5EkiRJUh6rrJJmmb38MpxzTu40KjcWDtJyamiAm2+G886DH/84dxpJkiQpr913hx490myzv/41dxqVkxBjzJ1hmerr62NjY2PuGBLPPAM77QS77Qb33w8tW+ZOJEmSJOX38cew9dbQti1MmQLt2uVOpFIJIUyOMdYv6T5XOEhf0wcfpLkNa6+dtlJYNkiSJElJ+/ZpJfBLL8HAgbnTqFxYOEhfQ4zwi1/AjBlpSGRdXe5EkiRJUnnZay/o1g2uuAKeeCJ3GpUDCwfpa7jiCrj7bhg6FHbeOXcaSZIkqTwNHQqdOsHxx8OcObnTKDcLB2kZHn0U+vWDww6D3r1zp5EkSZLKV4cOcN11MG0anH9+7jTKzcJB+gr//CcceWQ6W/j66yGE3IkkSZKk8vbjH0PXrnDZZfDUU7nTKCcLB2kpFi6ELl3g/fdh4sTU1kqSJElatiuugI4d09aKuXNzp1EuFg7SUgwcCJMmwTXXwHe+kzuNJEmSVDlWWw3GjoXnn4cLL8ydRrlYOKjyNTWltVpNTQX7lvfeC4MHp6Vgv/hFwb6tJEmSVDP23z/9Lj1kCDz9dO40ysHCQZVt/HjYcEPYZ590O358s7/la6/BMcfAttvCyJEFyChJkiTVqGHD4BvfSMXDvHm506jULBxUuZqa0hKE2bNh1qx027Vrs1Y6zJ0LRxyR5jdMnAjt2hUwryRJklRj1lgDrr0Wpk6Fiy/OnUalZuGgyjV9OrRp88VrrVun6yuob9+0O+PGG2GTTZoTTpIkSRLAT34CRx+dtixPmZI7jUrJwkGVq3Pn/1yXNX9+ur4CbrsNRo+G00+Hn/602ekkSZIkLXbVVbDWWunUivnzc6dRqVg4qHLV1UFDQ9r30KFDum1oSNeX0wsvQLdusMsuLvWSJEmSCm3NNWHMmLTCYciQ3GlUKiHGmDvDMtXX18fGxsbcMVSumprSNorOnVeobPj4Y/j+9+Ff/4JnnknnBUuSJEkqvKOOSrPSJk+GbbbJnUaFEEKYHGOsX9J9rnBQ5aurgx12WKGyIca0suHvf09bKiwbJEmSpOIZMSINkjz+eFiwIHcaFZuFg2ramDHpJM0LLoC99sqdRpIkSapua6+d5qZNngxDh+ZOo2KzcFDNeuop6NMH9t8fzj47dxpJkiSpNhx2WDqKftAgeP753GlUTBYOqknvvZd+yK27LtxyC7TwX4IkSZJUMqNGpbnvJ5zg1opq5tss1ZxFi+DYY+Gtt2DChHQ8jyRJkqTS+cY3Uunw5JMwbFjuNCoWCwfVnEsvhXvvhSuvTKdTSJIkSSq9n/0MDj0UzjsPXnwxdxoVg4WDasqkSXDOOfDzn0OPHrnTSJIkSbUrhDRAcpVV0qkVCxfmTqRCs3BQzXjrrVQ0fPvbMHZs+gEnSZIkKZ911oGRI+Hxx+Gqq3KnUaFZOKgmLFiQyoaPP4aJE2HVVXMnkiRJkgTQpQscfDAMGAAvvZQ7jQrJwkE1YcAAePTRtLJhq61yp5EkSZL0mRBgzBhYaaV0aoVbK6qHhYOq3t13w2WXwa9+Bf/937nTSJIkSfqy9dZLWyr+8pd0eoWqg4WDqtorr8Bxx8H3vpdOpZAkSZJUno45Bg48EM4+G15+OXcaFYKFg6rW7Nlw+OFpidaECWmJliRJkqTyFAJcey20aQNdu8KiRbkTqbksHFS1eveGKVPg5ptho41yp5EkSZK0LJ06pZXJjzySjsxUZbNwUFW66Sa47jro1w9+8pPcaSRJkiR9Xb/4Bfz4x+l3+VdfzZ1GzWHhoKozdSp07w677w4XXpg7jSRJkqTlEUL68LBFCzjxRLdWVDILB1WVDz+Eww6D1VaD8eOhVavciSRJkiQtr/XXh2HDYNKkdLS9KpOFg6pGjKkBffVV+M1vYN11cyeSJEmStKK6doV99oEzzoDXXsudRivCwkFVY+TIdBrFxRfDbrvlTiNJkiSpOT7bWgHpg8UY8+bR8rNwUFV47DE47TQ4+GA4/fTcaSRJkiQVwoYbwtCh8NBDMG5c7jRaXhYOqnjvvgs/+1na53XjjWm4jCRJkqTq0K0b7Lln+oBxxozcabQ8fGumihZjOjanqQkmToQ11sidSJIkSVIhtWiRVjcsWpTKB7dWVA4LB1W0O++Ee++FIUNg++1zp5EkSZJUDBttBJdeCn/8Y1rVrMoQYgXUQ/X19bGxsTF3DJWZjz+GLbaAtdaCxkaPwJQkSZKq2aJFaWvFlCnw/PPQqVPuRAIIIUyOMdYv6T5XOKhiXXghvPEGjB5t2SBJkiRVuxYtoKEB5s2DX/7SrRWVwMJBFWnaNBg2DI4/HnbeOXcaSZIkSaWw8cZwySVpW/Utt+ROo2WxcFDFiRF69YL27dPsBkmSJEm14+STYZddoHdvePvt3Gn0VSwcVHF+8xv405/g4ovhG9/InUaSJElSKbVoAddfD3PmwK9+5daKcmbhoIry0UfQt286kaJbt9xpJEmSJOWw6aYweDD87ncwfnzuNFoaCwdVlEGDYObMNCiyZcvcaSRJkiTl0rs37LRT2mIxc2buNFqSLIVDCGH1EMLEEMKLIYRpIYSdcuRQZXnuObjqKjjxRNhxx9xpJEmSJOXUsmXaWvHJJ9Cjh1srylGuFQ5XAX+IMW4ObAtMy5RDFSJG6NkTVlstzW6QJEmSpM03hwsugN/+Fu64I3cafVnJC4cQQgdgN6ABIMY4L8b4QalzqLL8+tfwyCPpVIq1186dRpIkSVK56NsXvv/99AHlO+/kTqPPy7HC4VtAE3BDCOGZEMK4EMIqX35QCKFbCKExhNDY1NRU+pQqG7Nmwemnpx8iXbvmTiNJkiSpnLRqBTfckAbM9+qVO40+L0fh0ArYHrgmxrgd8AnQ78sPijGOjTHWxxjr6+rqSp1RZeS881JTOXp0OgJHkiRJkj5vyy3TgPkJE2DixNxp9Jkcb9/eAN6IMT6x+OuJpAJC+g9/+xuMGpXO1/3e93KnkSRJklSuzjgjvWfo0QPefTd3GkGGwiHGOBN4PYSw2eJLewEvlDqHyt+iRemHxZprwkUX5U4jSZIkqZx9trXigw/SUZnKL9cC9ZOBX4cQngW+C3jugP7DzTfDX/8Kl12WSgdJkiRJ+irbbAPnngu3355OrlBeIVbAYaX19fWxsbExdwyV0Pvvw2abwSabwJ//7OwGSZIkSV/P/Plp4Pzbb8Pzz8Naa+VOVN1CCJNjjPVLus+3cSpL55wD//qXgyIlSZIkLZ/WreHGG9P7iT59cqepbb6VU9mZPBmuuSado/vd7+ZOI0mSJKnSbLstDBgAt94K99yTO03tckuFysqiRbDzzjB9Orz4Iqy+eu5EkiRJkirRvHmwww7Q1JS2VqyxRu5E1cktFaoY118PTzwBQ4daNkiSJElacW3apFMr3nkHTj01d5raZOGgsvGvf0G/frDrrnD00bnTSJIkSap0228PZ58NN90E992XO03tsXBQ2ejfP52Ze/XVEELuNJIkSZKqwTnnwFZbQbduMGtW7jS1xcJBZeHJJ+G66+CUU9LZuZIkSZJUCG3bplMrZs6E007Lnaa2WDgou4ULoUcPWHddGDQodxpJkiRJ1aa+Hs44Axoa4I9/zJ2mdlg4KLvrrktHYV5xBXTokDuNJEmSpGo0cCBssQWcdBJ8+GHuNLXBwkFZNTWl2Q177AE//3nuNJIkSZKq1UorpVMr3nwzrXZQ8Vk4KKt+/eCjj2DUKAdFSpIkSSquHXdMcxzGjoWHHsqdpvpZOCibv/4Vrr8+nYm75Za500iSJEmqBeefD9/+Npx4YvrwU8Vj4aAsFiyAnj2hUyc477zcaSRJkiTVinbt0taKGTPSimsVj4WDsrjmGpgyBa68Etq3z51GkiRJUi3ZeWfo0wdGj4ZJk3KnqV4hxpg7wzLV19fHxsbG3DFUIP/8Z1rCtOOO6UgaZzdIkiRJKrVPP4Vtt4WFC+HZZ/0gdEWFECbHGOuXdJ8rHFRyZ54Js2fDyJGWDZIkSZLyWHnlNFNu+vR0cp4Kz8JBJfXoo3DzzekYms02y51GkiRJUi3bdVc4+eT0Yegjj+ROU33cUqGSmT8ftt8ePvwQXngBVlkldyJJkiRJte6TT+A730mrr599Nq180NfnlgqVhVGj4LnnYPhwywZJkiRJ5WGVVaChAV55BQYMyJ2mulg4qCTeegsGDoT994f/+q/caSRJkiTp33bfHXr0gKuugr/8JXea6mHhoJI44wyYNw9GjHBQpCRJkqTyc+mlsOGGcMIJaci9ms/CQUU3aRLcdhucdRZssknuNJIkSZL0n9q3h3Hj4KWX4LzzcqepDhYOKqr586FnT9hoI+jXL3caSZIkSVq6vfaCX/4Shg2Dxx/PnabyWTioqIYPh2nT0laKdu1yp5EkSZKkr3bZZdCpExx/PMyZkztNZbNwUNG88Qacfz785Cdw0EG500iSJEnSsnXokLZWvPgiDBqUO01ls3BQ0fTtCwsXpkmvkiRJklQp9t0XTjwRhg6Fp57KnaZyWTioKB58ECZMgP790/wGSZIkSaokl18OHTvCL34Bc+fmTlOZLBxUcHPnQq9esPHG6ThMSZIkSao0q60GY8fCCy/AhRfmTlOZLBxUcMOGpaNkRo2ClVbKnUaSJEmSVsz++6cVDkOGwOTJudNUHgsHFdSMGan9++lPYb/9cqeRJEmSpOYZNgzWWSedWjFvXu40lcXCQQXVp0+6HT48bw5JkiRJKoQ11oBrr4WpU2Hw4NxpKouFgwrm/vvht7+Fc8+FDTbInUaSJEmSCuOgg+CYY+Dii2HKlNxpKkeIMebOsEz19fWxsbExdwx9hTlzYOutoWVLePZZaNs2dyJJkiRJKpz33oOttoJ114Unn4TWrXMnKg8hhMkxxvol3ecKBxXE0KHwyitpUKRlgyRJkqRqs+aaMGZMWuEwZEjuNJXBwkHN9o9/pKVFRxwB++yTO40kSZIkFcchh0CXLmlQ/rPP5k5T/iwc1Gy9e6etFMOG5U4iSZIkScU1YkQaJHn88TB/fu405c3CQc1yzz3pz8CB8M1v5k4jSZIkScW19towejQ8/XTaWq6lc2ikVtjs2WloSrt2aR+TQ1MkSZIk1Yojj4S77krFw1Zb5U6Tj0MjVRRDhqT5DaNGWTZIkiRJqi2jRkGHDmlrxYIFudOUJwsHrZCXX4ZLL00DU/bYI3caSZIkSSqtujq4+mp46im44orcacqThYOWW4xwyinQpg1cfnnuNJIkSZKUxxFHwKGHppl206blTlN+LBy03O6+G+6/H84/Hzp2zJ1GkiRJkvIIIQ2QXGUVOOEEWLgwd6LyYuGg5fLJJ+kYzG22gZNPzp1GkiRJkvJaZx0YORIefxyGD8+dprxYOGi5XHwxzJiR9iq1apU7jSRJkiTl16ULHHIInHMOvPRS7jTlw8JBX9vf/57OmT32WNh119xpJEmSJKk8hADXXAPt2rm14vMsHPS1xJi2ULRrB5ddljuNJEmSJJWX9daDq66Cv/wlbbGQhYO+pokT4cEH4aKL0h4lSZIkSdIXHX00HHgg9O8PL7+cO01+Fg5apo8/hlNPhe9+F7p3z51GkiRJkspTCHDttdCmDXTtCosW5U6Ul4WDlunCC+HNN9NxLw6KlCRJkqSl69QpnVbxyCPpPVQts3DQV3rhBRg2LA0+2Wmn3GkkSZIkqfwddxzsvz+cdRa8+mruNPlYOGipYoRevWDVVWHIkNxpJEmSJKkyfLa1olWr2t5aYeGgpbr9dpg0CS6+GOrqcqeRJEmSpMqx/vpwxRXw8MOpfKhFIcaYO8My1dfXx8bGxtwxasqHH8Lmm0PHjvDEE9CyZe5EkiRJklRZYoQf/xgeewymToXOnXMnKrwQwuQYY/2S7nOFg5Zo0CCYOTMNObFskCRJkqTlFwKMG5f+ftJJqYCoJRYO+g9Tp8KIEekfxPe/nzuNJEmSJFWuDTaAyy+Hhx76d/lQKywc9AUxQs+esPrqaXaDJEmSJKl5unWDPfeE006DGTNypykdCwd9wa23wqOPplMp1lordxpJkiRJqnyfba1YtKi2tlZYOOj/fPABnH467LgjnHBC7jSSJEmSVD022gguvRQeeABuuCF3mtKwcND/Oe88aGpKgyJb+MqQJEmSpILq3h1+9CM49VR4443caYrPt5UCYMoUuPrq9A9g++1zp5EkSZKk6tOiBTQ0wIIF8MtfVv/WCgsHsWgR9OiRZjZcdFHuNJIkSZJUvTbeGC65BO67D26+OXea4rJwEDfdBI89BpddBmuskTuNJEmSJFW3Xr1gl12gTx94663caYrHwqHGvf8+nHkm7LwzHHts7jSSJEmSVP1atIDrr4c5c+BXv6rerRUWDjVuwAB47z0HRUqSJElSKW26KQweDPfcA7fdljtNcfgWs4ZNngxjxqTlPNtumzuNJEmSJNWW3r1hp53g5JNh5szcaQovW+EQQmgZQngmhPD7XBlq2WeDIr/xDbjggtxpJEmSJKn2tGyZtlZ8+il0P2Eu8cmnoKkpd6yCybnCoTcwLePz17SGBnjySbj8clhttdxpJEmSJKk2bb45XHjoM9x1f1t+s/s1sOGGMH587lgFkaVwCCF8EzgQGJfj+Wvdu+9Cv36w227w3/+dO40kSZIk1bCmJvr+dleO5HY6zn4ZZs+Grl2rYqVDq0zPOxw4E1h1aQ8IIXQDugFssMEGJYpVG/r3h1mz4OqrIYTcaSRJkiSphk2fTsu2rbh9Tpd/X2vdGqZPh7q6bLEKoeQrHEIIBwHvxBgnf9XjYoxjY4z1Mcb6ugr/H7mcPPEEjBuXhpNsvXXuNJIkSZJU4zp3hnnzvnht/vx0vcLl2FLxQ+DgEMJ04HZgzxDCrRly1JyFC9OgyPXWg0GDcqeRJEmSJFFXl4bstWsHHTqk24aGil/dABm2VMQYzwbOBggh7A6cHmM8utQ5atHYsfD002n+yKpL3cwiSZIkSSqpLl1g773TNorOnauibIB8MxxUYu+8k2Y37LknHHlk7jSSJEmSpC+oq6uaouEzWQuHGOPDwMM5M9SKfv3g449h1CgHRUqSJEmSii/LsZgqrb/8BW64Afr2hS22yJ1GkiRJklQLLByq3IIF0LMnfPObcO65udNIkiRJkmqFMxyq3DXXwN/+BhMmQPv2udNIkiRJkmqFKxyq2MyZcM45sO++cNhhudNIkiRJkmqJhUMVO/NMmDMHRo50UKQkSZIkqbQsHKrUI4/ALbfAGWfAt7+dO40kSZIkqdZYOFSh+fPToMgNNoD+/XOnkSRJkiTVIodGVqGRI+G55+C3v4WVV86dRpIkSZJUi1zhUGXeegsGDoQDDoBDDsmdRpIkSZJUqywcqszpp6ctFSNGOChSkiRJkpSPhUMV+dOfYPx46NcPNt44dxpJkiRJUi2zcKgS8+ZBr16w0UZw1lm500iSJEmSap1DI6vE8OEwbRr8/vfQrl3uNJIkSZKkWucKhyrw+utwwQVw8MFw4IG500iSJEmSZOFQFU47DRYuTKscJEmSJEkqBxYOFe6BB2DCBBgwIM1vkCRJkiSpHFg4VLC5c+Hkk2GTTdJxmJIkSZIklQuHRlawK66Al16C+++HlVbKnUaSJEmSpH9zhUOFeu01uOgiOPRQ2G+/3GkkSZIkSfoiC4cK1acPhABXXpk7iSRJkiRJ/8ktFRXovvvgrrvgkktggw1yp5EkSZIk6T+5wqHCzJkDp5wCm20GffvmTiNJkiRJ0pK5wqHCXHYZvPIKPPggtGmTO40kSZIkSUvmCocK8uqraRvFz34Ge++dO40kSZIkSUtn4VBBeveGli1h2LDcSSRJkiRJ+mpuqagQ99wDv/89DB0KnTrlTiNJkiRJ0ldzhUMF+PTTNChyyy3TKgdJkiRJksqdKxwqwJAhMH06TJoErVvnTiNJkiRJ0rK5wqHMvfwyXHopHHUU7L577jSSJEmSJH09Fg5lLEY4+WRo2xYuvzx3GkmSJEmSvj63VJSxu+6CP/wBrrwS1lsvdxpJkiRJkr4+VziUqU8+SQMit9kGevXKnUaSJEmSpOXjCocyNXgwvP463HYbtPL/JUmSJElShXGFQxl68cU0s+G442CXXXKnkSRJkiRp+Vk4lJnPBkWuvHI6nUKSJEmSpErkYv0yM3EiPPQQjBwJ66yTO40kSZIkSSvGFQ5l5KOP4NRTYbvtoHv33GkkSZIkSVpxrnAoIxdeCG++CRMmQMuWudNIkiRJkrTiXOFQJp5/Hq68Erp2hZ12yp1GkiRJkqTmsXAoAzFCr16w6qpwySW500iSJEmS1HxuqSgD48fDww/DmDFQV5c7jSRJkiRJzecKh8w+/BBOOw3q6+HEE3OnkSRJkiSpMFzhkNmgQfDPf8LvfuegSEmSJElS9XCFQ0ZTp8KIEdCtG+ywQ+40kiRJkiQVjoVDJjFCjx6w+uoweHDuNJIU6deeAAAH6klEQVQkSZIkFZZbKjK55Rb4859h3DhYa63caSRJkiRJKixXOGTwwQdwxhnwgx/A8cfnTiNJkiRJUuG5wiGDc8+Fd9+F+++HFlY+kiRJkqQq5NvdEnvmGRg9Grp3h+23z51GkiRJkqTisHAooUWLoGfPNLPhootyp5EkSZIkqXjcUlFCN94Ijz2WbldfPXcaSZIkSZKKxxUOJfLee3DWWfDDH8Ixx+ROI0mSJElScVk4lMiAAfD++2l+g4MiJUmSJEnVzre+JdDYCNdeC716wXe+kzuNJEmSJEnFZ+FQZAsXQo8esM46cP75udNIkiRJklQaDo0ssoYGeOopuPVWWG213GkkSZIkSSoNVzgU0bvvwtlnw49+BEcdlTuNJEmSJEmlY+FQRGefDbNmwdVXQwi500iSJEmSVDoWDkXy+OMwbhz06QNbbZU7jSRJkiRJpWXhUAQLF0LPntCxIwwcmDuNJEmSJEml59DIIvjb32DaNLjhBlh11dxpJEmSJEkqPQuHIth+e3j5ZVhvvdxJJEmSJEnKw8KhSDp2zJ1AkiRJkqR8nOEgSZIkSZIKzsJBkiRJkiQVXMkLhxDC+iGESSGEaSGE50MIvUudQZIkSZIkFVeOFQ4LgNNijFsAPwB6hhC2zJCjuJqa4Kmn0q0kSZIkSTWm5IVDjPHtGOPTi//+ETAN6FTqHEU1fjxsuCHss0+6HT8+dyJJkiRJkkoq6wyHEEJnYDvgiZw5CqqpCbp2hdmzYdasdNu1qysdJEmSJEk1JVvhEEJoD9wJ9IkxfriE+7uFEBpDCI1NlfRmffp0aNPmi9dat07XJUmSJEmqEVkKhxBCa1LZ8OsY4/8s6TExxrExxvoYY31dXV1pAzZH584wb94Xr82fn65LkiRJklQjcpxSEYAGYFqMcVipn7/o6uqgoQHatYMOHdJtQ0O6LkmSJElSjWiV4Tl/CBwDTA0hTFl8rX+M8b4MWYqjSxfYe++0jaJzZ8sGSZIkSVLNKXnhEGP8MxBK/bwlV1dn0SBJkiRJqllZT6mQJEmSJEnVycJBkiRJkiQVnIWDJEmSJEkqOAsHSZIkSZJUcBYOkiRJkiSp4CwcJEmSJElSwVk4SJIkSZKkgrNwkCRJkiRJBWfhIEmSJEmSCs7CQZIkSZIkFZyFgyRJkiRJKjgLB0mSJEmSVHAWDpIkSZIkqeAsHCRJkiRJUsFZOEiSJEmSpIKzcJAkSZIkSQVn4SBJkiRJkgouxBhzZ1imEEIT8FruHCprawPv5g4hFYGvbVUjX9eqVr62VY18XWtZNowx1i3pjoooHKRlCSE0xhjrc+eQCs3XtqqRr2tVK1/bqka+rtUcbqmQJEmSJEkFZ+EgSZIkSZIKzsJB1WJs7gBSkfjaVjXyda1q5Wtb1cjXtVaYMxwkSZIkSVLBucJBkiRJkiQVnIWDqk4I4fQQQgwhrJ07i9RcIYShIYQXQwjPhhB+G0JYPXcmqTlCCPuFEP4eQng5hNAvdx6puUII64cQJoUQpoUQng8h9M6dSSqkEELLEMIzIYTf586iymPhoKoSQlgf2AeYkTuLVCAPAlvHGL8DvAScnTmPtMJCCC2Bq4H9gS2BLiGELfOmkpptAXBajHEL4AdAT1/XqjK9gWm5Q6gyWTio2lwJnAk4nERVIcb4QIxxweIvHwe+mTOP1EzfB16OMb4aY5wH3A4ckjmT1CwxxrdjjE8v/vtHpDdmnfKmkgojhPBN4EBgXO4sqkwWDqoaIYSDgTdjjH/LnUUqkhOA+3OHkJqhE/D6575+A9+YqYqEEDoD2wFP5E0iFcxw0od5i3IHUWVqlTuAtDxCCA8B6y7hrgFAf2Df0iaSmu+rXtcxxrsXP2YAadnur0uZTSqwsIRrrkhTVQghtAfuBPrEGD/MnUdqrhDCQcA7McbJIYTdc+dRZbJwUEWJMe69pOshhG2AjYC/hRAgLTt/OoTw/RjjzBJGlJbb0l7XnwkhHAccBOwVPctYle0NYP3Pff1N4K1MWaSCCSG0JpUNv44x/k/uPFKB/BA4OIRwALAS0CGEcGuM8ejMuVRBgr+7qhqFEKYD9THGd3NnkZojhLAfMAz4UYyxKXceqTlCCK1Iw0/3At4EngKOijE+nzWY1AwhfdJxE/BejLFP7jxSMSxe4XB6jPGg3FlUWZzhIEnlbRSwKvBgCGFKCGFM7kDSilo8ALUX8EfSYL07LBtUBX4IHAPsufjn9JTFnwhLUs1zhYMkSZIkSSo4VzhIkiRJkqSCs3CQJEmSJEkFZ+EgSZIkSZIKzsJBkiRJkiQVnIWDJEmSJEkqOAsHSZJUdCGE9UMI/wghrLn46zUWf71h7mySJKk4LBwkSVLRxRhfB64Bhiy+NAQYG2N8LV8qSZJUTCHGmDuDJEmqASGE1sBk4HrgJGC7GOO8vKkkSVKxtModQJIk1YYY4/wQwhnAH4B9LRskSapubqmQJEmltD/wNrB17iCSJKm4LBwkSVJJhBC+C+wD/AA4NYSwXuZIkiSpiCwcJElS0YUQAmloZJ8Y4wxgKHB53lSSJKmYLBwkSVIpnATMiDE+uPjr0cDmIYQfZcwkSZKKyFMqJEmSJElSwbnCQZIkSZIkFZyFgyRJkiRJKjgLB0mSJEmSVHAWDpIkSZIkqeAsHCRJkiRJUsFZOEiSJEmSpIKzcJAkSZIkSQVn4SBJkiRJkgru/wMTV6wsXBCZ+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1296x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.array([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5], dtype=np.float64)\n",
    "Y = np.array([2, 7, 9, 12, 13, 14, 14, 13, 10, 8, 4], dtype=np.float64)\n",
    "A = np.array([[1, x, x ** 2] for x in X], dtype=np.float64)\n",
    "\n",
    "Q, R = np.linalg.qr(A)\n",
    "x = np.linalg.solve(R, Q.T @ Y)\n",
    "print(x)\n",
    "\n",
    "newY = list(map(lambda point: x[0] + (x[1] * point) + (x[2] * (point ** 2)), X))\n",
    "\n",
    "plt.scatter(X, Y, s=20, c=\"red\")\n",
    "plt.plot(X, newY, c=\"blue\")\n",
    "plt.rc('axes', axisbelow=False)\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "plt.legend([\"Approximated function\", \"Original points\"])\n",
    "plt.title(\"Approximated function\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jak widać rozwiązanie w dużym stopniu przybliża dokładniej kształt, który zakreślają pojedyncze punkty. Powyższą metodę można wykorzystać np. do dopasowania funkcji do wyników pomiarów. Jest bardzo skuteczna, gdyż aby rozwiązać układ wystarczy pomnożyć go przed odwrotnosć macierzy \\\\(R\\\\), a skoro jest to macierz ortogonalną to odwrotnością jest po prostu \\\\(R^T\\\\), co jest dużo szybsze niż obliczanie odwrotności w przypadku innych macierzy oraz redukuje błędy zaokrągleń, gdyż nie wykonujemy żadnych operacji arytmetycznych. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.3 64-bit",
   "language": "python",
   "name": "python38364bite4e77c707e6e47caa7d4a0acd81d8735"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}