05-Matrices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b596c56e",
   "metadata": {},
   "source": [
    "# 5 - Matrices\n",
    "- 5.1 Interpreting matrices\n",
    "- 5.2 Matrix terms and notation\n",
    "- 5.3 Matrix dimensionalities\n",
    "- 5.4 The transpose operation\n",
    "- 5.5 Matrix zoology\n",
    "- 5.6 Matrix addition and subtraction\n",
    "- 5.7 Scalar-matrix multiplication\n",
    "- 5.8 Shifting a matrix\n",
    "- 5.9 Diagonal and trace\n",
    "- 5.10 Exercises\n",
    "- 5.11 Answers\n",
    "- 5.12 Code challenges\n",
    "- 5.13 Code solutions\n",
    "\n",
    "Notes, code snippets, and the end of chapter exercises from the book _Linear Algebra: Theory, Intuition, Code_ by Mike X Cohen. \n",
    "\n",
    "Find more information about the book on [github](https://github.com/mikexcohen/LinAlgBook) and [amazon](https://www.amazon.com/Linear-Algebra-Theory-Intuition-Code/dp/9083136604)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6f9ec8a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92c4b406",
   "metadata": {},
   "source": [
    "## 5.1 Interpreting Matrices\n",
    "Applications of matrices\n",
    "- representing a linear transform or mapping\n",
    "- storing partial derivatives of a multivariable system\n",
    "- representing a system of equations\n",
    "- storing data eg features $\\times$ observations\n",
    "- representing regressors for statistical models\n",
    "- applying geometric transformations for computer graphics\n",
    "- applying kernels used in filtering or convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01981945",
   "metadata": {},
   "source": [
    "## 5.2 Matrix Terms and Notation\n",
    "Terms\n",
    "- row\n",
    "- column\n",
    "- element\n",
    "- block\n",
    "- diagonal\n",
    "- skew-diagonal\n",
    "- off-diagonal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5992f004",
   "metadata": {},
   "source": [
    "## 5.4 The Transpose Operation\n",
    "Transpose is defined for any rectangular matrix.\n",
    "$$\n",
    "A_{i,j} = A^T_{j,i}\n",
    "$$\n",
    "\n",
    "Proofs often make use of this fact.\n",
    "$$\n",
    "(A^T)^T = A\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f53582",
   "metadata": {},
   "source": [
    "## 5.5 Matrix zoology\n",
    "\n",
    "### Square Matrix\n",
    "A matrix with the number of rows equal to the number of columns.\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "d & e & f \\\\\n",
    "g & h & i \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "### Fat or Wide Matrix\n",
    "A matrix with number of colums $\\gt$ number of rows.\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & b & c & d \\\\\n",
    "e & f & g & h \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "### Tall or Skinny Matrix\n",
    "A matrix with number of rows $\\gt$ number of columns.\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & b \\\\\n",
    "c & d \\\\\n",
    "e & f \\\\\n",
    "g & h \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "### Symmetric Matrix\n",
    "A square matrix with elements $a_{i,j} = a_{j,i}$.\n",
    "$$\n",
    "A = A^T\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "0 & a & b \\\\\n",
    "a & 0 & c \\\\\n",
    "b & c & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Entries along the diagonal may not be zero.\n",
    "\n",
    "Examples\n",
    "* Covariance matrix\n",
    "    * $\\Sigma = E[(X - \\mu) (X - \\mu)^T]$\n",
    "    * Variance appears along the diagonal\n",
    "    * Covariance between non-diagonal elements is symmetric\n",
    "    * Imagine sample points from a bivariate Gaussian distribution in which $\\sigma$ in each dimension is different\n",
    "\n",
    "### Skew-Symmetric Matrix\n",
    "A square matrix with elements $a_{i,j} = -a_{j,i}$.\n",
    "$$\n",
    "A = -A^T\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "0 & a & b \\\\\n",
    "-a & 0 & c \\\\\n",
    "-b & -c & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Entries along the diagonal must be zero.\n",
    "\n",
    "### Identity Matrix\n",
    "A square matrix often written as $\\textbf{I}$ with 1 on the diagonal and 0 everywhere else.\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "0 & 0 & 1 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Column and row vectors are orthonormal $=>$ identity matrix is an orthogonal matrix.\n",
    "* Identity matrix is multiplicative identity $AI = A = IA$\n",
    "\n",
    "### Zero Matrix\n",
    "A matrix often written as $\\textbf{0}$ with elements $a_{i,j} = 0$.\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Zero matrix is additive identity $A + 0 = A = 0 + A$\n",
    "\n",
    "### ATA Matrix\n",
    "Matrix product $A^T A$ is one of the most important matrices in linear algebra.\n",
    "\n",
    "Properties\n",
    "- Produces a square matrix, even when $A$ is non-square.\n",
    "- Produces a symmetric matrix, even when $A$ is not symmetric.\n",
    "- If $A$ is full column-rank, then $A^T A$ is full-rank and the inverse $(A^T A)^{-1}$ exists.\n",
    "- Has the same row space as $A$.\n",
    "- Has orthogonal eigenvectors.\n",
    "- Is positive semi-definite.\n",
    "- Has non-negative, real-valued eigenvalues.\n",
    "- If $A$ has zero-mean, then produces covariance matrix.\n",
    "- Has distinctive appearance when rendered as an image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "204d4c90",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEwCAYAAAB4/k+CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABfIElEQVR4nO2debyVU//+r9VE0pw0U4aUJqmHnkIyZawMETImQyGZEhkfhEdkyhgRyixjopLoodlQigqVFDIUadD6/bF3v+85+7p27bvh1Mn1fr3Oq8773HvvdU9rr7PPtT4rxBhhjDHGGGNyp8jmboAxxhhjTGHDAyhjjDHGmIR4AGWMMcYYkxAPoIwxxhhjEuIBlDHGGGNMQjyAMsYYY4xJiAdQxhhjjDEJ8QDKGGPMWgkhXBhCGB1C+CaE8Hn6/z03d7uM2ZwEF9I0xhiTCyGERwE8EmP8eHO3xZjNjT+BMsYYkyv1AUzbkCcIIXwRQmi9UVqzFVAYjkdhaOPmwAOoQkZhuJBDCE+EEP6zuduxsSkMx96YDSGEUCyEMD2E8FWWTUrHGJes52MBADHGPWOMoxO2q1UI4aMQwm8hhMUhhA9DCM2TPMfmIP0nz4PXts36HI88z70ihFApw08OIcQQws6bu41bOx5AbSFsys4n/fzfhBCWhRCWhhB+SA9ytl/vBm8BFGQHArgTMf8IzgNQGUCdEEKpvD8IIdQEMHd9HrshhBDKAHgdwL0AKgCoDuAGAMs31mtsDkIIxTbC08wB0CnPczYEsN1GeN41z7cx2rjV4gHUlsMm6XwyODrGuD2AJgD2AnDVJnqdgmSTdiDp53QnYrZ6QghlAVwHoDuAvwE0yNhkTwBfrOdj8277/39hCSFcGUKYH0JYEkKYEUI4SDxkdwCIMT4bY/w7xrgsxvhOjPHT9HNcHkJ4MeM17gkh9M/zepeHED4NIfwRQngshLBjCOGt9Ou+G0Ion9G+nLYPIVQLIbwYQvgxhDAnhHBRnud5CkAtAK+lf3G9Iv3cV4YQPgXwR/oX53y/wIUQaoYQXko/588hhPuyHUsATwE4Lc/3pwN4MuNY9AohzEq3fVoIocOGtDGEsEv6U8CmeY7Bj//IT+djjP7azF8AygL4EamBwAoA+6xl228AHJz+/5UA5gNYAmAGgINyeVz6+9sBvJHn+14AZqWfaxqADhmPvQzApwB+AzAUwLZ5fr4XgEnpxw4FMATAf9I/qwdgNIBfkep8jxHtujz93H8AeAzAjgDeSj/fuwDKr2WfrgEwPo/7L4CrAUQAO69r/5DqgFYDWAZgKYAr8jz3lel2LQdQbM0xBLALgMUAmqa3rZY+f60397XkL3+t71f63vk4/f/PAZyT8fNLAZy5Po/N2HbNfVQXqU+0qqX9zgB2EduXAfAzgEEADs/sDwBUTfcd5dLfFwOwCMDeeV7vf+l+pXr6Z5PS/da2AEYCuC6jfevcHqkPICYCuBZACQB1AMwGcFjmvmZ8PwVATQAlM7cBUBTAVAB3ASiVfr1W6ziOM5DqZ4sCmAdgJ+Tp/wCckO6jigA4MX2sqq5vG9Pfn4NUP7odgOEA/ru5r9/N8eVPoLYM+gCYHWN8FsBMAI3W9YAQQl2kfttrHmMsDeAwpC7ydRJCqIFUR/R1Hj0LwH5IDeZuADA4hFA1z887AmgLoHa6fWekn6sEgFeQGohUAPA8gOPSPysO4DUA7yD16dqFAJ5Otz0vxwE4BKnfNI9GavDUG8AOSN30FyE7/wNQJoRQL4RQFMBJAAaL7eT+xRg7A/gO6U/nYoy353lMJwBHItUxr1ojY4yzkBpcDQ4hbAfgcQCDov+8ZwopIYQ6SN2fV6bVF8joh2KMd8YYH1+fx2bhbwDbAKgfQigeY/wmfW/lI8b4O4BWSA0KHgHwYwhhWAhhx/TPFwAYg9RAAUj1Uz/FGCfmeZp7Y4wLY4zzAXyA1GBvcozxLwAvIzU4QsLtmwPYIcZ4Y4xxRYxxdrp9J61jv++JMc6NMS4TP/sXUoOdy2OMf8QY/4oxjl3H8635FOoQANOR+qX6/xNjfD7G+H2McXWMcSiAr9Kvs75tRIzxEaTePz5GagB79Tqeb6vEA6jNzKbufDJ4JYSwBKnf+hYh9VsUgJxusnvSP1+M1KCoSdrvC6A4gLtjjCtjjC8AGJ/nZ9sD6JvuYEYilWXohPwk7dwyWWsHkuP+KbJ2Iu5AzFbG7QDezfNLwBcAGm/Kx8YYvwbQA8D1ABaFEIaEEKpl2XZ6jPGMGGMNpP48WA3A3Xk2GQTg1PT/T0WqT8jLwjz/Xya+z8yD5rL9TgCqhRB+XfOF1C9+O6p9yMPacmQ1AXyb9xe2HHgKwMlI/VL7ZOYPQwinhRCm5GljAwCVMrdL0MY1PJJ+rntjjIU6j7a+eAC1+dnknU8e2qc/rWoNYA/kuYlyuMl+yPP/P/F/HU41APNjjHkLin2b52dzY4yrM35WPaNdSTu3TNbagQCbrBP5x3cgpvATQmiF1KfArUJqgskPSP25ruGmfCwAxBifiTG2wv/92em2HB7zJYAnkD9n9QqARiGEBgCOAvB0Lq+/gcwFMCfGWC7PV+kY4xF5myset7bii3MB1EqSu4wxfotUFvQIAC/l/VkIYSek+qnuACrGGMsh9SfWsAFtREhNQLobqcjF9SGECrm2d2vCA6jNyObofNKPex+pDui/6XbkcpNlYwGA6iGEvNvWSv/7PYCaIYQiGT+jT4g2hLV1IEBO+5ets8jaibgDMVsD6fu2H4AHkcokNUl/HQqgXAih1qZ4bPrxdUMIbUII2wD4C6lfllaL7fYIIVyajh6smQ3YCak/3wMA0p9WvwDgGQCfxBi/W/febzCfAFiSDlyXDCEUDSE0CPnLKyxEKhuV5DkXAOgbQigVQtg2hNAyh8edDaBNjPGPDF8KqX7sRwAIIZwJDvgnbSMA9AcwIcbYBcAbSF0D/zg8gNpMFFTnsxbuBnBICKExcrvJsjEOwCoAF4UQiocQjsX//WnsY6Q+rboi/bPWSGWchiRoZ65k60CAde+fOxDzT+UUpP7kdHmM8Yc1X0gPDrD2OMGGPBZIRRD6AvgJqU+4K0PPDF4CYB8AH4cQ/kBq4PQ5Ur9s5mUQUr98Zv75bpMQY/wbqU+7miD1C9xPAB5FKme5hlsBXJP+5PuyHJ/zaAC7IpXNnIdU8Htdj5sVY5wg/DQAdyLVTy9E6vh8mLFZojaGENohlTM7P616AmgaQjhlXY/d6tiQBLq/1v8Lqb/Tfwtg+wwfAPwO4Kgsj/sGqZkXjfB/HdVipLJF1dbyet8gzyyKtBsA4MX0/29OP89PSA3s3gfQRT0WqT8bDs7zfTMAk/F/s/CG4v9m4e2Zfq7fkDG7L8tzDwZwfZ7vuyD1J86c9inti4Fn4a1t/9oh1Vn9CuCytRyvNce+HVKfolVI++2RykOdsrmvK3/565/6hdSn238CKLO52+Kvf8aX18IzxhhTqEnHBPohNXg6a3O3x/wzcIFAY4wxhZaQKjy8EKlP9Ntu5uaYfxD+BMoYY4wxJiEOkRtjjDHGJMQDKGOMMcaYhGxQBiqE0Bap6dxFATwaY+y7tu0rVaoUd9ppp3zuzz//pO2++eYbcpUrVyb3119/katWjWtJfvXVV7I95cuXJ1esGB+SVau4KGzx4sVz2q5MmTLk8pdMSvHbb7+RK126NDm1z7NnzyZXv359coA+tjvuyIVzlyxZQq5GjRrkZs6cSW7nnXcmp87B6tVcdUEd17p1M1d+AT7//HNy223HawirYwgAy5fnVveyVi2uJrF06VJyf/zB1ROKFi1KbuHCheQqVeJ6nt9//z25IkX49x11bL77TpfAUef5559/JrfDDjvk1J4qVaqs8/n+/PNPrFixIpd6YpuFJH1Y2bJlY+Y+L168mLZT1/XKlSvJbbPNNuTUNaMeCwClSvGa4+q6Vm7bbbclp/oltS8q9qH6JXW9qntU9Zvqsdnao15bUbZsWXLq2KjjrfqlZct4lZMSJUqQK1myZE6PVe892c69unbUeVmxYkVOr6O2U/usttt+e65zrLZT17a6DtU9Beh9znX/1DWmzotq4/fff/9TjJE7RWxABiq97thMpJbPmIfU8h2dYqruhGTvvfeOH3/8cT43ceJE2u7ss88m1717d3LqzbtPnz7kDj/8cNmek07iJYvKlStH7tdffyWn3vTUdoceeig5dWG+/vrr5A488EByX375JblOnTJXRgEmT55MDgDOOOMMcj169CD3wQcfkLvtNq7Tecghh5B77LHHyB1xxBHkVCeiBsCjR48mpwYOjRpx2ZnWrVuTA/RAUnH//feTU8cm87oG9LXUr18/cueccw656667jpzqnFVbLrjgAnIAcOmlmWVzgEGDBpE777zzyN10003kevbsSW7w4PzLEH7wwQf49ddft8gBVNI+rG7dunHAgAH53JAhXNJMDbB/+OEHcrvtths59QvX/Pm67myLFi3IzZrFqzmpX7D22GMPcmrQogYn6k1L9cVqgKfu0Z9++olctl981C/cM2bMkNtmcuSRR+b0WPVLjvplYfr06eSqVq1KrnFjXlhC/QJYoQLX4lW/uAD62lHnRf0yVbFiRXLz5s0jp/pi9XzqOlTtVte2ug4z+5A17LLLLuTUvaH278cffyTXsCHXq1aD7D59+kyMMTZTbdqQP+H9C8DXMcbZMcYVSBVHbLcBz2eMMQWJ+zBjzHqzIQOo6si/Vtg88BpnCCF0DSFMCCFMUL9pGGPMZmKdfVje/kt9wmyM+eeyyUPkMcaHY4zNYozN1J+9jDFmSyVv/6X+JGuM+eeyISHy+QBq5vm+BtaxSOysWbNw7LHH5nMqxDpmzBhyKsCswmvqb5hNmzaV7cn2d/ZMMtsMAO+//z459TdjlUdReaBp0zh2of7mq8LKw4cPJ6dyGIAOCKuA95NPPknu+eefJ6cyPSoP0bYt17dTfzt/+mleRP2EE07IabszzzyT3D333EMO0NktlSXLnPQAAA8//DA5lQ1QwUV1XNVkBnXNqmySys5l+6S3W7du5B58kJfwU/fFRx99RE5dd5l5m1zD+puJRH3Y4sWLKfN00EEH0XZ33XUXORVYnTt3Lrk999wzp+0AnSdS954K6qrMi+ozVKbn999/J6euOTUxRk1kUNeI6scBHSJXnwyqvKDKO6ngvMrVqCzYL7/8Qk7dy2pCjjrW6jhkC9Or11aBc/XaKiulXkc9VmX5VJ5UnRN1jaj3cHXdAMDff/9NLtcJOCoDpULptWvXlq+djQ35BGo8gN1CCLVDCCUAnARg2AY8nzHGFCTuw4wx6816fwIVY1wVQugOYDhSU4AHxhi/2GgtM8aYTYj7MGPMhrBBdaBijG8CeHMjtcUYYwoU92HGmPXFlciNMcYYYxKyQZ9AJWXlypVYsGBBPqeKM6oA2kUXXUTuk08+IffZZ5+RU8XlAB3SPuuss8ipgJ0qVHn33XeTU6FhVWhShZpVwLd58+bkVJE3FZADgEmTJpFTQenOnTuTU8Htk08+mZwKH/bv35/chAkTyKkQrArOq+KaqiibqooO6IrIl19+OTkVfFTh3RdffJGcCg7fcccd5FQoXV3vt956KzlViFa9BgA8++yz5FTRujZt2pBTEwhUezKrAKuQbmFl9erVFP5VgXE16UQVLFUTY6ZMmUJOTUYAgC++4L82qnt51KhR5FSfoaruqwkmqviyKoqrmDNnDjkVcs9WgVtN1lDvF+o5v/32W3Kqn1QBdBVeV/uc6/6p11VhbFV8EtDHQQWyVf+lgvzq2KgJBKrdKmCvjpcKcqtzpx4L6GtRvTerSubVq1OFJfn+n2uB5f//+om2NsYYY4wxHkAZY4wxxiTFAyhjjDHGmIR4AGWMMcYYk5ACDZHvuOOOuPjii/O53r1703a77747ORVAU2FNFSTu0KGDbM8BBxxA7sMPPyR30kknkVPVap955hlyKnD57rvvkhs7diy5Vq1akbv++uvJqZDofffdRw4AbrzxRnKqeu7AgQPJqVXKVSXyQw89lNwll1xCTh0vdaxV+FmtZn7aaaeRO+yww8gBOmh49tlnk1OVpp944glyV111FTlVmX7RokXkVOXd22+/nZxa/V2F5FUlfwD49NNPyan7SlVJVo9VEx/eeOONfN+rqvuFlZUrV1IlZjVRQAXGzz//fHKZxwrQVaJVHwIARYsWJafui5o1a5JTFZxV4DjXMLYKEqvK0SqYrELyKugM6HD5hlTlLl68ODkVJFYBbdVudbzUsVb3kzqfagINoCcLqPOiXlut/qH2T52DESNGkFPVydU1p1bByHZtK9TxUdXb1QQhda6SBsYV/gTKGGOMMSYhHkAZY4wxxiTEAyhjjDHGmIR4AGWMMcYYk5ACDZGXKVOGQqX33nsvbVenTh1yr7/+OrmRI0eSGzduHDlVYRoATj/9dHIqzK2CzSqA/txzz5FTFbNVGx944AFy33//PbkPPviAnAqoqmMD6P1TFV67dOlCbvbs2eReffVVcirUecIJJ5BTgVlV7VxVvx0wYAC5UqVKkevVqxc5AHjzTV7+TE1AUOF3FdBXkyFUReN27dqRU0HuypUrk2vYsCE5dU7U6wLAa6+9Rk6d52HDhpHLtcJ1ZnVfFSQurGyzzTbYbbfd8rm5c+fSdqrCuAqMd+rUiVyuFe0Bfb2qys4qQKsqxKtqzSpwrILJpUuXJpdriPz3338np/YjSRvVc9auXZucOjZqYoaaEKIqkavwuqp4rQL7KoifrSK7mmiVucoHoIPXKhz+559/klMBb3X8latQoQI51d+r6zBbsLxJkybk1AQcdU+qCVDqGCbtr/wJlDHGGGNMQjyAMsYYY4xJiAdQxhhjjDEJ8QDKGGOMMSYhBRoiX758OWbNmpXPHX300bRdx44dyX3yySfkevbsSU6FsV9++WXZnkcffVS2MRNVzbp79+7kVEBYVXi9+eabyQ0ePJjcpZdeSk6FyFWF1yFDhpADgGbNmpFTlauPO+44cqrCcvv27cl99NFH5FRF3fPOO4+cOidXX301OVUtW1W/3WeffcgBOlCqzosKPqp2q6rc6rGqSrsKvKpg+dSpU8lNmzaN3Kmnnkou27ZqX1TINDM8DeiwbWZIWFV6LqwULVoUZcqUyef23HNP2m7KlCnk1LWgAuMXXnghudtuu022RwVwy5YtS05dS+XKlcupjeocq9BwZr8O6FC0ap8KDavjCujQsNo/db2qe14FjlV7VJhbvVeoCS8qMK5Qk2CyVctWr924cWNyaoKQOn/q+YoV4+GBcs2bNyenJkCpKu3KZZs0ofpJdT2oCQ0q8L969Wpy2SrgZ8OfQBljjDHGJMQDKGOMMcaYhHgAZYwxxhiTEA+gjDHGGGMSUqAh8jlz5qBz5875XOvWrWm70047jdy7775LrmvXruRUhelsIXIVYL/oootycqqi7r777kvuyy+/JKdCbv369SOnjs1DDz1ELjPYCgB33303OQDYf//9yakK6g8++CC566+/ntzFF19MbtKkSeRURVwVcLzsssvIqfDgnXfeSU5VO58xYwY5QJ8XVRX/2WefJVevXj1y3333Hbnx48eTU1X2y5cvT05VyV+8eDE5FTreb7/9yAE61KtCr++99x45VUn7xBNPJNe3b99836twamFl5cqVmD9/fj6nQsiqurUKJquwrAqMZzuf6npQTr22CieratuZ+wvolQZUiHy77bbL6XVVFWx1DAFdoVxNolGVp9W1rqr2qzaqlQHUxBgVTFbbqf1QYfiKFSuSA/T+qfMyb948cqovLlKEP0tRz6cm6qhzr8LduVamzzV0D+hrR1WmV9eICsSr47o2/AmUMcYYY0xCPIAyxhhjjEmIB1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJKRAZ+GVLFkS9evXz+eOOuoo2k7NSLvqqqvI9enTh9yCBQvIqZlUAPD222+T22OPPcipZQ/UEiZq+ROV9Fdl7h977DFyajkDtX9q1kL//v3JAcCxxx5LTi2To86BKnP/2muvkVPLKDRp0oTcJZdcQk7NEhw+fDi5k08+mdznn39Obvbs2eQAoFKlSuTeeecdcl26dCGnzrNakkDN6lPt6dChAzk141Et7TNs2DByAwcOJAcATz/9NDl1Pbzyyivk1PIPY8eOJTd9+vR836slFAorpUqVQosWLfK5Ro0a0XZffPEFOTXzSfU1ankWNbMO4BmPgL5GKlSokNNrq5mpaikXNUttp512IqeW8VF9iNpn1WZAz2BWs13VMfv3v/9NTs32U7MH1Qw+tX+q72vZsiU5NStMzdbLNotVLduiZtypWW5q1q9aYke9d6kZd02bNiWn7nu1f+o6zOxD1pB57wG5zzodN24cOdVutfTaww8/LNsD+BMoY4wxxpjEeABljDHGGJMQD6CMMcYYYxLiAZQxxhhjTEIKNES+ZMkSjBo1Kp9TwTK1FMi3335L7qmnniKnAmTdunWT7VEhzCuvvJLckCFDyH366afkDj30UHIqsHfDDTeQu/XWW8mpELIKeqow3IsvvkgO0CHA119/nZwKoTdv3pzcSy+9RE6F0tXSBT179iTXq1cvcmo5A7W8hQpJ77333uQAYMSIEeTUPqulc2bOnElu8ODB5Hr06EFOhUJvvPFGcm3btiWnAqFHHnkkuWyhx+eff57cokWLyKmgu9pntYzCmWeeme/7+++/X7alMLJ8+XK6f9RxUYFqNcFB3RMqzJstRK4C48cccww5dQ7UOVbhabVEkVqmQy23pK4PNQlm5cqV5Nq1a0cOAD744ANy6tg2bNiQ3A8//EBOnYM5c+aQU8vcqGU/1POpCSuqH1aBdjXhCNDHTJ0r1R51bNR7bpUqVcipa1H1S+pYq6Vc1LFRywcBwPvvv09OXbNqUpWaxKH6YtUfrg1/AmWMMcYYkxAPoIwxxhhjEuIBlDHGGGNMQjyAMsYYY4xJSIGGyOvUqUOBxnr16tF2Z5xxBjlVmVlVfVVh3jfffFO2R4XQVXj3gQceyKmNKhB/6qmnkqtcuTK5e+65h5wKr6vjpQLRAwYMIAfoSsAXXnghORW8P+GEE8jVqFGDnAozqqC0qm593333kfv555/JqTDpvvvuS06FwAFdYVyFHFWo/a+//iKnwsSHHXYYudtvv53c1KlTyanK7SoQr6r2jh49mhygg6ePP/44uQkTJpA77bTTyKl7LXPywi+//CLbUhhZvnw5ha9VBfrMiTIAULNmTXIqNPzjjz+SUxNjAF2tWwXG1eSWt956i5wKHJcoUYKcutYV6r5VE0JUGFtNGgL0vaecugcaNGhATlU2VyFydRzU66pjuGLFCnIqYK/C0+p1AR2AVtei6tNUJfkiRfizFLUvyqnK+2qigarmro51tkkTKqyuKoerKugqRK4mgpUqVUq+djb8CZQxxhhjTEI8gDLGGGOMSYgHUMYYY4wxCVnnACqEMDCEsCiE8HkeVyGEMCKE8FX6X64WaYwxWwDuw4wxm4JcQuRPALgPwJN5XC8A78UY+4YQeqW/5xLeGfz2229U9VqFmlXwVFXgPvvss8l1796d3GuvvSbbc9ZZZ5Fr1KgRORWoVgHvgw8+mJyqyq3C5qpy62233UbuiCOOINexY0dyKlAI6ND3HnvsQe7rr78mp6oXq+rdqqp6hw4dyKlKt1dddRU5FexXYWxVKb1OnTrkAF21vFy5cuRatGhBbvfddyfXuHFjcioc26xZM3L9+vUj17VrV3Lq2vzqq6/I1a9fnxwATJs2jZyqdKyOmQqrH3jggeRuvvnmfN9nC0AXME9gI/Rh2267Ld0r8+bNo+3Uda0qLocQyKlrUIV+AX3fqntUBcYvueQScqq/qVu3Ljm1uoK61osV47cXdRxUGF5VmM62rXI77LADOTUxQ1UxVxN1VOhbBbl33XVXciqoro6N2o9s948KT6s2qhC6msSkzqm6Fl999dWcnk/dA6tXryanKu+rcwfoCv9qsoG6LyZPnpzTduo6Xhvr/AQqxjgGQGYsvh2AQen/DwLQPtGrGmNMAeE+zBizKVjfDNSOMcY1ixr9AIA/RkoTQugaQpgQQpigphcaY8xmIKc+zP2XMSYbGxwij6nPGLN+Th9jfDjG2CzG2KxkyZIb+nLGGLNRWVsf5v7LGJON9R1ALQwhVAWA9L/JljA2xpjNi/swY8wGsb6VyIcBOB1A3/S/nCwTVKhQASeddFI+N3DgQNpu1qxZ5F544QVyjz32GDkVVFPBckBXq61duzY5VY1aVW5VldGnTJlC7o477iA3dOhQciqo3q5dO3IqcKkC9gCwzz77kHviiSfIqSqt6jfwGTNmkFNB8AsuuICcCofvsssu5FQleBUenDhxIrl//etf5ABdJVwdx759+5JTFezPPfdcciqEqcK7nTt3Jrd06VJyN910E7mxY8eSUxXZAR2OVedUVRg/8cQTyanwe2aFeBWW3UJI3IeFEOi+V+epVq1a5H766SdyqhK/qlpdrVo12R41sUZdw6p6tAqMX3PNNeTuvvvunF5D9dk777xzTo9VFcuzfdqnqlmr+0xV4FZ9tgosq5UUVPhdPZ8KY6tgv9oPtVKAmmgAAEuWLMnp8bkemzFjxpBT152qqK/aokLkquq4CtiriuwAcNBBB5FT+6cmoamgu6okn23yQjZyKWPwLIBxAOqGEOaFEM5GqtM5JITwFYCD098bY8wWh/swY8ymYJ2/HsYYO2X5EQ8HjTFmC8N9mDFmU+BK5MYYY4wxCfEAyhhjjDEmIQWa8Pzmm29wzjnn5HNFixal7VT120cffZTcZZddRu6zzz4jp6otAzogqQLn7du3J6fCdKratgpUq3Be06ZNyalQtKqmqx6rws+ADqmqSsUq/KtC3yqMrYLzqoKtCp4ecsgh5NTxUlW5r7/+enJXXHEFOUBX7VXHrEuXLuTUOVXV7tW1pCYQqGrWquK8qpqsKvmfcMIJ5AB9HFVFd1WJXFU8P/nkk8mdcsop+b5XgdXCyurVqyl4WrVqVdpOHSu1MoAKT6vgdbZQrdpWBaBVNWrVx6r+UK3W8Pbbb5NT15aq1K1WYVi1ahW53377jRyg3y9Un6+qTKuq3mrSkTpeaju1eoS6HlRQWt3zKjivAtqADjurYLoKShcvXpycWnFBHa/FizPr0QJ77rknufHjx5NT14iq3K6eD9D3gdo/dY2p8zd37lxyFStWlK+dDX8CZYwxxhiTEA+gjDHGGGMS4gGUMcYYY0xCPIAyxhhjjElIgYbIq1evTtWUmzdvTts1aNCAnKpCqqq+qqrh1157rWzPtGnTyFWoUIGcqsKsAnG33HILuW7dupEbOXIkuYYNG5JToXQV5L7rrrvIqQrhgA7YqaCoOgd9+vQhpyp9d+rEZXeOP/54cqq6/KmnnkpOVUVXVdpVgFNVWQd0qF1VCVcVbFX1YhWIV4HeyZMnk1OB9quvvprciy++SE6FNVXAHtAh9J49e5JTYVYV7FSTHDIndqjAamElxkiTMH7//Xe5XSaqcr6qbq2qk6vKykBqUk4mKnw9c+ZMcqrvVMFkFRhv06YNuenTp5NTwV1VWVsFqtWKEICuqq7C9Cp8rcLq6vypCuOLFvFKPypErkLuP/74I7lSpUqRU32ImiQCZJ8YlYm6dtR7jTquuT7fpEmTyG2//fbk1LFW17C6VwB9b6iVANR5zrUKetL+yp9AGWOMMcYkxAMoY4wxxpiEeABljDHGGJMQD6CMMcYYYxJSoCHyFStWUPVPFQwrV64cuRtvvDGn11CVsVWQG9DhvpdeeonciBEjyB188MHk+vXrR27vvfcmp0LpKrjYsmVLcqq69bBhw8ipyrKArmSugsS77bYbORUsV05VE1eBfRX+rFWrFjkViFf717VrV3Kq0i2gK9urSszqWly2bBm58847j9ytt95KToV8P/nkE3Jff/01uQULFpBTx0Zd1wDw008/kZs9eza5/v37k1MBdlVdO/M+VW0urPz1118UyFbHVF0fqmKyCrGqILFyALDTTjuRU0FdxcKFC3N6HTVZRgXGd999d3JqoocKY6tJC+q4AjrErO4VFYhXgWXF1KlTyal+QPXZ7733Hjn1HqCqtKsgt6r8Deh+UgWqV65cSU7d86oqt3oPUJNg1P6piQHq2lR9SLYQuQq/q+OoAvbq/lPHRoXN14Y/gTLGGGOMSYgHUMYYY4wxCfEAyhhjjDEmIR5AGWOMMcYkpEBD5MWKFaPKzioIqSomq+rPt99+O7mnnnqK3F577SXbU79+fXIqMP7www+TU9W/n3/+eXKqkrWq7qtCuvfffz85FcJU22Wrvq78G2+8QU5VDB49ejQ5FdZU50VViFcBaFXxWlUnV8chs0o0AJx++unkAGDChAk5bauqLqswowoEf/bZZ+QefPBBciqges0115B78803yakqwG+99RY5QF/vKoyvzp+avKAmQ2RWg9+aKpEXKVKEKkirSQGKv//+Oyenwq7qHgP0sVVV7X/++WdyxYpx168mdaiQrqowrgLjJ598Mrl33nmHnNqPbJW21UoKJUqUIKfOiwqmq4rgqpK1Cq+rEHLdunXJqeOlAu3qnKhJVtlQldFVVXy1f+q6Uf2Xeh9W50/tc8WKFcmpY6jaB+iAtwreq/dcdWzUta2Ow9rwJ1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJMQDKGOMMcaYhBRoiHzu3Lm46KKL8jlVNfn8888nN2XKFHIq8PXaa6+RU6FYABgzZgy5L7/8ktz+++9P7rLLLiN37LHHknv//ffJffDBB+RU2O+ss84iN3ToUHIDBw4kd/bZZ5MDdDBQVSC++eabyamApAoh77jjjuTOPPNMcioQ+vjjj5P77rvvyKnK3yqgeuCBB5IDgNNOO42cOi+jRo0iV6lSJXLquKoJDepaUkFKVaFXtUVV3h8+fDg5QFdJVufgggsuIPfII4+QU1WqM8P56vkLK9tttx0aNWqUz6lrXYVdM8P12VzZsmXJZaugrV5bVX5X4V11Dauwujp/aoKJmtSh7ke1ooTarmjRouQAXc1ahcNVZXQ1+UMdW7V6gQpPq/eumjVr5vRYdU5U8F0FnQE92UadU9XG8uXLk1MTcNR26v1RbafaoiqlV6lShZyqbA4Ae+yxBzk1QUKhzr1qt7puVL+5Bn8CZYwxxhiTEA+gjDHGGGMS4gGUMcYYY0xCPIAyxhhjjElIgYbIQwhUtfTjjz+m7VQA7dJLL5XPl4mq6H3bbbfJ9txzzz3k5s2bR06FAJs2bUpu+vTp5F544QVyqrK2CgPvt99+ObXviSeeIPfVV1+RA4BTTjmFXLNmzch17NiRnAo7X3311eTOOecccmqfVeXj6667jpyaQKDC9CoEqyqJA7oi+9SpU8k1aNCA3J133klOXccDBgwgp46DqpQ/f/58ct26dSOnQphNmjQhB+ggeO/evckNGzaMnAr//uc//yGXWTlZVSsvrKxatYoCyyrkqyYAqOrKKhCtJpNkq46s+j91jy5cuJBchQoVyKlArmq36oNKlixJTlWoVoHxgw8+mNy4cePIATrUro632u7rr78mp95rZs+endN26tioCtzLli0jp4LvKtSc7dznWl1evY7qM1Qb1Wurit4qeK2q7Kv+WVV4zzaBYMmSJeRU5f5cq/6r/VOVzdeGP4EyxhhjjEmIB1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJKRAQ+QVK1akCtAq0KYC440bNyZ37733klNVlFXIENCVTceOHUtu/Pjx5I477jhyo0ePJqeqk99///3kTj/9dHIffvghOVWBW1WlHTRoEDkAuPjii8mpML2qoq0ChGqfZ86cSU6FiVUAXVUI79evH7lrrrmGnKqUvnjxYnKArhyuguUq4NqpUydyKrRfr149cqrauWq3On9du3Yl16VLF3IqNAwAffr0Iacq1qvrSa0O8Mknn5B79dVX832vKiEXVooUKYLSpUvnc6pyuApyq6rhKrCqjr0KKwM6CN6uXTty3377LTkV6FVBcFUdu3bt2uRUWPn7778npwLCKjCuJjcAwLPPPktOBcbLlClDTk38UcFkdV7UeVYhfnUMDz30UHJqwpGaOFKuXDlyAFCqVClyal/UShHNmzcnp97jVFV1dZ4rV65MrkgR/mxGXXMqgJ5tApTqJxVz584lpyaC1ahRg5w63o899ljW1/InUMYYY4wxCfEAyhhjjDEmIR5AGWOMMcYkxAMoY4wxxpiEeABljDHGGJOQAp2F99dff2HGjBn5nJqhoGadjBgxgtzgwYPJtW/fntz1118v26NmcqmZSmpJlQcffJDclVdeSU7NBlEzmk4++WRyZ5xxBjk1C0XNhvriiy/IAXpZFLXkglpS5bXXXiN37rnnkrvvvvvIqZl+H330ETl1XCtWrEjuoYceIqdmMKnXBfSMnJdffpmcOg5qBswbb7xBbsWKFeTUzM/MmWsA8OWXX5JbunQpuZtuuolctllbrVu3JqdmcqnZfn379iWn7qvMGT5jxoyRbSmMrF69mpZ/UEt3fPfdd+TUzLzMZW8APest2/ISmTMCAT3LUy13oc77NttsQ07NmlPtUUtyqNlxavkatZ3q5wBg3333Jaf6pRIlSpBT50odG7VMikK1W72GupfVTFnVZrX8DKCPg1piR80UVMvpqPOn3rtUn6auBzXjTr3PqHtFHVdAL72SbYZ9Jmpmn2pPtlmP2fAnUMYYY4wxCfEAyhhjjDEmIR5AGWOMMcYkZJ0DqBBCzRDCqBDCtBDCFyGEi9O+QghhRAjhq/S/5Td9c40xJnfcfxljNhW5hMhXAbg0xjgphFAawMQQwggAZwB4L8bYN4TQC0AvAJyizkPRokUppKUCxyoY+80335BTy3G0adOGXLFiejdfeOEFcvvssw+5t99+m9yTTz5J7plnnsmpPcOGDSOnAu0qrFytWjVyxxxzDDkVpgd04E8FGs866yxyammS4cOHkzvllFPIqWUd1BITl112GTkV4lfL+KjwZ7ZAqArHPv300+SWL19OTk0W6N69O7kOHTqQU9ecmgzRsGFDcrvtthu59957j5xacgcAXnzxRXKvv/46uWnTppFTS7Jcfvnl5B555JF836ugZgGz0fovNQlG3U+5hqzVtamOs1qGAtDLb3z++efkcg2RqwCt6mN33nlncirsrELRaikQFQTOFiRWgfFmzZqRU32xWuJDTVBR13+lSpXIqXC4miSijmuuS7lkO/cqZJ9rwFv1LZnXNaBD6eoYLliwgJxa7ke9RsuWLcmpJcwA/X6oJmeoe0iF1dXzjRo1Sr52Ntb5CVSMcUGMcVL6/0sATAdQHUA7AGsW7BoEoH2iVzbGmE2M+y9jzKYiUQYqhLAzgL0AfAxgxxjjmqHnDwB4tczUY7qGECaEECao34SMMaYg2ND+S/22a4z555LzACqEsD2AFwH0iDHmK+YRY4wAonpcjPHhGGOzGGOzXGs2GGPMxmRj9F9bwJ8jjTFbEDkNoEIIxZHqfJ6OMb6U1gtDCFXTP68KYNGmaaIxxqw/7r+MMZuCdYbIQ6qU6WMApscY++X50TAApwPom/6XyynnwEsvvUROhfNq1qxJbvbs2eS6dOlCToWVAeDSSy8lp0K5qjr5c889R059wqZCjyoofcQRR5BTVXvnzJlDrn79+uRU1WlAV2Q99thjyanQpPoTbJ06dch17tyZnApeqxCfCh+qfXn//ffJnX766eTUdQPoiQpqYoAK5aqwZr9+/cipSvIq1Nm/f39y6tyr4LwKhB5++OHkAOCSSy4hpyYBTJw4kZy6L9QqApltXLRo845LNnX/pYLSqh9Qf/5T51hNFFCTPABg8eLF5NT1MHr0aHIqnPz333+TU6HoJUuWkFOVp1VV9d13352cCoyrKtiAPt4qMK7uPXUc1XuImnSkzqlazUBVIl+2bBk5NYlll112IZdtRQlViVzdtzVq1CCnjreaGKD2RV0jagUI1c+pMH2RIvwZjgr2A/qYKafOVdWqVcmpfVHnVK2YsYZcZuG1BNAZwGchhClp1xupjue5EMLZAL4F0DGH5zLGmILE/ZcxZpOwzgFUjHEsAF5QJ8VBG7c5xhiz8XD/ZYzZVLgSuTHGGGNMQjyAMsYYY4xJSC4ZqI1GmTJlcMghh+RzDzzwAG336aefkttxRy7TogLoL7/8MjkVkAN0NddatWrl9NozZ84kp0LIqkrriSeeSE6FdLt160ZOVUq9+OKLyT344IPkAB24POyww8hdf/315FR1chUWVKF0FRY84YQTyDVp0oScqqCtws8dO3KM5ZxzziEH6BCtqhKurk8Vjr3uuuvIderUidysWbPIqWr8yt1yyy3k1IQLFewHdGC2cePG5NS12Lx5c3JVqlQhl3ndqarAhZWyZcviyCOPzOdUdWU1aUUFr9WECeW22WYb2Z5///vf5FQV+gYNGpBT95m6R1MVHvKzatUqcmqlCFV1XIXcVRVzFeYFdLBZVcdWgfGrrrqKnFp9QB0HFYBWk2rU+0fZsmVzej41qaZ06dLkAL0qhHq8unbU+VNO9Q1qEoCaxKT6OTUBQF2bqmo4oFekUJOi1MQH1ffttdde5NT9p/riNfgTKGOMMcaYhHgAZYwxxhiTEA+gjDHGGGMS4gGUMcYYY0xCCjREXqxYMapwrUKFKhSogskqRK6qLd91112yPSrsWb16dXI9evQgt88++5AbOnQouSlTppDr27cvOVU9WgW0zzvvPHIqOJ8tiKcYOXIkOVV9VQX2VBVaVbldBbTVsXn00UfJXXPNNeRUoFpNClBhS0BXu/3ggw/IqbB6yZIlyXXv3p3c5MmTyd16663kVDCzUaNG5FQQedy4ceRUWBbQEwPeffddcueffz65Dh06kFOTKzIrv6truLCyfPlyCo2nCp3nRwV31Tp6Krir+j4Vvs32eHW8VRBZVdhXEytUaFtVUFeoULMKm6tVD9S1DujgtrqXVWhYBcbbtm1LbtiwYeTUPb906VJyatJEhQoVyKmgunqsqqAN6IC+qsCtjpeaBKOuWXUtqvPy/fffk1PXjTpequq7qmqfrT1qn9VECnX+fvnlF3JqX9aGP4EyxhhjjEmIB1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJKRAQ+TffPMNzjjjjHxOhb5ViHXIkCHkKleuTE5Vgj3++ONle0aPHk1OhftU2ExVLL3xxhvJnXrqqeRU8FeFPy+//PKcXleFyG+//XZy2TjmmGPIqerkc+bMIfff//6X3EknnUSudevW5Nq3b09OBS4zqz8DwLRp08ipSrVvvfUWOQAYM2YMuRo1apBTIU712GbNmpEbMGAAORVkHTRoEDlViVkFf9X+qSrAgA7Ub7/99uQeeughciqYrsK7xx13XL7vVZi3sLJy5UqqcD1//nzaToW+VXharZCgKn+rxwI6bKvuUeXUtaAmnqh9UX3k1KlTyamQ+6677kpOXSPZJh/8/PPP5FRfoMLJqt0qMK4mCL3//vvk1H22yy675NQWtc9///03uSTkGsguX748OTVBSFVzVytwqErkqi0qvK6C5RMnTiQH6OtJBe/VPqvXUcFytXLI2vAnUMYYY4wxCfEAyhhjjDEmIR5AGWOMMcYkxAMoY4wxxpiEFGiIvHr16hS07tWrF2133333kfvjjz/I7b///uRUsE8FuQFd/VsFwY844ghyquq4CmOXLVuW3L333ktuwoQJ5FT4XVXvVoFvFbwGdFX28ePHkxs8eDC5+++/n5wKRT/22GPkVJXwm2++mZyqOK8CzC1btiR31FFHkVNBSAC44ooryKlq0So0qar7qury++67L7lXX32VnApHPvHEE+TUfdGtW7ecHgvoYHq1atXILVu2jNy8efPIqWPbqlWrfN9vTZXIixcvjipVquRz6lpQVZhVWF8d+5UrV5JTk2UAXT1aPacKgv/555/kVEVwtS+LFi0iV65cOXJq8o2qdq5eV/Wb2VCPVxWqVRtVkFgFxlV/8cADD5BTEysyrxlA3ztqAkE26tSpQ04FpVXoXu2zOt5qX1RAW7VFVfRW50kF0LMFudVxVO1Wk7SaN29OTgXd1SSHteFPoIwxxhhjEuIBlDHGGGNMQjyAMsYYY4xJiAdQxhhjjDEJKdAQeYwRK1asyOdUoFpVXFYVXi+66CJyHTp0IPf888/L9tx0003kqlevTq5Hjx7kVDivaNGi5Dp16kROVSAeMWIEOVXR+7zzzsvp+TIrvq9BVcfea6+9yKlAddeuXcn98MMP5B588EFyKox98cUXk6tZsya5HXfckZwKNd9xxx3kDj74YHIAMHLkSHIqFKpe54ILLiCngvN33nknuY8++oicClKqc/ree++Re+SRR8ipCvvZnlNVRFbV7lXlarWKQGZF4zZt2si2FEaWLVuG6dOn53MqKK1C+GqigNpu+fLl5NQEGkDfoyqcrILu6nVUpW51fyunHqsC8apytDo2IQRygD4WX375JTk1YUIdBxW8VhORVN/QuHFjciqAroLO6r1CtU9Vhwd0wFtVZFeTOFTAW7FkyRJy6vwp1LlX4XB1/NV+APqaUEFwFZJX1fx/+ukncg0aNJCvnQ1/AmWMMcYYkxAPoIwxxhhjEuIBlDHGGGNMQjyAMsYYY4xJSIGGyMuVK0ch74kTJ9J2qvK3CkyqCtwqdKcCwwDQu3dvcirwd+KJJ5JT4eRhw4aRU0G81q1bk1MVqr/99ltyKvBdpkwZco8++ii5bI9X+7zffvuRU9W/3377bXJdunQhp87fWWedldNrqNC2CuCee+655FSFcECHw/fcc09yqpq4Co+qiQpq4sNuu+1G7qWXXiKnwvQffvghOVUx/tprryUHAL///js5VbFehTBVJXkVks/cP3X9F1ZKlCiBqlWr5nPly5en7ebMmUNOVcZWQVtVXVxVAwf0taRCw8qpistqXzL3F9ABaDXBoW7duuTUda32T12D2bbdY489ctquVq1a5FQoWt23KrStAuOqT5s8eTI5de5VpXR13QC6sr2q1K2eUx0HFWBXgX31fKVKlSKn+ntVhV69d6n2AcDuu+9OTr3fq8kyKsCujq1q99rwJ1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJMQDKGOMMcaYhBRoiHz16tVUefSoo46i7T755BNyKoyqgtcqoD137lzZnoYNG5Lr169fTu62224jp4J9KmB3xRVXkFNVUVX11FtuuSWn7VTlb0CHoitVqkROhayPPvpocips/tBDD5FTgb0bb7yRXJ8+fci98MIL5FT4UIUMVWV5ANhpp53IlShRgpw6L6ravboW27dvT27RokXk3nzzTXIqEK8mXKiqvaqSOACcfPLJ5C6//HJygwYNIvfkk0+Su+eee8hddtll+b5XFfsLKyVLlqQJBKpfUte6usdUWFldb9kqkasJF+p1Mld/AIDSpUuTU1Xp1XYqeL333nuTU0FuFSRWE0IOPfRQcoCuOq6Czeo51WSZChUqkFPVrVVAWz2fCoz/97//JXffffeRW7x4MbnKlSuTy+bVcVDXTp06dcipCSbqOIwZM4acmpCgqoarY6j6ZzXpAdD3iwqMq8lXajUL9V643XbbydfOhj+BMsYYY4xJiAdQxhhjjDEJ8QDKGGOMMSYhHkAZY4wxxiSkQEPkIQSq9Pm///2Ptrv66qvJqXCkCs1dddVV5M477zzZHlVh/KKLLiLXoEEDciqsrgLVY8eOJXfAAQeQ23///cmpSsO9evUip6oFd+rUiRwArFq1ipwKwqpqwyo0PGnSJHLq/L388svkVLvfeecdcir8/NRTT5FT4e7XX3+dHABceuml5NTx7t+/PzkV9B0wYAA5FaDOtVrw8OHDyY0fP57cXXfdRU7dA4A+PqpquZqUoIKdbdq0IZc5eaFFixayLYWRZcuW4fPPP8/nMifFALqa+KeffkpOXUeKbKFaFcieNWsWuT///JOcCgirdqvJJGpizPLly8mp60hV9FbHYfr06eQAYOHChTltq/ovdZ8VKcKfIagwvXpd9RqqwrgKjDdt2pTcc889R27+/PnkAB2KVudZVeBWAfvvv/+eXP369cmp611dIzVq1CA3ZcoUcn///XdOz5dtWzUuUCFyNbFDvXepSvlrw59AGWOMMcYkxAMoY4wxxpiEeABljDHGGJMQD6CMMcYYYxKyzhB5CGFbAGMAbJPe/oUY43UhhNoAhgCoCGAigM4xRi55m4elS5figw8+yOdeeeUV2u6CCy4gpwLHKkj39NNPk8tW1faYY44hp4K2KrB83HHHkTvssMPIdenShVz37t3JqYCdqm59ww03kFNB4DPOOIMcoCvOqpD2HXfcQa5Hjx7kVEhRBbTPPPNMcqoqtwoFquD77bffTq5z587kLrnkEnKADvI3a9aMnKqU/9dff5FT4f6RI0eSU9XXp06dSm7IkCHkVEVjFax85JFHyAE6uPrFF1+QU0Fk9di2bduS+9e//pXv+2wV8QuKjdl/FStWjCpX5xqeVoFjdR2pgLaq/A3oKtPqXlFBZFWBWwWgS5YsSU5dHyroroLqv/32G7ny5cuTyxaeVqsFqG1V1epcj626p2KM5NT5U0F1VWFcBcabN29O7uGHHyYH5F7FXp1TNWlI7YsKlqvzlzkxDNBBdTVJQb1/qHMM6GtMnRfVnl9++UU+ZyZqnLE2cvkEajmANjHGxgCaAGgbQtgXwG0A7oox7grgFwBnJ3plY4zZ9Lj/MsZsEtY5gIop1sxdLJ7+igDaAFizSNkgAO03RQONMWZ9cf9ljNlU5JSBCiEUDSFMAbAIwAgAswD8GmNc83nxPABy1dYQQtcQwoQQwgT15yNjjNmUbKz+S/1ZwhjzzyWnAVSM8e8YYxMANQD8C8Aeub5AjPHhGGOzGGOzMmXKrF8rjTFmPdlY/ZfKAxlj/rkkqkQeY/w1hDAKQAsA5UIIxdK/xdUAoFN/eV+sWDEKfpctW5a2U+HIs846i5yqwqyCxGPGjJHtOemkk8ipquUq/KbCh6qa+EcffUTu7bffJqeC3A0bNiR3/PHHkxs0aBC5Pn36kMv2eBVgv/POO8nNmTOHnApm1qpVi9wpp5xCbsKECeRU5XBVrf7YY48lpyYfqIA9ABx++OHkxo0bR06FclXw+rHHHiM3evRockcffTS5Pfbg93NVXbl3797k1EQDFRAGgBdffJHcsGHDyF188cXkOnbsSG7o0KHk+vbtm+97Vdl/c7Gh/dfKlSspWKsqWatfFFWfpgK+uVbvBnSAXQV/VShXhW9DCORU4FhN3lFB6erV+UM91WZVQVuFsQHg66+/Jqcqsqt7tHTp0uSShobz8t1335FTFa8rV65MTgXfVWA8W/+lQugq/K7O/fbbb09OhbnV9akqdatrRF3HaqUH1VepFRwAoHHjxuRWrOB5H+q9ebvttiOnjsPEiRPla2djnZ9AhRB2CCGUS/+/JIBDAEwHMArAmnfj0wFwT26MMZsR91/GmE1FLp9AVQUwKIRQFKkB13MxxtdDCNMADAkh/AfAZAD8K7gxxmxe3H8ZYzYJ6xxAxRg/BbCX8LORyhMYY8wWifsvY8ymwpXIjTHGGGMS4gGUMcYYY0xCEs3C21DULBY1E0sl+NXsi3PPPZfcnnvuSe6dd96R7VHLnQwYMICcmunUoEEDciNGjCCnZpCNGjWKnJr1dt9995FbuHAhObX8TLaZWGpmxYwZM8ipWRl77UV/CZHLu6jjpWYeduvWjZxaNufWW28lt/fee5NTswnr169PDgAOOuggcs8++yw5NbNSzYy56qqryNWoUYPchx9+SG7u3Lnk1LJAaumbJ598kpya1QcAO++8M7nMpZUAfe20bNmS3NKlS8llXrOqzYWVbbbZhmYSqSUi1OysnXbaidzuu+9OTs1SU/c8oGclqftbPaeagaRmU6nZS2qmn5p5q1D9s5p1pWb3AsC+++5LTt23art69eqRU8ufKOrUqUNOzY5Ux1XNwlMzGdXyLGq2HaD7oMcff5yc2j/VRtW3qxl8X375JTl1bNT7o5r5qWbhq2MDALvuuis51QepGZxqSSn1Ovvttx+5N954Q7YH8CdQxhhjjDGJ8QDKGGOMMSYhHkAZY4wxxiTEAyhjjDHGmIQUaIh8u+22Q9OmTfO54cOH03ZqWYBTTz2VnFoeQYV5sy0CqsLq5cuXJ6fK16uQr1qupG7duuTU8i5//fUXORXmVYE9FXyvUqUKOQA48cQTyaklBNRSMCqQqJa5UcuadOrUidyBBx5ITgXnX3vtNXIqXHnvvfeSu/zyy8lla49a/kEFr88//3xybdu2JaeWGqpWrRo5FWifNGkSObV8jVqfLdvyKWqZgkaNGpG7//77yalj88ADD5DL3L/ixYvLthRGYoy0dIQKbauQ7rfffktOTdRQwXAVVAeAefPmkct1eRG1fIraF3V/K1SY94cffiBXqVIlcsWK8duQCrQDep/Vch7qWlfBdHVsVHvU/k2bNo2c6ndVqFkdf9W/quVZAB0YV0Hw559/npy6PtVr53ps1HFV16EKpavlg1T7AB1gV9eJWn5LLa+kJi+osPna8CdQxhhjjDEJ8QDKGGOMMSYhHkAZY4wxxiTEAyhjjDHGmIQUaIj8559/pvCbqgb6xBNPkPvss8/IqbCyCmirCq+ADnEOHDiQnKoSfskll5BTFca7du1K7qijjiLXrFkzcqqStQqlqwC6CvEBuiJyixYtyKnAX8eOHcnVrl2bXK9evcip4zB9+nRy55xzDrmRI0eSU8fwlltuIZetErkKYapjpkKhqpq4CoWqAPuQIUNyel21f6qqffPmzcmpID4AjBs3jpy6ntQ92bp1a3IqdJxZpVpVwS6srFixgs6VCrGqQK4KTxctWpTcV199RU4FpwGgatWq5H799VdyKois7m8VdlZ9iwq1q2tBVZ7+8ccfyanjpfokQIeO1T6rvlMF4tX+qerrKnCsXlcdQzXZSb33qGOozhOgj5kKjKuJOup9Sh0HNeFF9dkqWD5r1ixyqnK76mvUygyArkSujpk6L2pymNpO3ZNrw59AGWOMMcYkxAMoY4wxxpiEeABljDHGGJMQD6CMMcYYYxJSoCHyxYsXY+jQofmcCrllbgMA//73v8n17t2b3CeffELujjvukO3p378/ORUQvu2228hdeeWV5K6++mpyqjq5qkqswocqcKnap1y2EKaqoH7AAQeQU1W0VQhQBRJV4FiFlXv27ElOHWv1uioQvXDhQnLHHHMMOUCHr1XFevX4448/npwKpe+///7kJk+eTE7tX6tWrcgNHjyYnLqGr732WnIAsPvuu5NTx1tNaHj33XfJqdBqZsXfo48+WralMFKsWDE6V6oKtgrGqhCy2q569erksgVbVTC2YcOG5NSKC5UrVyanwuaqkrx6jdmzZ5NTFaFVm9VEDXV/AsA777yTU3tUcFu157fffsupjarif+aECQCoVasWObV6hOpr1ISEbFW5lVdO9c/qffOhhx4ip64Hdc2q1RpU5W/lVLXzbKhrVk1SUU4dG9XubJXfs+FPoIwxxhhjEuIBlDHGGGNMQjyAMsYYY4xJiAdQxhhjjDEJKdAQeZEiRSjMNWLECNpOBQVPP/10clOnTiWnKjj37dtXtkdVqX777bfJPfXUU+Tuvvtuci+88AK5evXqkWvQoAG5W2+9lZwKHA8fPpzcnDlzyP33v/8lB+hQrwoxL168mJyqAqwmAcyfPz8nd9BBB5FTAe1BgwaRe/PNN8l16NCBnDo2ADBz5kxyKrg6duxYcmeeeSY5FdDu0qULOTWhQe3zaaedRm716tXkunXrRu7www8nBwA33HADuSeffJLckUceSU6dA3V9Zu6LqlJcWFmxYgVVBS9ShH8HVZNEVIhVVa9XAVj1GgBQtmxZcpMmTcrp8SooPWbMGHJqMoqq4KyqR6v9U6F7FageP348OUBXQZ8xYwY5FRBWIXJVnVxVg1fHWqEqev/+++/kVABdPVZVLAeAvfbai5yqTq6eUwXG1SQtdQ5yrQSvXlcF+1XVd9XPAcCiRYvIqeC9aqOa5KD6JrVSx9rwJ1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJMQDKGOMMcaYhBRownPlypUUNuzRowdt95///Cen51PB5AcffJDcsGHD5OPvvfdecipoqwLsKnitAmhTpkwhp/avU6dO5FRgfL/99iP3ww8/5OQA4LXXXiPXpk0bct27dyd30003kWvXrh25p59+mpwKK6vAf9euXcl17NiRnAriq3C3CvYDOlw+evRochdccAG5zp07k3vvvffIHXzwweRUmPiZZ54hp6ruqqrJKsD51ltvkQOAWbNmkatQoQI5FSZu2rQpuVKlSpG75ppr8n2v7tHCSvHixVGtWrV8ToVYq1atSk5VcI4xklPBVhWczratquqtrhG18kHmvmVro2K33XYj99lnn5FT4XUVNq9Zs6Z8HVXNOjPYD+ggcuPGjcmp8Lta0UD19+rcq6C0Ok9qApMKlq9cuZIcAJQoUYKcquqtzqmqtq0C423btiWn+k1VKV/1K6p96jiovhTQ4fBKlSqRU9esmgSg7slsld+z4U+gjDHGGGMS4gGUMcYYY0xCPIAyxhhjjEmIB1DGGGOMMQkp0BB5lSpVcOGFF+Zzl19+OW2nqmirSqmqgu2ee+5JTgWBAaB3797kBg4cSE6FIZs0aUJOBaVVFXNVBVsFmNWxefjhh8lNmzaNnArIAcBLL71EToXudtllF3KqinCdOnXIHXLIIeT22WcfcirsrCopT5w4kVz16tXJNWzYkJwKtwLAp59+Su6AAw4gp87p888/T04dwwceeICcCuefcMIJ5O677z5y/fv3J/fqq6+SU5XNAaBVq1bkXnzxRXLq3KvJFSoQfN555+X7XlWML6ysWLEC3333XT6nJmuo6tZqxQV1DavAcbZAtZoUoFYQUCFydd2o11HPp0Ltqir9DjvsQO7LL78kp46hqqoN6D5IBcFVUFqtFqCC4Kp/VhWzf/zxR3IqnK8mZSxdupScmiSS7dyr46iunenTp5NT7w0qoK0C46qPVNeD6l9V2FxVyVfHFdCTM37++WdyalKCCoerSTDq+lob/gTKGGOMMSYhHkAZY4wxxiTEAyhjjDHGmIR4AGWMMcYYk5ACDZEXKVKEglvHHHMMbVe3bl1yqjqsCtUqVMVxIFVZOJORI0eSO/nkk8mpKuFqOxX2U6HHRo0akatRowa5li1bklPVXAcPHkwOAE466SRyKiitwnndunUjd+mll+bURhWiVcfm448/JqeOjarmrqpeq0AhoIO1mQFhQIdRDz30UHKqQvLNN99MTl1ztWvXJqdCnSo0rML52arQq6r/qjK9Cow/9NBD5FR4dNSoUfm+V5WZCyvbb789WrRokc998803tN1ff/1FTp0TFSJv3rw5OXXeAV0dXgVtv/jiC3Kq0r0KVKtJOWqix957701OXetq4oEKd6v2Afq+XbBgATnV7lyrf6vt1GQZhQom77rrruRU9XR1vNQ5ydYe9doqWL7zzjuTU8dVhb7VPa/eF66//npyKsitJvnsscce5ACgVq1a5FQ/qarGq3tSTQRT1dLXhj+BMsYYY4xJiAdQxhhjjDEJ8QDKGGOMMSYhHkAZY4wxxiQk5xB5CKEogAkA5scYjwoh1AYwBEBFABMBdI4xrrWMZ8mSJalatKqEfNlll5Hr2bMnufbt25NTIbdnnnlGtueNN94gl2vgUgW8VbBPVVVVwWsVkr/yyivJdenShZw6NqoKMKCDq+eeey45FQzs1KkTuRdeeIHcfvvtR05VeD3jjDPIffXVV+Qyq9cDwJAhQ8gdffTR5FRVewB4++23yakJBKrC8j333ENOXQ+qirC6Hm688UZyqmK8CldmhrYB4NtvvyWX7XVUFfRXXnmFnAqRn3jiieT233//fN+ryvmbg43Rf61YsYJCxyp8u3r1anLqWqhQoQI5FWpWFaoBvRKDCqurx1epUiUnN378eHIqDKxC0eo4VKpUiZyaqKEqVAO6IriahKHOiwrjqzaq4LU6LytXriSnguAq6Kz6C3U+1SQFAGjQoAE5NQlG7bOa8KJC1ur6VBXGVWBcTVh59NFHyamQfLYq9KrCf7169cjNnj2b3KJFi8ipALt6jbWR5BOoiwHkrQt/G4C7Yoy7AvgFwNmJXtkYYwoO91/GmI1KTgOoEEINAEcCeDT9fQDQBsCajx8GAWi/CdpnjDEbhPsvY8ymINdPoO4GcAWANZ93VgTwa4xxVfr7eQD4b0MAQghdQwgTQggTsn0UbYwxm5C7sRH6L/VnDmPMP5d1DqBCCEcBWBRjnLg+LxBjfDjG2CzG2Kxs2bLr8xTGGLNebMz+K2mRPWPM1k0uIfKWAI4JIRwBYFsAZQD0B1AuhFAs/VtcDQBcBjqDpUuX4v3338/n3nzzTdru2muvJXfccceRUwHfXXbZhVyHDh1ke1RgrFWrVuRU4E+FCkePHk1OBctVMFNVNm/dujU5FfjOrI4MAM8//zw5QIeqVQjwiiuuIKeCmf/73//IqcrCp59+Ojl1/Hv16kXulltuIacqNqvrQYXSAT2xQAVcVSVzdWxvvfVWch999BE5VeFdTaRQwfk777yTnLqWVDgfAO644w5yV199NTkV0FfH4ZFHHiGXOUFCBW0LmI3WfxUtWpQCuCrwqiYKqEkdamKFGqSpkDWgq7yrbdV9psLTc+bMIVetWjVyv//+e05OVYReunQpORUMVw7QIe0ZM2aQK1euHLkSJUrk1B71lxLVN6gK6Oo8q3D+lClTyKnzpILOABBjJKf2T50DFSJX15K6FlV1cjWpQAXG1aSTXN/DAf2eq869Cu2r86cm22S77rKxzk+gYoxXxRhrxBh3BnASgJExxlMAjAJwfHqz0wG8muiVjTFmE+P+yxizqdiQOlBXAugZQvgaqUzBYxunScYYs8lx/2WM2SASLSYcYxwNYHT6/7MB/GvjN8kYYzY+7r+MMRsTVyI3xhhjjElIok+gNpSyZcviqKOOyueaNWtG26kqpqoCt6pCqp5PhZABHe4bMGAAuZYtW5IbOHAgucMOO4zcc889R65q1ark2rZtS+6mm24iV7p0aXLXXXcdufPPP58coI+PCt1lhv0BoHfv3uQ+//xzcmq6t6r8rY6hqsY7ZswYcirUPHEiT7Q6/PDDyQE6bHvRRReRU4FqtZ0Ktb/++uvk1ESD22+/XbYxExXCVMFKFV4HdNjz3XffJaeC/GpigKrQnznT9sgjj5RtKYxsu+22VA1eHVM12UKFftV9p4K7qr8AdGV6FZZV4XA1I1r1Lbvuuis5VR1btVtNIFB9rlo9Ql1vgK62rfpnVclcVe9WVcdVsLx+/fo5bacC2mr1B9X/qFC6qgYO6POnrqcdd9yRnDpXKjyt9lkdVxV0VxXGVWBcrTrSuXNncgDQuHFjcup6Un2iOg5qgkS2CvjZ8CdQxhhjjDEJ8QDKGGOMMSYhHkAZY4wxxiTEAyhjjDHGmIQEFW7cVJQoUSJmhr4eeugh2m7+fC4KrCpZn3nmmeRUsFK9BqCr56rwoqoYvHDhQnIvvPACuTZt2pCrVasWORUuVvuiwvQ33HADuWyVyI855hhyKuyswsUqHN69e3dyKui5atUqcvfeey85VQVYhSNVdXh1vMqXL08O0MFa1R4Van/ggQfI3X333eRGjRpFTk00UGFs1T4VeFUV0MeNG0cOAA444ICcnLpmR4wYQW7mzJnkMgPsQ4cOxaJFi/gEFkLKlCkTMydhqOtVVflWqL63ZMmS5FS1c0CHalXfuXjxYnKqYrZ67J577klOVXBWTt3ze++9NzlV3VpV5QZ0xewPP/yQnKrA3aRJE3IqPK0mo6j3hWnTppFTfbuqrD1v3jxyqpK4WrkA0KFoNaFh7ty58vGZqGtWBapVlX3V76prVh0Htc/Zzv3HH38sfSZqUsJOO+1ETr2vq+tm+PDhE2OMPPsK/gTKGGOMMSYxHkAZY4wxxiTEAyhjjDHGmIR4AGWMMcYYkxAPoIwxxhhjElKgS7mEECjZ/+yzz9J2araXKoffvn17ciNHjiTXtGlT2R61zIdamkTNUFAl5Dt27EjumWeeIXfaaaeR++STT8ip2SlqGQU1+0/NgMnG2WefTU7NcFQzd9QyK4MGDSKn9nns2LHkLr74YnLqWNepU4ecmimjzicALF++nNzJJ59MbtGiReTUTKB27drl9NiePXuSe+KJJ8idddZZ5NSSFWpmnpptBOjlI+rWrUtOzSR65ZVXcnKdOnXK933x4sVlWwoj22yzDc0kUsf0q6++IqfuZXWe1MzgbLRo0YKcWoJJLQWjZiUddNBB5NTMPNUXN2zYkJyaSapmbKllP9R1CeglnEqVKkVOHe/99tuPnJqlpvpO1Qep2cGqD1Kzz9R1o5YAUjMtAb3EzpdffpnTdmrWr+qr1P6pa0n1F2omnVpiRe2fuoYBPVbo168fOTV+UNdxtmW+Mhk+fHjWn/kTKGOMMcaYhHgAZYwxxhiTEA+gjDHGGGMS4gGUMcYYY0xCCjREvv3226NVq1b5nAoAqpCbWq5ElbNXQUi1tAigQ3LDhg0jp4LNgwcPJqdCfOo1VND9wgsvJKeCi926dSOngn3jx48nBwBPPfUUuQ4dOpA79thjyallGPbaay9yKoCuAtD3338/ubvuuoucWn5GhW0zl9kAdLgV0Of0ueeeI6cCiWrygno+tTyPOjYXXXQRuX//+9/kPv30U3Jqn9UxBPQ99PXXX+f02n369CH3yCOPkNt///3zff/kk0/KthRGVqxYQWHUSpUq0XZqCQw1aUEtA6OcCk4DeokWtTSJWuJDLb1Srlw5cur6X7p0KTm1fyoUrZYmUa+RjW233ZbcypUryakQuQr3q+dTy5WULVuWnJrAod5rdt55Z3LqulFL+6jlRgB9DlQYXx0bda7UY1Ub1flT73H16tUjN2PGDHLqPS4bKjCulqJS783qOKoJWWqJsLXhT6CMMcYYYxLiAZQxxhhjTEI8gDLGGGOMSYgHUMYYY4wxCSnQEHm1atVw3XXX5XMqDNe1a1dyvXv3JnfeeeeR+89//kPusMMOk+0ZNWoUuR49epBTwU4VflfVqFWVcBXiU9VvP/74Y3KzZs0iN3ToUHL33HMPOSD3cPI555xDToWsP/jgA3KHHnoouSOPPJKcqlKtgqfly5cnp46rqjZ76aWXkgP0uVcVxgcMGEBOVW5X26nAvqqg/ssvv5Dbd999yalqwapq77x588gBwLRp08ip8zJixAhyp556KrmffvqJXGYoXV3rhZVixYqhYsWK+ZyqrpxruFtVwS5dujQ5FdwFdOhYnZNly5aRUxXB1XWo7gnVV6njoCYDKVRVbnVssqFC3yocrvb5hx9+IKfO1bhx48iVLFmSnAqMq8C+cqo/VPsG6MC/Oi+qP1X3pAqCq+tLTWjIvCcAYPbs2Tm9hrpu1HUI6MkGKjCuguUvv/wyOXUMVR+7NvwJlDHGGGNMQjyAMsYYY4xJiAdQxhhjjDEJ8QDKGGOMMSYhBRoiX7hwIYWbVeB4l112IafCdKra7OOPP06uZ8+esj0tWrQgp0K1qkrrww8/TE6FD1VYWe3zgQceSE6F+K688kpyKhiuQu6ArhyuqkX379+f3PDhw3N67f3224+cCqWrcL4KWavq8Ntssw05FTx94403yAG6Cu0ll1ySU3tUBXwV4Bw7diw5FYRUVXLV86nq8KqSeK1atcgBwMyZM8mpyrsq/P7oo4+SU2HpzJUAkoYyt2RWrVpF/ZDql6pXry4fm8nChQvJqcD3DjvsINujgs2qT1Ru8uTJ5FToe/Xq1eTmzJlDTl3XKqCtVgZQIfJs+6yqiavq32p1BhVsVkFwNbGpefPm5FQfsmDBAnKqv1AhdxWeVu8pgL6Xy5QpQ05NwFGV6dWxUdXX1WPVuVf3vapsro7NTjvtRA7Qk4RU36kC4126dCE3depUckmq4gP+BMoYY4wxJjEeQBljjDHGJMQDKGOMMcaYhHgAZYwxxhiTkAINkVeqVImqOGeGTgHg1ltvJadCj6oC6hdffEHu2Wefle1RITIVYFOBxgcffJBcvXr1yB1++OE5tUdVyVXVrXfffXdyFSpUIHf99deTA3Q16vr165NTVVqHDBlCTgXBVXVydRzUvqi2tGzZkpyqVq/2TQV1Ad1udRwbNWpETk1yUJMFWrVqRU4Fy9u0aUNOBZEV559/PrnOnTvLbVVQVIVZ1USFBx54gFzdunXJZVZB32effWRbCiMlS5akELSazKCuDxU4Vte/Cm0XK6a76aZNm5JTVabVSgOqv1EB2lyrpa9cuTKntqhQs6purSY3ALpat5owo4LpagKNCm6rYLk6B6rqu7qf1KoCkyZNIqdQoW1Ah69VxXpViVxVS1f7p15DHesmTZqQ22233cip6uvVqlUjp4LhgH4PUfeVeu9S7/WXXXYZObXCxdrwJ1DGGGOMMQnxAMoYY4wxJiEeQBljjDHGJMQDKGOMMcaYhBRoiLxYsWJUjVRVv61SpQq5QYMGkfvoo4/IqUBuu3btZHtUEPzaa68lN2bMGHJ//PEHOVX1+u677yanwuaqmm5mIBfQ1X1VqHPkyJHkAKBt27bkVPBeVQzef//9yfXt25ecqpKrQtsqoK0qsl911VXkVNBZhXf79OlDDgCuueYacgMHDiR34oknklOh19tuu42cCsyqKrkdO3Ykd/vtt5NT94UKwaqq7wDQo0cPcqpysgqmH3300eQ6depE7qyzzsr3vQp5FlaKFi1K1Zlr165N2+W6zyrorIK72SZCqICxqgCtgr/qOVXVazVRRwW5VXVydZ+o4LwKOqsgcDavJnCogL06V99//z05dU/suuuu5Bo0aEBOHRt1XGvWrElOVU+fOHEiOUCv9qD6XVWFXl13avKUur5UwFsF59XkCnUcihThz3DUNZcNNdlA3QNqgoQKjF9wwQXk1PvPGvwJlDHGGGNMQjyAMsYYY4xJiAdQxhhjjDEJ8QDKGGOMMSYhIcZYcC8Wwo8AvgVQCQCXcS2ceF+2TLaWfSns+7FTjJFL+RdC8vRfQOE/L2vYWvYD8L5sqRT2fcnahxXoAOr/v2gIE2KMzQr8hTcB3pctk61lX7aW/dja2FrOy9ayH4D3ZUtla9qXTPwnPGOMMcaYhHgAZYwxxhiTkM01gHp4M73upsD7smWytezL1rIfWxtby3nZWvYD8L5sqWxN+5KPzZKBMsYYY4wpzPhPeMYYY4wxCfEAyhhjjDEmIQU+gAohtA0hzAghfB1C6FXQr78hhBAGhhAWhRA+z+MqhBBGhBC+Sv/LK8huYYQQaoYQRoUQpoUQvgghXJz2hXFftg0hfBJCmJrelxvSvnYI4eP0dTY0hFBic7c1F0IIRUMIk0MIr6e/L5T7sbXi/mvLYGvpw7a2/gv4Z/VhBTqACiEUBXA/gMMB1AfQKYRQvyDbsIE8AaBthusF4L0Y424A3kt/v6WzCsClMcb6APYF0C19HgrjviwH0CbG2BhAEwBtQwj7ArgNwF0xxl0B/ALg7M3XxERcDGB6nu8L635sdbj/2qLYWvqwra3/Av5BfVhBfwL1LwBfxxhnxxhXABgCoF0Bt2G9iTGOAbA4Q7cDMCj9/0EA2hdkm9aHGOOCGOOk9P+XIHWxV0fh3JcYY1ya/rZ4+isCaAPghbQvFPsSQqgB4EgAj6a/DyiE+7EV4/5rC2Fr6cO2pv4L+Of1YQU9gKoOYG6e7+elXWFmxxjjgvT/fwCw4+ZsTFJCCDsD2AvAxyik+5L+yHgKgEUARgCYBeDXGOOq9CaF5Tq7G8AVAFanv6+IwrkfWyvuv7ZACnsfthX1X8A/rA9ziHwjElM1IQpNXYgQwvYAXgTQI8b4e96fFaZ9iTH+HWNsAqAGUp8S7LF5W5ScEMJRABbFGCdu7raYfyaF6Z5fw9bQh20N/Rfwz+zDihXw680HUDPP9zXSrjCzMIRQNca4IIRQFanfIrZ4QgjFkep4no4xvpTWhXJf1hBj/DWEMApACwDlQgjF0r/5FIbrrCWAY0IIRwDYFkAZAP1R+PZja8b91xbE1taHFfL+C/gH9mEF/QnUeAC7pVP5JQCcBGBYAbdhYzMMwOnp/58O4NXN2JacSP9d+jEA02OM/fL8qDDuyw4hhHLp/5cEcAhSeYhRAI5Pb7bF70uM8aoYY40Y485I3RcjY4ynoJDtx1aO+68thK2lD9ta+i/gH9qHxRgL9AvAEQBmIvV33qsL+vU3sO3PAlgAYCVSf8s9G6m/8b4H4CsA7wKosLnbmcN+tELqo+1PAUxJfx1RSPelEYDJ6X35HMC1aV8HwCcAvgbwPIBtNndbE+xTawCvF/b92Bq/3H9tGV9bSx+2NfZf6fb/I/owL+VijDHGGJMQh8iNMcYYYxLiAZQxxhhjTEI8gDLGGGOMSYgHUMYYY4wxCfEAyhhjjDEmIR5AGWOMMcYkxAMoY4wxxpiE/D81px1Fy0FDAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.random.random((50,50))\n",
    "ATA = A.T @ A\n",
    "\n",
    "fig1, axes = plt.subplots(1, 2, figsize=(10,6))\n",
    "axes[0].imshow(A, cmap='gray')\n",
    "axes[0].set_title('$A$ is Random Matrix')\n",
    "axes[1].imshow(ATA, cmap='gray')\n",
    "axes[1].set_title('$A^T A$ is Symmetric Matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42dcc460",
   "metadata": {},
   "source": [
    "## 5.5 Matrix zoology (cont'd)\n",
    "\n",
    "### Diagonal Matrix\n",
    "A matrix often written as $D$ with zeros along the non-diagonals.\n",
    "\n",
    "$D$ can be square (top) ...\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "d_1 & 0 & 0 \\\\\n",
    "0 & d_2 & 0 \\\\\n",
    "0 & 0 & d_3 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "... or non-square (bottom).\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "d_1 & 0 & 0 & 0 & 0 \\\\\n",
    "0 & d_2 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & d_3 & 0 & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Transforming a matrix to diagonal form can be achieved through eigendecomposition or singlular value decomposition.\n",
    "* When diagonal elements all have the same value, then often written as scalar multiple of identity matrix $\\lambda I$.\n",
    "* Hollow matrices are opposite of diagonal matrix and have zeros on diagonals and non-zero off-diagonals.\n",
    "    * Euclidean distance matrices are often hollow since self-distance is 0.\n",
    "\n",
    "### Triangular Matrix\n",
    "A square or rectangular matrix containing non-zero values along the diagonals and upper triangle $U$ or lower triangle $L$.\n",
    "\n",
    "$$\n",
    "U =\n",
    "\\begin{bmatrix}\n",
    "a & b & c \\\\\n",
    "0 & d & e \\\\\n",
    "0 & 0 & f \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "L =\n",
    "\\begin{bmatrix}\n",
    "a & 0 & 0 \\\\\n",
    "b & c & 0 \\\\\n",
    "d & e & f \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "a, b, c, d, e, f \\neq 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1b3bcd41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEvCAYAAACKfv/MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkaklEQVR4nO3df7xVc9738fdHpXL6oWJSMdJISkgSujF+zETiJmHEUEO6HoRpZlwjhRpug3vGhWbCVCJuFaHBNSYiEtGouRo6dSE06HeEFCZ87j/2znUcnfbne84+5+yzez0fjx5nn7Xfe63vOvvsT5+z9trfZe4uAAAAxO1Q2wMAAACoa2igAAAAEtFAAQAAJKKBAgAASEQDBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UHWAmZWa2dG1PY6aVij7bWbLzOxHtT0OoC6qzOu4Lr/mCqFu1eWfX11CA1ULzOwIM/vYzKzc8mfM7IryeXffz92fq7EBVrPoi7vY9hsoJma2i5m5mbXbVq5YXsfULZRHA1U7DpK00L97HZ1ukv5e88OpPmZWvyYeU2iKYR+AHLpJWuvuy2t7IPlU2dduMbzmi2EfahINVO04SNJ/lV1gZu0ltSy/PHvfN3/5mNkVZrbczDaY2etmdtzWNpD9y3DvMt/fY2b/p9w6rzSzxWa23szuNrNGkfvNrK2ZPWxma83sHTO7bCvjvcLMXpW0seyL0szuk/R9SY+b2adm9uuKHlNuv4eb2VvZ/V5sZv22ss3LzezV7NG9B8qMt7uZ/Vf2sdOy95X9WWzzZ1VuOxWOY1v7DRShbpL+kStUmfq1lXV0NrPnzOyj7Ftk/zu7/Gdm9niZ3JtmNq3M9++ZWbfs7QrrVq7XbnXUrW3VrOz91K0CRwNVO7rpu0eaDpL0vruvq+hBZtZJ0iWSDnH3ppKOl7SsCuM4J7uOH0jaR9JVue43sx0kPa5M4Wwn6ThJw8zs+HKPHSCpr6Sd3f3LLQvd/VxJ70o62d2buPv/zfWYrLckHSmpuaTfSPp/ZtamXOZMSSdI2kvSAZIGmdmOkqZLukeZBnWKpH6qvFzj2NY+AMWkm6SF0XBl65eZNVCm5jwl6XuSLpV0f3Z9syUdaWY7mFlbSTtKOjz7uA6Smkh6NVi3KnztVmPd+k7Nyo6dulUH0EDVsGwx2E/fPdLUXbnfvvtKUkNJXcysgbsvc/e3qjCcP7r7e+7+oaTrlXkR5br/EEm7uvu17v4vd39b0nhJZ5V77JjsYz9LGE+Fj3H3ae6+wt2/dvcHJL0pqedWHr8iO97HlSnwh0mqn71vs7s/IulvCWNKHUdl9huoi7opoYFS5evXYco0Qjdma84sSf8paUC2/mzIjuUoSU9KWmFm+0r6oaQ57v61YnWrsq/dqtStrdWsLftM3SpwNFA1r3P265Jyy3tKemVbD3T3pZKGSRotaY2ZTc3+1VVZ75W5/U9J5de1tfv3lNQ2eyj9IzP7SNIISa238djKjOdbzOw8M1tYZptdJe1SLraqzO1NyhTdtpKWlzvfrDJji46j0usG6orsW02dtJVTDipShfrVVtJ72UZoi38qcyRJyhyFOlqZBmq2pOeUaZ5+mP1eitWtyr52q1K3tlazJOpWnUADVfNaS/qs7GFSM9tFmRf7Y7ke7O6T3f0IZQqCS7qpgugmSTuV+X63rWT2KHP7+5JWBO5/T9I77r5zmX9N3f3E8kPd1m6kLDezPZX5a/ESSa3cfWdJiyTZ1vLlrJTUzuxbn3jco1wm8rOKjmNb+w0Ui66SNkt6PeVBCfWrrBWS9si+DbfF9yVtOXl9SwN1ZPb2bH23gYrUrVyvXeoWvoUGqua9JqmhmV1oZo3NbB9l3t/+T3d/dVsPNLNOZnasmTWU9LmkzyR9XUF8oaSzzayemZ2gTDEpb6iZ7W5mLSWNlPRA4P6/SdqQPemwcXb9Xc3skNjuS5JWS+qQkC9R5gW+VsqcOKpMAY94SZm3Di7JnuB5ir771t9C5f5ZVXUcQDHpJmmxpAZm1ij7r+G2HpBYv8qap0yz8Gsza2CZOZZOljQ1e/9sScdIauzu70uao8x5Ra30P0fIqFvUrbyjgaph7r5KmRMHL5P0oaQZypz7NCjw8IaSbpS0TplDv9+TdGUF2Z8rU2Q+UuZk8D9vJTNZmRMz31bmJMPyn+D4zv3u/pWkk5QpoO9kxzJBmZMTo25Q5oT0j8zs8lxhd18s6WZlispqSftLejGyIXf/l6TTJF2gzM/ip8qcP/FFmVjkZ1WlcQBFppukg5Vpgrb8W5zjMSn16xvZ1/DJkvpkH3u7pPPc/b+z978h6VNlGie5+yfK1KwXs/VK1C3qVnWw705FhO2BmS2TNNjdn67M/XWZmc2TdKe7313bYwGACOpW4eEIFIqemf3QzHbLHgofqMzHhWfU9rgAoCLUrcLHhFnYHnSS9KAy5wK8Lel0d19Zu0MCgG2ibhU43sIDAABIxFt4AAAAiWigAAAAEtXoOVAlJSW+884758ytWrUqZ0aSdtxxx1CubduqTNa9dStXxt6K7tSpUyi3eHGuTwBnRH5+klRSUhLKSVLDhtucvuUbK1aUn2dz61q1ahXKff7556HcLruUn3B86957Lz6ZbvR3Yt26Ci9N+C1t2pS/LN/WLVq0KJRr0KBBXrcb/VlL0ocffhjORnzxxRfr3H3XvK60FpSUlHiLFi1y5pYvX54zA6DOqLB+VamByk7edZukepImuPuN28rvvPPOGjp0aM71Xn/99aHt77XXXqHc6NGjQzlJ+vrryLxu0g033BDKzZo1K5Tr0aNHKNe3b99QrlevXqGcJO25556h3HXXXRfKnX322aHc0qVLQ7nBgweHcpdddlnuUNa1114byo0bNy6UGzVqVCjXsWPHUK5du3a5Q5JGjhwZypWWloZykjR16tTcIUn16tUL5V5//fV/hjdew1JqWIsWLUL1a8SIEfkbIIDaVmH9qvRbeGZWT9JYZSY36yJpgJl1qez6AKAmUcMAVEVVzoHqKWmpu7+dnTV1qqRT8jMsAKh21DAAlVaVBqqdvn0F5/f1P1fHBoBCRw0DUGnV/ik8MxtiZvPNbP7GjRure3MAkDfULwAVqUoDtVzSHmW+3z277FvcfZy793D3HimfDAOAapazhlG/AFSkKg3UK5I6mtleZrajpLMkPZafYQFAtaOGAai0Sk9j4O5fmtklkp5U5iPAE909/llpAKhF1DAAVVGleaDc/QlJT+RpLABQo6hhACqrRmciX7FiRWjiv0mTJoXW9/LLL4dy0UkvJWnYsGGh3LJly0K5N954I5SLTvb5+OOPh3IDBw4M5aT4zzE6E/mGDRtCuYMPPjiUO/TQQ0O5xo0bh3KStHr16lBuzJgxoVx0hvjoBKf9+/cP5Q455JBQbsKECaGcJN144zbnw/3G7NmzQ7nXX389vO1C1qhRI3XpknuaqHvvvTe0vvPOO6+qQwJQi7gWHgAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJCIBgoAACARDRQAAECiGp2JPCo6i/a4ceNCub/+9a/hbdevH/uR9OnTJ5R7+OGHQ7nTTz89lGvevHko16hRo1BOknbddddQLjrTd3Qm66uvvjqU69y5cyj31ltvhXKS9Nxzz4Vy69evD+Wiz/PMmTNDuegM6Pvvv38oF50dXpJGjRoVykX35Q9/+EN424Vsp512Uvfu3XPmevbsGVrf+PHjQ7kLL7wwlANQszgCBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJCoIGcijxoyZEgoN3/+/PA6L7/88lDu6aefDuWmTJkSyg0dOjSUKy0tDeUefPDBUE6SXnrppVDujDPOCOWiM4f/6Ec/CuWiz8mtt94ayknxMc6bNy+UW7RoUSg3fPjwUG7hwoWh3P333x/K7bTTTqGcJJ155pmhXK9evcLrLAYbN27U3Llzc+ZWrFgRWl/0Nb9gwYJQTpIOPvjgcBZA1XAECgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJCIBgoAACARDRQAAECiOn0pl6gePXqEs9HLlVx66aWh3JtvvhnKzZ49O5RbvHhxKNeqVatQTpK6d+8eynXp0iWUu+SSS0K5pk2bhnIdO3YM5SZMmBDKSZnLckT07ds3lDv22GNDuWuuuSaU+9Of/hTK9e7dO5Rr0KBBKCfFf3duueWWUK5fv37hbReyr776Sp988knO3LvvvhtaX/S5e/7550M5KX4ZnpRLPQHYOo5AAQAAJKKBAgAASEQDBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACQyd6+5jZnV3Maq2eDBg0O56GzDpaWlodz1118fyvXq1SuUk6TLLrsslLvqqqtCuebNm4dy55xzTij32muvhXLR2bslacyYMaHc6tWrQ7lBgwaFckcffXQo16FDh1Bu7NixodwNN9wQyklSy5YtQ7mXXnoplFu4cOECd49fDqBA7bDDDl6/fu6LN0yePDm0vvHjx4dyJSUloZwkjRgxIpQbN25cKBcdI1DEKqxfVbqUi5ktk7RB0leSviyGIglg+0ENA1BZ+bgW3jHuvi4P6wGA2kANA5CMc6AAAAASVbWBcklPmdkCMxuSjwEBQA2ihgGolKq+hXeEuy83s+9Jmmlm/+3uz5cNZIsShQlAIdpmDaN+AahIlY5Aufvy7Nc1kqZL6rmVzDh378HJmQAKTa4aVrZ+mVltDBFAgap0A2VmJWbWdMttSb0lLcrXwACgOlHDAFRFVd7Cay1pevavsvqSJrv7jLyMCgCqHzUMQKVVuoFy97clHZjHsQBAjaGGAaiKfMwDtV2aMGFCKBedyTo6e/fEiRNDuRNPPDGUk6Rp06aFctEx9u3bN5TbaaedQrkVK1aEcrNmzQrlJOmmm24K5ZYtWxbKPf7446HcsGHDQrlf//rXodzixYtDuSeffDKUk6SuXbuGcitXrgyvsxjsu+++uvfee3PmolcLiM5Yfu6554ZykjRy5MhQ7rbbbgvlVq1aFcpFf/+BYsI8UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJCIBgoAACARDRQAAEAiGigAAIBENFAAAACJmIm8mkVn6G3SpEko9/LLL4dy1157bSgnSc2aNQvl2rZtG8r17t07lGvatGko17Bhw1Bu6dKloZwUnwX9mGOOCeXmzp0byrVq1SqU++lPfxrKde7cOZQbMWJEKCdJhx9+eCj35ptvhnLnn39+eNuF7KuvvtLGjRtz5tavXx9aX3Tm/CVLloRyknTmmWeGcjfccEMo16JFi1AuWr8+/fTTUA6oCzgCBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJCImcgLRHSG3uis4R999FF42wsXLgzlJk+eHMoNGjQovO2IRYsWhXL16tULr/Pjjz8O5f74xz+GctOnTw/ldt5551Cuf//+odw111wTyrVp0yaUk6RHHnkklDvggAPC6ywG9evXD83MPXLkyND6ojP2z5w5M5STpIkTJ4Zy0d/XU089Na/rGz58eCi3YMGCUA6oTRyBAgAASEQDBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhk7l5zGzOruY1t5zZs2BDO3nHHHaHckiVLQrn169eHcmPHjg3lOnTokNf1SZKZhXJPPfVUKBedLTq63aOOOiqUO+igg0K5vn37hnKS9K9//SuUO+SQQ0K5+vXrL3D3HuEBFKiuXbv6Qw89lDN36623hta3zz77hHKXXnppKCdlZkuP6NOnTyj35ZdfhnLR10m/fv1CuSlTpoRyklRSUhLOApVQYf3iCBQAAEAiGigAAIBENFAAAACJaKAAAAAS0UABAAAkooECAABIRAMFAACQiAYKAAAgEQ0UAABAIhooAACARLF5/1HnNG3aNJydOnVqKPfss8+GcnfddVcod/jhh4dy3bp1C+UeffTRUE6SPvjgg1AuehmZt956K5Rr27ZtKDd8+PBQ7i9/+Uso98QTT4RyknT++eeHcpdffnl4ncVghx12UKNGjXLm7r333tD6or+D+++/fygnSU8++WQod9JJJ4Vyn332WSjXv3//UK5Vq1ah3KBBg0I5SXrwwQdDuTPPPDO8TiAi5xEoM5toZmvMbFGZZS3NbKaZvZn92qJ6hwkAlUMNA1AdIm/h3SPphHLLhkt6xt07Snom+z0AFKJ7RA0DkGc5Gyh3f17Sh+UWnyJpUvb2JEmn5ndYAJAf1DAA1aGyJ5G3dveV2durJLXO03gAoCZQwwBUSZU/hefuLskrut/MhpjZfDObX9VtAUC+bauGla1f0ZO+AWwfKttArTazNpKU/bqmoqC7j3P3Hu7eo5LbAoB8C9WwsvUr+gkyANuHyjZQj0kamL09UFL88+MAUPuoYQCqJDKNwRRJL0nqZGbvm9kFkm6U9GMze1PSj7LfA0DBoYYBqA45J9J09wEV3HVcnscCAHlHDQNQHZiJHDrrrLNCuXnz5oVy0VmEDz300FBu7dq1odxuu+0WyknSu+++G8odeOCBoVx0n3fddddQbtOmTaFc8+bNQzkzC+Uk6aqrrgrlevfuHV5nMfjss8+0ZMmSnLmWLVuG1nf33XeHcvfdd18oJ0kzZ84M5WbNmhXKRWebb9++fSjXs2fPvOak+Ozrf/vb3/K+bWzfuBYeAABAIhooAACARDRQAAAAiWigAAAAEtFAAQAAJKKBAgAASEQDBQAAkIgGCgAAIBENFAAAQCJmIkdYdObw1157LZS76KKLQrnOnTuHcnvvvXcoJ0lTpkwJ5aZPnx7KjR8/PpTr2rVrKPfEE0+Ecs8++2wod8cdd4RyknT11VeHcr169Qqvsxhs3rxZy5cvz5mLziI/Y8aMUC7l59yuXbtQbsCAiq5u821PP/10KNesWbNQ7vbbbw/lorO5S9LYsWNDuejvdWlpaSi33377hXIoXhyBAgAASEQDBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhk7l5zGzOruY2h4K1YsSKU+9WvfhXK1atXL7zt9u3bh3L9+vUL5aL7csUVV4Ry69atC+U+//zzUO78888P5SSpY8eOoVx0tugBAwYscPce4QEUqJYtW/rxxx+fMzdmzJjQ+tavXx/KRX8XJOmYY44J5aIz4l933XWhXHQm8vvuuy+U27hxYygnSXPmzAnlRo4cGcpNnTo1lItezWC33XYL5VCwKqxfHIECAABIRAMFAACQiAYKAAAgEQ0UAABAIhooAACARDRQAAAAiWigAAAAEtFAAQAAJKKBAgAASFS/tgeA7Vfbtm1DudGjR4dyS5cuDW+7efPmodyVV14Zyk2fPj2Ue/HFF0O56D7ffffdoVz37t1DOUm6+eabQ7lOnTqF11kMGjVqpL333jtnrrS0NLS+6Kzhv/3tb0M5KT5r+axZs0K5xo0bh3KDBw8O5fbaa69QbsaMGaGcJPXp0yeUu+mmm0K5I444IpTbZ599Qrnf/OY3odyoUaNCORQOjkABAAAkooECAABIRAMFAACQiAYKAAAgEQ0UAABAIhooAACARDRQAAAAiWigAAAAEtFAAQAAJDJ3r7mNmdXcxrDd2WOPPcLZF154IZSLzjbcvn37UG7PPfcM5T788MNQLjqb9QMPPBDKSdLhhx8eykVns167du0Cd+8RHkCBat++vUdmi77//vtD67v66qtDuWnTpoVykjRz5sxQLjoj+DnnnBPKzZkzJ5SLvkZbt24dyknSkCFDQrm5c+eGcqeddloo165du1DusssuC+WuueaaUE6S3n333XAWVVZh/cp5BMrMJprZGjNbVGbZaDNbbmYLs/9OzOdoASBfqGEAqkPkLbx7JJ2wleW3uHu37L8n8jssAMibe0QNA5BnORsod39eUuz9BAAoMNQwANWhKieRX2Jmr2YPj7fI24gAoGZQwwBUWmUbqDsk/UBSN0krJd1cUdDMhpjZfDObX8ltAUC+hWpY2fr16aef1uDwABS6SjVQ7r7a3b9y968ljZfUcxvZce7eoxg+hQOgOERrWNn61aRJk5odJICCVqkGyszalPm2n6RFFWUBoNBQwwBUVf1cATObIuloSbuY2fuSRkk62sy6SXJJyyT9W/UNEQAqjxoGoDrkbKDcfcBWFt9VDWMBgLyjhgGoDlzKBQAAIFHOI1BAXfHee++Fs8OGDQvlli5dGsrtu+++odwZZ5wRyvXsWeHnMr5lv/32C+U6deoUyknSyJEjQ7kOHTqEcmvXrg1vu5A1btxYXbt2zZmbPXt2aH3R362LL744lJPil4eJXtpnwICtHbyr/Pqil0kZOHBgKCdJRx55ZCj3wQcfhHKnn356KNe4ceNQ7oADDgjlUp7nefPmhXLTp08PrxPpOAIFAACQiAYKAAAgEQ0UAABAIhooAACARDRQAAAAiWigAAAAEtFAAQAAJKKBAgAASEQDBQAAkMjcveY2ZlZzGwPyYMOGDaHciy++GMr17ds3lPv3f//3UG7u3Lmh3KZNm0I5STr77LNDueuvvz6U++CDDxa4e4/wAApUo0aNfPfdd8+Ze+6550LrW716dSh3113xy/ZFZ7NevHhxKFe/fuxiFccdd1wod8cdd4RyK1euDOUk6aSTTgrlxo4dG8rNmDEjlJs0aVIo9/LLL4dyP/nJT0I5SXrhhRdCuejvQ0lJSXjb26EK6xdHoAAAABLRQAEAACSigQIAAEhEAwUAAJCIBgoAACARDRQAAEAiGigAAIBENFAAAACJaKAAAAASMRM5kAfRWaXNLJTr1atXKHfkkUeGcg0aNAjlJGn27Nmh3FVXXRXKnXvuuUUxE/nuu+/uQ4cOzZn74osvQuvbf//9Q7nLL788lJOkd955J5SbOHFiKHfqqaeGck2aNAnlTjnllFDutttuC+VS1jly5MhQLnq1gE8++SSUq1evXigXvfqAJB1//PGhXPSqAtF92XXXXUO5IsNM5AAAAPlCAwUAAJCIBgoAACARDRQAAEAiGigAAIBENFAAAACJaKAAAAAS0UABAAAkooECAABIVL+2BwAUg9atW4dyzZo1C+U2bdoUynXu3DmUO+yww0I5Sfrzn/8cym3cuDG8zmKwefNmrV27NmcuOkP7QQcdFMqtWbMmlJOkUaNGhXK/+MUvQrlbb701lFu5cmUoF5nJXZLGjBkTyklSv379Qrl99tknlBswYEAoF93n0tLSUC76WpakG2+8MZQ76qijQrl99903lIvWr+jM5nUdR6AAAAAS0UABAAAkooECAABIRAMFAACQiAYKAAAgEQ0UAABAIhooAACARDRQAAAAiWigAAAAEjETOVCDojP07r777qFc7969Q7lnnnkmlJOkgw8+OJSbNGlSeJ3FoE2bNho+fHjOXHT27r59+4ZyP/vZz0I5SZoxY0Yo9/HHH4dy1157bSi3atWqUO7HP/5xKBedvVuSunXrFsr17NkzlJswYUIot3nz5lBu4sSJodzVV18dyknSBRdcEMr1798/lHv44YdDuTPOOCOUe/LJJ0O5999/P5QrVDmPQJnZHmb2rJktNrNSM/t5dnlLM5tpZm9mv7ao/uECQBz1C0B1ibyF96WkX7l7F0mHSRpqZl0kDZf0jLt3lPRM9nsAKCTULwDVImcD5e4r3f3v2dsbJC2R1E7SKZK2HMOfJOnUahojAFQK9QtAdUk6idzM2ks6SNI8Sa3dfcvlqFdJil2OHgBqAfULQD6FGygzayLpYUnD3P1bZ8K6u0vyCh43xMzmm9n8Ko0UACopH/Xrgw8+qIGRAqgrQg2UmTVQpvjc7+6PZBevNrM22fvbSFqztce6+zh37+HuPfIxYABIka/61apVq5oZMIA6IfIpPJN0l6Ql7v4fZe56TNLA7O2Bkh7N//AAoPKoXwCqS2QeqP8l6VxJr5nZwuyyEZJulPSgmV0g6Z+SzqyWEQJA5VG/AFSLnA2Uu78gySq4+7j8DgcA8of6BaC6MBM5UICiM/QeeuihodzXX38d3nb0ZOmU2aKLwbp160KzSj/6aOzdwJNPPjmUe/HFF0M5SRo9enQod8ghh4Ryt912Wyj3xhtvhHKNGjUK5aIzqktSixaxOVDXr18fXmfEkUceGcpFx1evXr3wttu1axfK/fKXvwzlBg8eHMr94x//COXmzJkTys2fH/tsWY8ehXkKNdfCAwAASEQDBQAAkIgGCgAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJDI3L3mNmZWcxsD8I3bb789nN20aVMot3Tp0lDuzjvvXODuhXkthgQHHnigP/XUUzlzHTp0CK1v+PDhody6detCOUlq3LhxKDdmzJhQbsiQIaFc9JIhO+wQ+5u9pKQklJOk3/3ud6Hcddddl9dtRy9Xcsstt4Ryd999dygnSd27dw/lopeRif4Mr7zyylDu+eefD+UWLlwYys2dOzeUk6SLLroonA2qsH5xBAoAACARDRQAAEAiGigAAIBENFAAAACJaKAAAAAS0UABAAAkooECAABIRAMFAACQiAYKAAAgUf3aHgCA6nfxxReHsw899FAod/bZZ4dyd955Z3jbhczMVK9evZy5vffeO7S+Cy+8MJT7/e9/H8pJ0hlnnBHKRWeRj85s3qxZs1Cuffv2oVx0BnRJql8/9t/YW2+9FcpFZ7KO7svo0aNDuSuuuCKUk6Tx48eHcuedd14ot3bt2lAuOnN4dKb7Y489NpQrLS0N5SSpdevWodxpp50WXmdFOAIFAACQiAYKAAAgEQ0UAABAIhooAACARDRQAAAAiWigAAAAEtFAAQAAJKKBAgAASEQDBQAAkMjcveY2ZlZzGwNQrWbNmhXKHXvssQvcvUc1D6famZlHZr1u165daH0lJSWh3KWXXhrKSVLbtm1DueiM0p07dw7l+vTpE8pNmzYtlFu/fn0oJ0nR/8OiM79v3rw5lJs8eXIoN3HixFBuwoQJoZwkvf/++6HcoEGDQrlXXnkllDv33HNDuSVLloRyu+22WyjXsGHDUE6Sdtxxx1AuOsP/0UcfXWH94ggUAABAIhooAACARDRQAAAAiWigAAAAEtFAAQAAJKKBAgAASEQDBQAAkIgGCgAAIBENFAAAQKKanol8raR/llu8i6R1NTaI6lUs+1Is+yGxL4VgT3fftbYHUVXUrzqFfSk8dXU/KqxfNdpAbXUAZvOL4TIPUvHsS7Hsh8S+oHoV03PCvhSmYtmXYtmPsngLDwAAIBENFAAAQKJCaKDG1fYA8qhY9qVY9kNiX1C9iuk5YV8KU7HsS7Hsxzdq/RwoAACAuqYQjkABAADUKbXaQJnZCWb2upktNbPhtTmWqjCzZWb2mpktNLP5tT2eFGY20czWmNmiMstamtlMM3sz+7VFbY4xqoJ9GW1my7PPzUIzO7E2xxhhZnuY2bNmttjMSs3s59nldfJ5KVbFUr8kalghKJb6JW0/NazWGigzqydprKQ+krpIGmBmXWprPHlwjLt3q4Mf07xH0gnllg2X9Iy7d5T0TPb7uuAefXdfJOmW7HPTzd2fqOExVcaXkn7l7l0kHSZpaPa1UVefl6JThPVLoobVtntUHPVL2k5qWG0egeopaam7v+3u/5I0VdIptTie7ZK7Py/pw3KLT5E0KXt7kqRTa3JMlVXBvtQ57r7S3f+evb1B0hJJ7VRHn5ciRf0qEMVSw4qlfknbTw2rzQaqnaT3ynz/fnZZXeSSnjKzBWY2pLYHkwet3X1l9vYqSa1rczB5cImZvZo9RF6nDhmbWXtJB0map+J7XuqyYqpfEjWskNXZ+iUVdw3jJPL8OMLduytzOH+omR1V2wPKF898TLMuf1TzDkk/kNRN0kpJN9fqaBKYWRNJD0sa5u6flL2vCJ4XFBZqWGGqs/VLKv4aVpsN1HJJe5T5fvfssjrH3Zdnv66RNF2Zw/t12WozayNJ2a9rank8lebuq939K3f/WtJ41ZHnxswaKFN47nf3R7KLi+Z5KQJFU78kalihqqv1S9o+alhtNlCvSOpoZnuZ2Y6SzpL0WC2Op1LMrMTMmm65Lam3pEXbflTBe0zSwOztgZIercWxVMmWF2tWP9WB58bMTNJdkpa4+3+UuatonpciUBT1S6KGFbK6WL+k7aeG1epEmtmPZN4qqZ6kie5+fa0NppLMrIMyf7FJUn1Jk+vSfpjZFElHK3Ol7NWSRkn6s6QHJX1fmavPn+nuBX9yYwX7crQyh79d0jJJ/1bmPfiCZGZHSJoj6TVJX2cXj1DmHII697wUq2KoXxI1rFAUS/2Stp8axkzkAAAAiTiJHAAAIBENFAAAQCIaKAAAgEQ0UAAAAIlooAAAABLRQAEAACSigQIAAEhEAwUAAJDo/wOYF4dBa8ZP4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.random.random((25,25))\n",
    "U = np.triu(A)\n",
    "L = np.tril(A)\n",
    "\n",
    "fig1, axes = plt.subplots(1, 2, figsize=(10,6))\n",
    "axes[0].imshow(U, cmap='gray')\n",
    "axes[0].set_title('$U$ is upper triangular')\n",
    "axes[1].imshow(L, cmap='gray')\n",
    "axes[1].set_title('$L$ is lower triangular')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8371307c",
   "metadata": {},
   "source": [
    "## 5.5 Matrix zoology (cont'd)\n",
    "\n",
    "### Orthogonal Matrix\n",
    "A matrix is orthogonal if it satisfies the following 2 criteria:\n",
    "- All columns are pairwise orthogonal eg dot product between any two columns is 0.\n",
    "- Each column has unit magnitude.\n",
    "\n",
    "A more compact expression for these criteria is:\n",
    "$$\n",
    "Q^T Q = I\n",
    "$$\n",
    "\n",
    "Notes\n",
    "- The standard mutually orthogonal unit vectors form an orthogonal matrix when concatenated along rows or columns.\n",
    "- The identity matrix $I$ is an orthogonal matrix.\n",
    "\n",
    "### Block Matrix\n",
    "A matrix composed of other matrices is referred to as a block matrix.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "A & B \\\\\n",
    "C & D \\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "abde9bd5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAEuCAYAAADVxIWbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAt60lEQVR4nO3de7xvVV3v/9fbzc1EBNyGCMglyfKKug9ScooUFTkGnoegeFLBMMoitTolZCFRJnp6aMfUlJQANVDQbGf4MxDJ/BkIGnIV2CIGOxTjpqhB4Of8McfCyeK7bvv7XWtN1no9H4/5WHOOMeYcY87vnGN9vvPynakqJEmSNCwPWe4GSJIk6YEM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0jSjJFck2W+52/Fg4fbSpkpySpI/GXMZxyf50KTaNClJHpvkziRrlrstQzfOZ5jkWUmubdv6RRNu2pKaxL6c5Igkn59Um2ao4/FJLkny3SSvXYw6DNJWoSSbJbkqybWzlauqJ1bV+Qtc9pZJPpDkG23HvSTJC2You3WS65P8Ui/t4Un+Lckh4yy7lb8+yd1J1k5L/9cklWS3ea7T9Un2n6vcpmwvrQ5tH/pB+wd6W5J/SLLLMrbnlHYMHDwt/R0t/Yh5LmfOY6Oq/q2qtq6qezexrc9P8rl2zH87yT8lOWhTlrXCnQC8q23rT/QzNqHv/M0klyfZopf2+tZ3bjai/D5Jzklya/uMzkyy4yzLPz/Jf7bj4Y72+T5501Z7Wf0e8NmqenhVvXOmQu14u2e2bTITg7TV6deAHwf2SPKwCS97M+AG4OeBRwB/AHx0VEBUVXcCvwr8eZJHteS3ARdX1VnjLLvn68DLpiZaR/BjC1ul2Y3qtKQRfrGqtgZ2BL4F/MUyt+ca4JVTE20/fgnwtUlVMO6x0b6snQmcBuwM7AAcB/zi+K1bcXYFrpghb6F957uB24E3AiTZA/gj4MiqumdE+e2Ak4DdWju+C/z1HO09uh0P2wPnAx+co/wQzbbNAWj/Y18M3AG8fME1VJXDKhroDtBv0wUudwPPnKXs9cD+bfwNwEa6g+9q4DkLqPNS4MWz5J8CnA7sB9wCPHoSy27t/wPgol7an9F1PAXs1ks/hu6f03eBK4H/2dI/CPwQ+AFwJ/B7vWW/odV/F10neD2wP/ATwK3A01vZx7Rtvt9yf/4OyzP0j6U2fSBwTW/6FOBPetO/Amxo+9F64DG9vCcC57S8bwG/39KPBz7Uxjdvx9THgC1GtOeUdix8C9iupb0Q+BTweeCIlvYTwHntuPwP4MPAti3vAccG3T/pAo4E/g34XC9tM7p/yDfSBawAW7f1fOWINqYt43dn2a4Pacf4N4Cb6YK5R7S8qXpfRReg3Eb3BfW/teP2drozT1PLOgL4/4F30f1D/Sq9fq4dx+vbdt8A/Eovb2/gYuA7bZu+vZe3D/CFVt9X+v0AsDvwT3T9zjmt7g/Nsr4j9wu6vqv/WWw5j31yrn758W19ngx8BnjLAvb3pwPfnSX/fODVveknAHf3po/vbwfgILpg6PY270/38nYBPk7Xx94y9Zm2z/PzvXL/h27ffgTwuLbd76Dbrz8yS1tH1k13XNwL/Gfb5j85w/yvbPvf64DLF9x3jNv5ODy4BrqO+cI2fnm/oxlR9nq6oOPxbSeb6hB2A35invXt0Hbin5qlzHbATe1gedUC1mXWZffafzXw08Aaun8Qu/LAIO1Quk74IcBLge8BO/aXM2LZl7QO4qHTy9F1plfSnbX7NPBny/3ZOyzfMG3f+DHgVOC0Xv4ptCANeHY7Fp4ObEl3xu1zLe/h7Vj5HWCrNv3Mlnc88CHgocA/tGWumaE9pwB/Qnf24zUt7aN0X976QdrjgOe2djyKLuj681Hr1aZ3a8fWacDDWlum0jZrZZ4HfJPubP5fAWfN0MafavPtPst2/WW6oGUPuoDv48AHp7XlvW1bPY+uv/hEq3snusDu51v5I4B7gN+iC3JfSvdPfPuW/zngPW1Ze9EFBc9uef8CvKKNbw3s08Z3ogscDqTrW57bph/Vm+/tbfv+HF2wNjJIm22/GPVZzLE/ztkvt3LHtjqvBrZawP7+euCCWfLPpwVpwBbAm6ety/H86AvHT9L1x89tn8vvtc98C7o+/SvAO9r+thWwb+/z/Hzb7n9F1w//WMs7ne7L+kP684xo54x1T1+PWdb1M3RXiHZo+9czFtR3LFan5DC8ga4ju4v2TQ74CPAXs5S/ni7IeRxdZ7Y/sPkC6tscOBd43zzKngt8n/YteBLL7rX/D4C3AAfQfVvdjGlB2oh5LwEO7i9nxLJ/eVR9ven1wGV031jn/GbrsHKHtm/cSfdt/L+Afwee3Ms/hR8FaR8A3tbL27rNsxtdEPWvM9RxfNvn/gl4J5BZ2nMKXZC2L12gsC3dGaCH0gvSRsz3on79I/b53dqxtceItM16aX/Rjo2NwCNnqOtZbb4ZgwO6f4C/3pt+fNtWm/Xq3amXfwvw0t70x4DXt/Ej2ueSXv4XgVfQfRm7F3h4L+8twClt/HN0lwPXTmvfG2hBYy/t08DhwGPp/mk/rJf3N8wcpM24X4z6LGbZZgvpl/dt2/DNC9jXn0J3pu+/z1LmfLr+/na6/0l3cP+zlsfzoyDtD4GP9vIe0vab/YCfoQuWNxtRxxHAhXT/5+53RpnuS8RJwM5zrMuMdffWY8YgrX3GPwT26n32/3e+27KqvCdtlXkbcG796Ob2K4CnzjVTVW2g+2Z0PHBzkjOSPGa2eZI8hO5yyN3A0XOUfTldh3ou8Na52rOQZTcfBP4X3UF72gzLfGW7mfb2JLcDTwLWjirbc8Mc+X/VlvMXVXXXPNqple1FVbUt3Tf3o4F/SvLoEeUeQ3f5Drjv3s1b6M7K7MLs94ztQ/dP8sRq/xVmU1WfpztD9kbgk1X1g35+kh3a8b4xyXfoztTNdVzA3MfGSXTHxilVdcsMZabSZ7vZ+n7bqo1vRnfWYsq3euM/GDG9dW9647Tt9o1Wx2OAW6vqu9PydmrjR9KddflqkouSvLCl7wocOtWvtL5l37ZOjwFuq6rvTVvmvNZ12n4xLwvsl7cA3kcXUB/d7kuba/mPo7tk/rqq+uc5ir+2HQ8PpbvUflaSp4woN329f0i3f00dD9+o0ffJQXeC4WDgj6rq7l7679FdTv9iuqfyf3mG+Werez5eAVxVVZe06Q8D/yvJ5vOc3yBttUiyL93Ni/sm+WaSb9JdMpnXEzVV9TdVtS8/ulQ4YzCVJHTf+nagu+fhv2Yp++N0p6p/he4hgpck+e+TWHav7d+ge4DgQLrLIdOXuStdQHU03bf6bekuBWdqETMtepZ2bg38eWvr8Um2n6udWh2q6t6q+jjdmZl9RxT5d7rjDLjvxuNH0n2Dv4HujPhM/pHuDM9nkuwwS7m+D9H1BaO+wPwp3X7+5Krahu7G5/TyN+XYWEMXpJ0G/Hr7xz7K1XTr++JZ2n6/bcWPzk59a3TxOe3U+pj+8v69Ddsnefi0vI0AVXVtVb2M7jLqW+kCjoe19n+wqrbtDQ+rqhPpLltvN+3hrcfO0rbZ9os5bULf+Yd0V1BeR3fJ+H1zLH9Xui/af1xV834IoKp+2AK6DXSXpKebvt6hC86mjofHzvKAylV09yR+Ksnje3V+s6p+paoeQ/d/5z0z7Iez1T0fr6R7QG/qf+7b6b7kHDjP+Q3SVoO2Y72d7kB7PN39FHvRHRDbJpmtY5j6LZhnJ9mS7j6GH9Cdwp3JX9LdA/aL07+Zj/Au4BNV9dmquonuG85ftbrGXXbfkXT3j3xvRN7D6P6pfBsgyavovuVP+Raz/2Mc5f/SPaX6arr7g967wPm1QqVzMN29mFeNKHI68Koke7Xj4E/p7iO9HvgksGP7OYQt0/1kzTP7M1fV2+gum30m035+ZgbvpLvn5nMj8h5Od5n2jiQ7Ab87LX9Tjo3fpzvefpnuZu7TMuI31NoZrd8G/jDJq5Jsk+QhSfZNclIrdjrwW0l2b1+M/pTuJvCZzqzM5ceB1ybZPMmhdH3N2VV1A93N/29JslU743MkXYBLkpcneVQ703J7W9YPW/4vpvsZkTVt3v2S7Ny+PF4M/FGSLdoX6dmeWp1tv5iPefedSZ4KvJbunuWiu4qyW+sbR5Xfie5G+ndV1YL7uiQ/Q/fwwKgnJT8K/I8kz2lnoH6H7hLpF+guR98EnJjkYW37Pqs/c1WdTrfPnZvkJ1p9hybZuRW5jW5/HPU/bba657NOP0H3UMlebXgS3bH5yhlnnG4h10YdHpwD3bffbwBbT0sP3dM7L5xhvuvp7ul6Ct3B8F26ew0+Se9ps2nzTJ1pm3riZWr4pRFlX0T3TWXbaennMeIeiIUsu9/+EekPuCeN7sbVW+lukn073X09Uze2Hkz3lNntwP+eadm97XUw3TetqRuOp55gG9lOh5U/tH1j6sm779Kdqf2lXv4p3P/pzl+ju6w5dbzt3Mt7Et29WLfR3YB/TEs/nvs/EfcndPdWbj+iPferb1pe/8GBJwJfau2+hO6f1I29svc7Nhh9/9l9acAzWrsf1/LW0D1R+cZZtt0BwD+3Nnyb7j6g/9HyHkL3kxw3tLwP8aOnVUe15Ubu/3Tlh4A/aONHcP+nO68Bntcru3P7LG5tn82vTVvOza2NV9Bd2p7KeyZdf3Jra+M/AI9teXv01m0+T3fOtl9czwz3pLGwfnkNXfD4e9PS96PrH3cYMc+b2vL7y75zlvU4f1pbNgC/1cs/nvvvy/+T7kGsO9q2fGIv77F0D4NMPYH8zt7n2X+681fo/g/uRnfrz8ZW99eAo2Zp62x1n88M96TRfTH/2Ij0vekCvQccl6OGtJkkSVq10v2A76uru61DGgQvd0qSJA3QWEFaku3TvQri2vZ3uxnK3ZvuyblLkqzvpe+e5MIkG5J8JL1XUEjSchm3b5OkSRjrcmeSt9E9lnxikmPo7gV4w4hyd1b3+ofp6R8FPl5VZyR5L/CVqvrLTW6QJE3AuH2bJE3CuEHa1XQ3Yd6U7sWh51fV40eUe0BH1p44/DbdK4DuaU9CHF9Vz9/kBknSBIzTt0nSpIx7T9oO1f1sAnRPGc30uzxbJbk4yQVJXtTSHgncXj96VPpGFvCjfJK0iMbp2yRpImb6Abj7JDkXGPWr2G/sT1RVJZnptNyuVbUx3S8Wn5fkMrrHWectyVHAUW3yGQuZVyvHdtuNvDVIK9z3vvc97rrrrsxdcv4Wq2+rqge8EcD+q7Paj9/bbrttuZuwbPzsb/uPqnrUQuebM0irqv1nykvyrSQ79i4J3DzDMqZ+lfm6JOcDT6N7l9a2STZrZ9N2ZpZf8a2qk+h+pZpZOkytcPvvP+PuqBXs3HPPnfgyF7Fve0CQZv/VWe3H75lnnrncTVg2fvZnzvbKrxmNe7lzPd2LYml//256gSTbtV9Hpv369bOAK6u7Ge6zwCGzzS9Jy2CT+7Yla6GkFW/cIO1E4LlJrqX7pfUTAZKsS/L+VuangYuTfIUuKDuxqqY6sjcAv51kA909ah8Ysz2SNAnj9m2SNLY5L3fOpqpuAZ4zIv1i4NVt/AvM8BLvqrqO7hUJkjQY4/ZtkjQJvnFAkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmAxgrSkmyf5Jwk17a/240os1eSf0lyRZJLk7y0l3dKkq8nuaQNe43THkmapCQHJLk6yYYkx4zI3zLJR1r+hUl2W4ZmSlqhxj2TdgzwmaraE/hMm57u+8Arq+qJwAHAnyfZtpf/u1W1VxsuGbM9kjQRSdYA7wZeADwBeFmSJ0wrdiRwW1U9DngH8NalbaWklWzcIO1g4NQ2firwoukFquqaqrq2jf87cDPwqDHrlaTFtjewoaquq6q7gTPo+ry+fh94FvCcJFnCNkpawcYN0naoqpva+DeBHWYrnGRvYAvga73kN7fLoO9IsuWY7ZGkSdkJuKE3fWNLG1mmqu4B7gAeuSStk7TibTZXgSTnAo8ekfXG/kRVVZKaZTk7Ah8EDq+qH7bkY+mCuy2Ak4A3ACfMMP9RwFFztVeShsb+S9KmmDNIq6r9Z8pL8q0kO1bVTS0Iu3mGctsA/wC8saou6C176izcXUn+Gvjfs7TjJLpAjtmCQUmakI3ALr3pnVvaqDI3JtkMeARwy/QF2X9J2hTjXu5cDxzexg8H/m56gSRbAH8LnFZVZ03L27H9Dd39bJeP2R5JmpSLgD2T7N76scPo+ry+fh94CHBeVRmESZqIcYO0E4HnJrkW2L9Nk2Rdkve3Mi8Bfg44YsRPbXw4yWXAZcBa4E/GbI8kTUS7x+xo4NPAVcBHq+qKJCckOagV+wDwyCQbgN9m9BPukrRJ5rzcOZuqugV4zoj0i4FXt/EPAR+aYf5nj1O/JC2mqjobOHta2nG98f8EDl3qdklaHXzjgCRJ0gAZpEmSJA2QQZokSdIAGaRJkiQNkEGaJEnSABmkSZIkDZBBmiRJ0gAZpEmSJA2QQZokSdIAGaRJkiQNkEGaJEnSABmkSZIkDZBBmiRJ0gAZpEmSJA2QQZokSdIAGaRJkiQNkEGaJEnSABmkSZIkDZBBmiRJ0gAZpEmSJA2QQZokSdIAGaRJkiQNkEGaJEnSAE0kSEtyQJKrk2xIcsyI/C2TfKTlX5hkt17esS396iTPn0R7JGkS5tG3HZHk20kuacOrl6OdklamzcZdQJI1wLuB5wI3AhclWV9VV/aKHQncVlWPS3IY8FbgpUmeABwGPBF4DHBukp+sqnvHbZckjWOefRvAR6rq6CVvoKQVbxJn0vYGNlTVdVV1N3AGcPC0MgcDp7bxs4DnJElLP6Oq7qqqrwMb2vIkabnNp2+TpEUz9pk0YCfght70jcAzZypTVfckuQN4ZEu/YNq8O02gTZI0rvn0bQAvTvJzwDXAb1XVDSPK3GePPfbgxBNPnFwrH0QOPfTQ5W7CsjrzzDOXuwnLZrV/9t15qYV70Dw4kOSoJBcnuXi52yJJzd8Du1XVU4Bz+NEVg/vp91/f+c53lrSBkh68JhGkbQR26U3v3NJGlkmyGfAI4JZ5zgtAVZ1UVeuqat0E2ixJc5mzf6qqW6rqrjb5fuAZoxbU77+22WabRWmspJVnEkHaRcCeSXZPsgXdgwDrp5VZDxzexg8BzquqaumHtac/dwf2BL44gTZJ0rjm7NuS7NibPAi4agnbJ2mFG/uetHaP2dHAp4E1wMlVdUWSE4CLq2o98AHgg0k2ALfSdXa0ch8FrgTuAX7DJzslDcE8+7bXJjmIrv+6FThi2RosacWZxIMDVNXZwNnT0o7rjf8nMPKuwap6M/DmSbRDkiZpHn3bscCxS90uSavDg+bBAUmSpNXEIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgCYSpCU5IMnVSTYkOWZE/m8nuTLJpUk+k2TXXt69SS5pw/pJtEeSxpXk5CQ3J7l8hvwkeWfr9y5N8vSlbqOklW3sIC3JGuDdwAuAJwAvS/KEacX+FVhXVU8BzgLe1sv7QVXt1YaDxm2PJE3IKcABs+S/ANizDUcBf7kEbZK0ikziTNrewIaquq6q7gbOAA7uF6iqz1bV99vkBcDOE6hXkhZNVX0OuHWWIgcDp1XnAmDbJDsuTeskrQaTCNJ2Am7oTd/Y0mZyJPCp3vRWSS5OckGSF02gPZK0FBba90nSgmy2lJUleTmwDvj5XvKuVbUxyR7AeUkuq6qvjZj3KLpLCpL0oNLvv9auXbvMrZH0YDGJM2kbgV160zu3tPtJsj/wRuCgqrprKr2qNra/1wHnA08bVUlVnVRV66pq3QTaLEnjmlffB/fvv7bZZpslaZykB79JBGkXAXsm2T3JFsBhwP2e0kzyNOB9dAHazb307ZJs2cbXAs8CrpxAmyRpsa0HXtme8twHuKOqblruRklaOca+3FlV9yQ5Gvg0sAY4uaquSHICcHFVrQf+D7A1cGYSgH9rT3L+NPC+JD+kCxhPrCqDNEnLLsnpwH7A2iQ3Am8CNgeoqvcCZwMHAhuA7wOvWp6WSlqpJnJPWlWdTddh9dOO643vP8N8XwCePIk2SNIkVdXL5sgv4DeWqDmSViHfOCBJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEATCdKSHJDk6iQbkhwzIv+IJN9OckkbXt3LOzzJtW04fBLtkaRxJTk5yc1JLp8hf78kd/T6teOWuo2SVrbNxl1AkjXAu4HnAjcCFyVZX1VXTiv6kao6etq82wNvAtYBBXypzXvbuO2SpDGdArwLOG2WMv9cVS9cmuZIWm0mcSZtb2BDVV1XVXcDZwAHz3Pe5wPnVNWtLTA7BzhgAm2SpLFU1eeAW5e7HZJWr0kEaTsBN/Smb2xp0704yaVJzkqyywLnlaQh+pkkX0nyqSRPXO7GSFpZxr7cOU9/D5xeVXcl+VXgVODZC1lAkqOAowDWrl3Le97znsm3UoN36KGHLncTtAzWrVu33E0Y5cvArlV1Z5IDgU8Ae44q2O+/AF7ykpcsSQOHZrUfv2eeeeZyN2HZrPbPflNN4kzaRmCX3vTOLe0+VXVLVd3VJt8PPGO+8/aWcVJVrauqddtss80Emi1Jm66qvlNVd7bxs4HNk6ydoex9/deSNlLSg9okgrSLgD2T7J5kC+AwYH2/QJIde5MHAVe18U8Dz0uyXZLtgOe1NEkatCSPTpI2vjddf3rL8rZK0koy9uXOqronydF0wdUa4OSquiLJCcDFVbUeeG2Sg4B76G7EPaLNe2uSP6YL9ABOqCpv1JW07JKcDuwHrE1yI92T6JsDVNV7gUOA1yS5B/gBcFhV1TI1V9IKNJF70tqp/rOnpR3XGz8WOHaGeU8GTp5EOyRpUqrqZXPkv4vuJzokaVH4xgFJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAEySJMkSRoggzRJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAEySJMkSRoggzRJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAEySJMkSRoggzRJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAGaSJCW5IAkVyfZkOSYEfnvSHJJG65Jcnsv795e3vpJtEeSxpVklySfTXJlkiuSvG5EmSR5Z+v7Lk3y9OVoq6SVabNxF5BkDfBu4LnAjcBFSdZX1ZVTZarqt3rlfxN4Wm8RP6iqvcZthyRN2D3A71TVl5M8HPhSknP6fRvwAmDPNjwT+Mv2V5LGNokzaXsDG6rquqq6GzgDOHiW8i8DTp9AvZK0aKrqpqr6chv/LnAVsNO0YgcDp1XnAmDbJDsucVMlrVCTCNJ2Am7oTd/IAzsyAJLsCuwOnNdL3irJxUkuSPKiCbRHkiYqyW50VwAunJY17/5PkhZq7MudC3QYcFZV3dtL27WqNibZAzgvyWVV9bXpMyY5CjgKYO3atUvTWkmrXpKtgY8Br6+q72ziMu7rvyRpviZxJm0jsEtveueWNsphTLvUWVUb29/rgPO5//1q/XInVdW6qlq3zTbbjNtmSZpTks3pArQPV9XHRxSZV//X778Wp6WSVqJJBGkXAXsm2T3JFnSB2AOe0kzyU8B2wL/00rZLsmUbXws8C7hy+ryStNSSBPgAcFVVvX2GYuuBV7anPPcB7qiqm5askZJWtLEvd1bVPUmOBj4NrAFOrqorkpwAXFxVUwHbYcAZVVW92X8aeF+SH9IFjCdOe3JKkpbLs4BXAJcluaSl/T7wWICqei9wNnAgsAH4PvCqpW+mpJVqIvekVdXZdJ1VP+24adPHj5jvC8CTJ9EGSZqkqvo8kDnKFPAbS9MiSauNbxyQJEkaIIM0SZKkATJIkyRJGiCDNEmSpAEySJMkSRoggzRJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAEySJMkSRoggzRJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAEySJMkSRoggzRJkqQBMkiTJEkaIIM0SZKkATJIkyRJGiCDNEmSpAGaSJCW5OQkNye5fIb8JHlnkg1JLk3y9F7e4UmubcPhk2iPJI0ryS5JPpvkyiRXJHndiDL7JbkjySVtOG452ippZdpsQss5BXgXcNoM+S8A9mzDM4G/BJ6ZZHvgTcA6oIAvJVlfVbdNqF2StKnuAX6nqr6c5OF0/dM5VXXltHL/XFUvXIb2SVrhJnImrao+B9w6S5GDgdOqcwGwbZIdgecD51TVrS0wOwc4YBJtkqRxVNVNVfXlNv5d4Cpgp+VtlaTVZKnuSdsJuKE3fWNLmyldkgYjyW7A04ALR2T/TJKvJPlUkicubcskrWSTuty56JIcBRwFsHbt2mVujaTVIsnWwMeA11fVd6ZlfxnYtaruTHIg8Am62zqmL+N+/dd73vOexW30QB166KHL3YRldeaZZy53E5bNav/sk2zSfEt1Jm0jsEtveueWNlP6A1TVSVW1rqrWbbPNNovWUEmakmRzugDtw1X18en5VfWdqrqzjZ8NbJ7kAd8i7b8kbYqlCtLWA69sT3nuA9xRVTcBnwael2S7JNsBz2tpkrSs0n31/QBwVVW9fYYyj27lSLI3XZ96y9K1UtJKNpHLnUlOB/YD1ia5ke6Jzc0Bquq9wNnAgcAG4PvAq1rerUn+GLioLeqEqprtAQRJWirPAl4BXJbkkpb2+8Bj4b6+7RDgNUnuAX4AHFZVtQxtlbQCTSRIq6qXzZFfwG/MkHcycPIk2iFJk1JVnwdmvZGkqt5F9/NDkjRxvnFAkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmADNIkSZIGyCBNkiRpgAzSJEmSBsggTZIkaYAM0iRJkgbIIE2SJGmAJhKkJTk5yc1JLp8h/5eSXJrksiRfSPLUXt71Lf2SJBdPoj2SNK4kWyX5YpKvJLkiyR+NKLNlko8k2ZDkwiS7LUNTJa1QkzqTdgpwwCz5Xwd+vqqeDPwxcNK0/F+oqr2qat2E2iNJ47oLeHZVPRXYCzggyT7TyhwJ3FZVjwPeAbx1aZsoaSWbSJBWVZ8Dbp0l/wtVdVubvADYeRL1StJiqc6dbXLzNtS0YgcDp7bxs4DnJMkSNVHSCrcc96QdCXyqN13APyb5UpKjlqE9kjRSkjVJLgFuBs6pqgunFdkJuAGgqu4B7gAeuaSNlLRibbaUlSX5Bbogbd9e8r5VtTHJjwPnJPlqOzM3fd6jgKMA1q5duyTtlbS6VdW9wF5JtgX+NsmTqmrkvbezsf+StCmW7ExakqcA7wcOrqpbptKramP7ezPwt8Deo+avqpOqal1Vrdtmm22WosmSBEBV3Q58lgfee7sR2AUgyWbAI4BbppWx/5K0SZYkSEvyWODjwCuq6ppe+sOSPHxqHHgesOBvqZI0aUke1c6gkeShwHOBr04rth44vI0fApxXVdPvW5OkTTKRy51JTgf2A9YmuRF4E91NtlTVe4Hj6O7TeE+7p/ae9iTnDnSXEKba8jdV9f9Nok2SNKYdgVOTrKH7QvvRqvpkkhOAi6tqPfAB4INJNtA9PHXY8jVX0kozkSCtql42R/6rgVePSL8OeOoD55Ck5VVVlwJPG5F+XG/8P4FDl7JdklYP3zggSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRABmmSJEkDZJAmSZI0QAZpkiRJA2SQJkmSNEAGaZIkSQNkkCZJkjRAEwnSkpyc5OYkl8+Qv1+SO5Jc0objenkHJLk6yYYkx0yiPZI0riRbJflikq8kuSLJH40oc0SSb/f6tlcvR1slrUybTWg5pwDvAk6bpcw/V9UL+wlJ1gDvBp4L3AhclGR9VV05oXZJ0qa6C3h2Vd2ZZHPg80k+VVUXTCv3kao6ehnaJ2mFm8iZtKr6HHDrJsy6N7Chqq6rqruBM4CDJ9EmSRpHde5sk5u3oZaxSZJWmaW8J+1n2mWDTyV5YkvbCbihV+bGliZJyy7JmiSXADcD51TVhSOKvTjJpUnOSrLL0rZQ0ko2qcudc/kysGu7bHAg8Algz4UsIMlRwFFt8q6XvOQlI+9/WyJrgf9YhXWv9vpX87ovd/2PX45Kq+peYK8k2wJ/m+RJVdXve/4eOL2q7kryq8CpwLOnL8f+azD1r+Z1X+76V/O6wyb2YamazNn7JLsBn6yqJ82j7PXAOrpA7fiqen5LPxagqt4yx/wXV9W6cdu8qZaz/tW87std/2pe9+Wuf7nXvbXhOOD7VfVnM+SvAW6tqkfMsZxV+zkud/2red2Xu/7VvO7j1L8klzuTPDpJ2vjerd5bgIuAPZPsnmQL4DBg/VK0SZJmk+RR7QwaSR5K94DTV6eV2bE3eRBw1ZI1UNKKN5HLnUlOB/YD1ia5EXgT3U22VNV7gUOA1yS5B/gBcFh1p/DuSXI08GlgDXByVV0xiTZJ0ph2BE5tZ8geAny0qj6Z5ATg4qpaD7w2yUHAPXQPTx2xbK2VtOJMJEirqpfNkf8uup/oGJV3NnD2Aqs8aYHlJ20561/N677c9a/mdV/u+pe87qq6FHjaiPTjeuPHAscucNGr+XNc7vpX87ovd/2red03uf6J3ZMmSZKkyfG1UJIkSQM02CAtyfZJzklybfu73Qzl7u29kmV9L333JBe21019pD2YMLG6k+yV5F/a62IuTfLSXt4pSb7ea9de86x31ldkJdmyrcuGtm679fKObelXJ3n+fNd1AXX/dpIr27p+JsmuvbyRn8GE65/x9TtJDm+f1bVJDl+k+t/Rq/uaJLf38sZa/8z9WrUkeWdr26VJnt7Lm8S6z1X/L7V6L0vyhSRP7eVd39IvSXLxItT9oH+l3Lh92Rj1bnJ/sgR1L9rrtMY5npao/hn36QnVv0uSz7b++ookrxtRZlG2wTzrXrT1z/xe5baY+/3kXyVXVYMcgLcBx7TxY4C3zlDuzhnSP0r3gALAe4HXTLJu4CeBPdv4Y4CbgG3b9CnAIQtc3zXA14A9gC2ArwBPmFbm14H3tvHD6F5HA/CEVn5LYPe2nDUTrvsXgB9r46+Zqnu2z2DC9R8BvGvEvNsD17W/27Xx7SZd/7Tyv0n3kMuk1v/ngKcDl8+QfyDwKSDAPsCFk1r3edb/s1PLBV4wVX+bvh5Yu4jrvh/dT/uM9Zkt5zCf/mQS+9FCt89M/ckS1T3yeJ5Q/Zt0PC1h/SP36QnWvyPw9Db+cOCaEdt/UbbBPOtetPVv67N1G98cuBDYZ1qZRdnvF1D/gvb9wZ5Jo3s91Klt/FTgRfOdMUnoflDyrE2Zfz51V9U1VXVtG/93ul8kf9QC6phuPq/I6rfrLOA5bV0PBs6oqruq6uvAhra8idVdVZ+tqu+3yQuAnRew/LHrn8Xz6X4J/taqug04Bzhgket/GXD6AuuYUc39WrWDgdOqcwGwbbqffpjEus9Zf1V9oS0fJvzZz2PdZ/JgeqXcJvdlYxinP1mKuhfNGMfTUtW/qKrqpqr6chv/Lt3Pwkx/k8+ibIN51r1o2vrM9Sq3xdrv51v/ggw5SNuhqm5q498Edpih3FZJLk5yQZIXtbRHArdX1T1teqGvm5pv3cB9v/22Bd23xylvbqeR35Fky3nUOZ9XZN1Xpq3bHXTrOu7rtRY6/5F038KmjPoMFmK+9Y96/c4kXi0272Wku8y7O3BeL3nc9d/U9i3Ha9Wmf/YF/GOSL6X7Vf3F8GB/pdw4fdmmGqc/Gdc4x/NSGMK+M2qfnrh2Ke9pdGd0+hZ9G8xSNyzi+mfuV7kt1n4/3/phAfv+Ur0WaqQk5wKPHpH1xv5EVVWSmaLRXatqY5I9gPOSXEa30Zei7qkfs/wgcHhV/bAlH0vXGW9B99jtG4AT5mrTg0GSl9O9LeLne8kP+Ayq6mujl7DJ5vX6nSVwGHBWda8LmrIU67/skvwCXZC2by9537buPw6ck+Sr7UzCpIz9SrmlsFh92Urcj5qhHM/LYUn26SRbAx8DXl9V35n08seoe1HXv+Z+lduimkf9C9r3l/VMWlXtX1VPGjH8HfCtqdOv7e/NMyxjY/t7HXA+XeR+C93p26kgdGdg46TrTrIN8A/AG9sp46ll39ROe94F/DXzu/S4EehH1A9oc79MW7dHtHWdz7zj1k2S/en+6RzU1g2Y8TNYiDnrr6pbenW+H3jGQto+bv09hzHtUucE1n9T2zeJdZ+XJE+h2+4HV9UtU+m9db8Z+FsWdpl9TlX1nanLB9X9puLmSdayhOs+H4vYl22qcfqTcY1zPC+FZd13ZtmnJybJ5nRB0oer6uMjiizaNpir7qVY/7bs24HP8sBbQBZrv59X/Qvd94d8uXM9MPW02uHA300vkGS7qUuJ7UN+FnBlVRXdxjlktvnHrHsLun9Kp1XVWdPypjrk0N1/Mp8ofj6vyOq36xDgvLau64HD0j21sjvdt5IvzqPOeded5GnA++gCtJt76SM/gwXUPd/6Z3r9zqeB57V2bAc8r6VNtP7Whp+iu0H/X3ppk1j/uawHXpnOPsAd7fLZJNZ9TkkeC3wceEVVXdNLf1iSh0+Nt/on+o01K+OVcpvcl41R5zj9ybjGOZ6XwkzH05KYZZ+e1PIDfAC4qqrePkOxRdkG86l7Mdc/83iVG4u338+r/gXv+7VIT5iMO9BdI/4McC1wLrB9S18HvL+N/yxwGd3TQ5cBR/bm34MuUNkAnAlsOeG6Xw78F3BJb9ir5Z3X2nM58CHa0x7zqPdAuqdhvkZ3dg66y6QHtfGt2rpsaOu2R2/eN7b5rgZesAnbe666zwW+1VvX9XN9BhOu/y3AFa2ezwI/1Zv3l9s22QC8ajHqb9PHAydOm2/s9ac7M3dT259upLuk+GvAr7X8AO9ubbsMWDfhdZ+r/vcDt/U++4t7x9hX2nDF1HabcN1H9z73C4Cfne0zG+LAmH3ZGPVucn+yBHXPeDxPoO5NPp6WqP4Z9+kJ1b8v3b2il/aO2QOXYhvMs+5FW3/gKcC/tvovB45b4v1+PvUvaN/3jQOSJEkDNOTLnZIkSauWQZokSdIAGaRJkiQNkEGaJEnSABmkSZIkDZBBmiRJ0gAZpEmSJA2QQZokSdIA/T+lXYQHkZSIbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([1,2,3,4]).reshape(2,2)\n",
    "block = np.block([[A,A],[A,A]])\n",
    "\n",
    "fig1, axes = plt.subplots(1, 2, figsize=(10,6))\n",
    "axes[0].imshow(A, cmap='gray')\n",
    "axes[0].set_title('$A$ is 2 X 2 Matrix')\n",
    "axes[1].imshow(block, cmap='gray')\n",
    "axes[1].set_title('Block Matrix Composed of 2 X 2 Blocks of A')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81f3c372",
   "metadata": {},
   "source": [
    "## 5.5 Matrix zoology (cont'd)\n",
    "\n",
    "### Toeplitz Matrix\n",
    "A matrix constructed from row and column vectors in which the descending diagonals are constant value.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & b & c & d \\\\\n",
    "b & a & b & c \\\\\n",
    "c & b & a & b \\\\\n",
    "d & c & b & a \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Notes\n",
    "* Toeplitz matrices are not always symmetric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1f8bd32b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3 4]\n",
      " [2 1 2 3]\n",
      " [3 2 1 2]\n",
      " [4 3 2 1]]\n",
      "[[1 4 5 6]\n",
      " [2 1 4 5]\n",
      " [3 2 1 4]]\n"
     ]
    }
   ],
   "source": [
    "def toeplitz(r, c):\n",
    "    \"\"\"\n",
    "    toeplitz returns the toeplitz matrix constructed from row and column vectors\n",
    "    \n",
    "    :param r: numpy.ndarray  Row vector.\n",
    "    :param c: numpy.ndarray  Column vector.\n",
    "    :return: numpy.ndarray   Toeplitz matrix.\n",
    "    \"\"\"\n",
    "    assert r.ndim == c.ndim == 1\n",
    "    assert r.dtype == c.dtype\n",
    "\n",
    "    nrow, ncol = c.shape[0], r.shape[0]\n",
    "    A = np.empty((nrow,ncol), dtype=r.dtype)\n",
    "    for i in range(nrow):\n",
    "        A[i,i:] = r[:ncol-i]\n",
    "        A[i:,i] = c[:nrow-i]\n",
    "    \n",
    "    return A\n",
    "\n",
    "\n",
    "a = np.array([1,2,3,4], dtype=int)\n",
    "T = toeplitz(a, a)\n",
    "print(T)\n",
    "expected = np.array([[1,2,3,4],[2,1,2,3],[3,2,1,2],[4,3,2,1]], dtype=int)\n",
    "np.testing.assert_equal(T, expected, \"toeplitz, symmetric\")\n",
    "\n",
    "a = np.array([1,4,5,6], dtype=int)\n",
    "b = np.array([1,2,3], dtype=int)\n",
    "T = toeplitz(a, b)\n",
    "print(T)\n",
    "expected = np.array([[1,4,5,6],[2,1,4,5],[3,2,1,4]], dtype=int)\n",
    "np.testing.assert_equal(T, expected, \"toeplitz, non-symmetric\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "017da782",
   "metadata": {},
   "source": [
    "## 5.8 Shifting a Matrix\n",
    "Shifted matrix $\\hat{A}$ is the result of adding a usually small multiple of the identity matrix to the matrix $A$.\n",
    "\n",
    "$$\n",
    "\\hat{A} = A + \\lambda I\n",
    "$$\n",
    "\n",
    "Properties\n",
    "- Only diagonal elements are modified.\n",
    "- Since shifting only applies to diagonal elements it can make previously identical rows in the matrix slightly different.\n",
    "  - This helps improve the conditioning of the matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4477c979",
   "metadata": {},
   "source": [
    "## 5.9 Diagonal and Trace\n",
    "\n",
    "### Trace\n",
    "Trace of the matrix $\\text{tr}(A)$ is the sum of the diagonal elements.\n",
    "\n",
    "$$\n",
    "\\text{tr}(A) = \\sum_{i=1}^M a_{i,i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "944003d1",
   "metadata": {},
   "source": [
    "## 5.10 Exercises\n",
    "\n",
    "3. Work out formulas for the count of elements based on matrix dimension $N$.\n",
    "- Number of elements: $N^2$\n",
    "- Number of elements along diagonal: $N$\n",
    "- Number of elements in upper / lower triangle: $N (N-1) / 2$\n",
    "    - ... combinatorics: form unique pairs from $N$ and $N-1$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "a & b & c & d \\\\\n",
    "e & f & g & h \\\\\n",
    "i & j & k & l \\\\\n",
    "m & n & o & p \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "- Number of elements: $4^2 = 16$\n",
    "- Number of elements along diagonal: $4$\n",
    "- Number of elements in upper / lower triangle: $4 (4-1) / 2 = 6$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d9e12a2",
   "metadata": {},
   "source": [
    "## 5.12 Code Challenges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6630d0e",
   "metadata": {},
   "source": [
    "> The goal of this exercise is to create a matrix that contains the dot products between all pairs of columns in tow other matrices. First create two $4 \\times 2$ matrices of random numbers. Then write a double-for loop in which you compute the dot products between each column of both matrices. The $i,j$ element of the resulting matrix will be the dot product between column $i$ of the first matrix and column $j$ of th second matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c0a1e425",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dot_pairwise_column(A, B):\n",
    "    \"\"\"\n",
    "    dot_pairwise_column returns the pairwise dot products of columns of matrices A and B\n",
    "    \n",
    "    :param A: numpy.ndarray  Matrix A of dimension m \\times n.\n",
    "    :param B: numpy.ndarray  Matrix B of dimension m \\times n.\n",
    "    :return: numpy.ndarray   Matrix C with pairwise dot products along diagonals.\n",
    "    \"\"\"\n",
    "    assert A.shape == B.shape\n",
    "\n",
    "    ncol = A.shape[1]\n",
    "    C = np.zeros((ncol,ncol))\n",
    "    for i in range(ncol):\n",
    "        for j in range(ncol):\n",
    "            C[i,j] = np.dot(A[:,i], B[:,j])\n",
    "    return C\n",
    "\n",
    "\n",
    "A = np.stack((np.array([1,2,3,4]),np.array([5,6,7,8])), axis=1)\n",
    "B = np.stack((np.array([4,3,2,1]),np.array([8,7,6,5])), axis=1)\n",
    "C = dot_pairwise_column(A, B)\n",
    "expected = np.array([[20,60],[60,164]])\n",
    "np.testing.assert_equal(C, expected, \"dot_pairwise_column\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f61a9b94",
   "metadata": {},
   "source": [
    "> Create a symmetric matrix by starting with a dense random numbers matrix and applying three matrix operations: \n",
    "> - convert to triangular form\n",
    "> - transpose\n",
    "> - matrix addition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "932e653d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def symmetric(A):\n",
    "    \"\"\"\n",
    "    symmetric creates a symmetric matrix from A using only addition and transpose\n",
    "\n",
    "    :param A: numpy.ndarray  Matrix A\n",
    "    :return: numpy.ndarray   Symmetric matrix A.\n",
    "    \"\"\"\n",
    "    S = np.zeros_like(A)\n",
    "    L, U = np.tril(A), np.triu(A)\n",
    "    return A + L.T + U.T\n",
    "\n",
    "\n",
    "A = np.random.random((4,4))\n",
    "S = symmetric(A)\n",
    "np.testing.assert_equal(S == S.T, True, \"symmetric\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce43e1b0",
   "metadata": {},
   "source": [
    "> Create a diagonal matrix of size $4 \\times 8$ without using the `diag()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1440879d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 2. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 3. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 4. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "def diagn(diag, n):\n",
    "    \"\"\"\n",
    "    diagn returns a diagonal matrix having n columns \n",
    "    \n",
    "    :param diag: numpy.ndarray  Vector of elements along the diagonal.\n",
    "    :param n: int               Number of columns.\n",
    "    :return: numpy.ndarray      Diagonal matrix.\n",
    "    \"\"\"\n",
    "    assert diag.ndim == 1\n",
    "\n",
    "    m = diag.shape[0]\n",
    "    D = np.zeros((m, n))\n",
    "    for i in range(m):\n",
    "        D[i,i] = diag[i]\n",
    "    return D\n",
    "\n",
    "\n",
    "D = diagn(np.array([1,2,3,4]), 8)\n",
    "v1 = np.array([1,0,0,0])\n",
    "v2 = np.array([0,2,0,0])\n",
    "v3 = np.array([0,0,3,0])\n",
    "v4 = np.array([0,0,0,4])\n",
    "z = np.zeros(4)\n",
    "expected = np.stack((v1,v2,v3,v4,z,z,z,z), axis=1)\n",
    "print(expected)\n",
    "np.testing.assert_equal(D, expected, \"diagn\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "425035c2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}