07_linear_algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Algebra\n",
    "\n",
    "This will cover linear algebra - a branch of mathematics you will make extensive use of in your degree. We will show how some linear algebra concepts can be implemented in python using `numpy`.\n",
    "\n",
    "We will cover the following:\n",
    "* Vectors & matrices in `numpy`\n",
    "* Dot products\n",
    "* Cross products\n",
    "* Matrix multiplication\n",
    "* Matrix inversion\n",
    "* Eigendecomposition\n",
    "* Singular Value Decomposition\n",
    "\n",
    "First, let's import the libraries we'll be using in this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we covered briefly in the numpy notebook, vectors and matrices are represented in numpy as arrays of the right dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "vector = np.array([1,2,3])\n",
    "matrix = np.array([\n",
    "    [1, 0, 0],\n",
    "    [0, 1, 0],\n",
    "    [0, 0, 1]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You could implement the dot and cross products yourself in code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n",
      "[0 0 0]\n"
     ]
    }
   ],
   "source": [
    "def my_dot_product(a, b):\n",
    "    return np.sum(a * b)\n",
    "\n",
    "def my_cross_product(a, b):\n",
    "    return np.array([a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]])\n",
    "\n",
    "# The dot product of a vector with itself is its length squared:\n",
    "print(my_dot_product(vector, vector))\n",
    "\n",
    "# The cross product of any vector with itself is always the zero vector:\n",
    "print(my_cross_product(vector, vector))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But thankfully numpy has functions implementing these already:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14\n",
      "[0 0 0]\n"
     ]
    }
   ],
   "source": [
    "print(np.dot(vector, vector))\n",
    "print(np.cross(vector, vector))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Dot products and cosines\n",
    "\n",
    "One very useful property of the dot product allows us to find the angle between two vectors. For two vectors a and b,\n",
    "\n",
    "dot(a, b) = ||a|| ||b|| cos(theta)\n",
    "\n",
    "Here, ||a|| is the length of the vector a.\n",
    "\n",
    "This property is so fundamental that it's considered an [alternative definition of the dot product](https://en.wikipedia.org/wiki/Dot_product#Geometric_definition). It is widely used in lighting calculations for computer graphics.\n",
    "\n",
    "Use this property to find the angle between the two vectors (1,2,3) and (-3, 2, -1). You'll probably want to make use of the inverse cosine function (also known as arcos), which in numpy is `np.arccos()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vector_a = np.array([1,2,3])\n",
    "vector_b = np.array([-3,2,-1])\n",
    "\n",
    "# Work out angle between these vectors, and print it.\n",
    "# Your code here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross product\n",
    "\n",
    "Recall that the cross product of two vectors is another vector perpendicular to both, following the [right-hand rule](https://en.wikipedia.org/wiki/Right-hand_rule).\n",
    "\n",
    "This property is also widely used in computer graphics, where it is common to need a vector perpendicular to a triangle (the \"normal\" vector).\n",
    "\n",
    "Below a cross product is calculated, and the 3 vectors are plotted using matplotlib's 3D plotting. Try changing `vector_a` and `vector_b` and seeing how the result changes. When setting these values, try using the right hand rule to predict the orientation of the blue line, and see if your prediction is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.92, 'Two vectors (red, green) and their cross product (blue)')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAECCAYAAABKcpRcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAB5JklEQVR4nO29d3hc5Zn3/3lGvYy6ZFmSLUuWi2w1FxnTAmmEErBDSAIpBAhZlgQC+y6/N2STJaRRNptNSNhA8pIlS0JCQjXFxiSUAAnGBtuSVSyr2eoaSaMy6lOe3x8z53A0mq4ZWbbP57p8WdKcOnPOd55zP/f9vYWUEh0dHR2d8GI42Qego6Ojczqii6uOjo5OBNDFVUdHRycC6OKqo6OjEwF0cdXR0dGJALq46ujo6EQAXVzDiBAiTgjRIITIDdP2LhRCdIVjW4uJEOIbQoj7FnF/vxVC/DBSy7ut+7AQ4t9DWfd0RwhxnRDi7UXa15x7zd9nKoSQQoiSMOz3v4QQ/xzIsj7FVQgxrvnnEEJMaX7/wkIPNNIIIe4WQvx+EXf5T8CbUsq+RdznUuTXwBeFEDkn+0DCfcNLKf9ZSvmDcG1PZz5CiONCiI/5Wexk3Ws/Br4thIj1t6BPcZVSJiv/gA7gcs3fHg/TwS5ZhBDRQa5yE/A7H9uLWtgRhY9IHouUchrYA1wbqX0sRUK4XsKybqRYSterB3zea5FCStkLHAWu8Lds0GEBIUS8awSb5fr9O0IImxAixfX7D4UQP3P9nCqEeEwIMSCEOOFadt4+hRB5rm1maP62SQgxKISIcf1+gxCiUQgxLITYK4Qo1Cy7UQjxFyGEWQjRL4T4NyHExcC/AZ9zjbRrNPt63rVsixDiq5rt3C2EeEoI8XshxBhwnRBimxDiPSHEmGvb/+XlfVkJrAbe1fztt0KIh4QQu4UQE8CHXft/2vWetAshvqFZPsG1zrAQogGoDvKz+b9CiF4hRI8Q4kbto1AIx2IQQtwphGgVQgwJIf6sfD5CiFWubX9ZCNHh+py+7XY4bwCX+TjWB4QQna739X0hxPma1+527e8xIYRFCFEvhNiqeX2TEOKg67U/AfFe9lEKPAyc7boGRjQvpwshXnJt410hxGrNeus111OTEOKzmtfUx0/hCtsIIb4phOgDHvVyHF91XbsW4XyU3ez6+3HXurXAhBAiWghxhet8R4QQb7jOQdnON4UQ3a7tNAkhPur6e6DXqHK8/+b6zI4LzROol2uk1HUcI67jukKzfKZw3ktjQoj9OK9/5TXlGonW/O0NIcSNvt4XIcTvgJXAC67P7P96OI9595qLLNfnZhFC/E1oNMJtfffjmPN04+vzd/EGPq5tFSllQP+A48DHXD+/CXza9fMrQCtwiea1T7l+fgzYBRiBVcAx4Ctetv8a8FXN7z8GHnb9vBNoAUqBaOA7wD9crxmBXuBfcd5kRuAs12t3A79328/fgF+6lq0CBoCPapa3uvZnABKAd4AvuV5PBrZ7Of7LgHq3v/0WGAXOdW0vEXgfuAuIBYqBNuATruXvA94CMoAVQB3QFeDnczHQB2x07ed3gARKQjyW24F9QAEQB/wK+KPrtVWubf8/13tUCcwApZrj2QyYfRzvF4FM1+f5r65jj9d8DtPApUAUcC+wz/VaLHAC+BcgBrjK9Zn90Mt+rgPe9vC5mIFtrv0/Djzhei0J6ASud722GRgENmrW/aHr5wsBG3C/6z1K8LD/zwDdOL8oBVACFGruqcOuzzoBWAtMAB93ndv/xXndxwLrXMeVp/kMVrt+DvQaVY73v1zHe4Frf+u8XCNG1/7/zXUMHwEsmuWfAP7ses/KXOf5tts1Eq3Z/xvAjQG+Lx/zce14u9cswIdc5/aA9nNn7r2gHof7NeLv83ctcyVw0O89GaK4/gD4uWvnfcBtOIUhHpgCsnDeFDPABs02bgLe8LL9G4HXXD8L1wl+yPX7HjSi7PrgJ4FC4BrgkJdt3o1GXHFexHbAqPnbvcBvNcu/6baNN4HvAVl+3p8v4BIAtw/8Mc3vZwEdbst8C3jU9XMbcLHmtX8icHH9H+Beze8lzBfXYI6lEdeXjuv35ThFLJoPbpwCzev7gas1v68B7EFcX8NApeZz+KvmtQ3AlOvnDwE9gNC8/g+CF9dHNL9fChx1/fw54C235X8FfFezrlZcZ3F9KXjZ/17gNh/31A2a3/8d+LPbdd7t2k8JYAI+BsSEeI1eiFNckzR/+zPw716ukfNx3t8Gzd/+6Pp8olzXw3rNa/cQuLj6e198iau3e+0Jze/JOO/1Fa7fAxVXn5+/6/ePA23+rulQswX+5vqgNgNHgL/g/BbcDrRIKQdxCqwyylA4AeR72eZTOB/f8nDeQBLnKA6cIvqA69FkBOeoQ7i2tQLnyDkQ8nCOpiw+jqnTbZ2v4BxRHBVCHBBCfNLLtodxftO7o91eIZCnnIfrXP4NWKY5Pu3y2vfOH+7rup9HsMdSCDyrea0R58W6TLMN7WTCJM4LWsGIcxTkESHEv7oeCUdd20/Fec1423a86xEzD+iWrqvcRTDvk79jLwTOcntfvgB4ywAZkM4Yszf8XZ/azyQPzblIKR2u1/OllC04nybuBkxCiCdc9woEfo0CDEspJzS/n3Dt19vxdLqOQ7t8PpCN84s21Os1mPvWHb/3mpRyHKdO5HlYzheBfP5GYMTfhkINov8D52PKp4C/SSkbXHGQy3AKLziH0lbXwTa4/rYS5zfxPKSUI0KIV4DP4nz8/6PmBuoEfiQ9TKK54irXeDlO6fZ7D5AhhDBqBNb9mOasI6VsBq4RzljxlcBTQohMtwsUoBYoFkJESyltXrbXCbRLKdd4Od5enBddvebYAqUX5yO8wgoPywRzLJ04R1V/d39BCLEqgOMpBWo8vSCc8dVvAh/F+XjnEEIM4/zC9EcvkC+EEJrrYyXeb1T3a8AfnTiv6Y8HuLy/7XeiiUX6Wb8HKFd+EUIInJ9jN4CU8g/AH4RzfuNXOMMRXwriGgVnrDlJ89pKnOEnb8ezQghh0AjsSpzhvQGco+AVOCd4lNcUlO0nAmOun7UC5et98feeervX1GteCJGMM7zW42H9CddxKbgfl7/P3+u1rSWkkauUchJnvO7rfCCm/8D52P831zJ2nI8cPxJCGF0i+H8AX6lRf8A5w/xp188KDwPfEkJsBHWi7DOu114EcoUQtwtn7ptRCHGW67V+YJXrokNK2ek6znuFc2KuAue3vtfMByHEF4UQ2a6La8T1Z7uH96QLaMYZx/PGfmDMNTGRIISIEkKUCSGUias/u84zXQhRANzqdiy/FUL81su2/wxc75qASMQZS/WFv2N5GOdnV+jad7YQYoefbWq5AGc4xxNGnDfmABAthLgLSAlwu++41v2GcE4AXYnv97wfKBABpM64eBFYK4T4khAixvWvWmgmloLkEeAOIcQW4aTE20QLzs/wMiHER4VzIvdfcYbW/iGEWCeE+IgQIg5nPHoK13UY6DWq4XtCiFjXl9wngSe9LPcuTiH6v6734ULgcpyP33bgGeBuIUSiEGID8GVlRSnlAM4vhS+6rq0bmCumvt6XfpxzAB7xca9dKoQ4z/VZ/wB413XPu3MYuNJ13CU4NUAhkM/f17WtspAigr/hDLrv1/xuxBn/UbgV54fTBryNUzD/x8c2n8cZq+uXUqrfDFLKZ3F+Sz8hnLP4dcAlrtcsOGMgl+N81GsGPuxaVblohoQQB10/X4MzHtQDPIszlvIXH8d0MVAvhBjHGSS/2sdj4K+AL3nbkOuCvBznRFo7ztH9IzgficEZNzvheu0V5qearADmjSRd296DMw7+Os5JiHdcL82EeCwP4Pw8XhFCWHBObp01f0vzEULE44xj/q+XRfbivDiP4TzfaTyHMTwd9yzO0dl1OB8PP4fzJvfGazifBPqEEIMBbN8CXARcjfMa6eODCaugkVI+CfwI57VvAZ7DOaLytGwTzom+X+D8PC7Hmf4469r/fa6/9wE5OMM4ENw12ofzfevBOaj4ZynlUU8LuvZ7Bc57bRDnRPC1muVvwRlO6cMZ83zUbRNfBf4/YAjnROs/Anxf7gW+43osv8PLeXi61/4AfBdnOGALzsd5T/wUZ6y8H+c1qg6u/H3+QojlOOcAnvOybRUxN3SlsxBco4pDOCeCesO87VicjyIVUkprAMuX4vwSinN7dIo4QohbcU4kzEuj0Tl5uEaev5dSFvhZdMkTyXvNz35/ArRKKX/pd1ldXE8fhBCfAl7CmU7yv4BDSrnzpB6UzpLhdBLXUwHdW+D04iacccxWnDG3m0/u4ejonLnoI1cdHR2dCKCPXHV0dHQigC6uOjo6OhFgyTnxLDJ6TERHJ7IEUhhyWqKPXHV0dHQigC6uOjo6OhFAF1cdHR2dCKCLq46Ojk4E0MVVR0dHJwLo4qqjo6MTAXRx1dHR0YkAurjq6OjoRABdXHV0dHQigC6uOjo6OhFAF1cdHR2dCKCLq46Ojk4E0MVVR0dHJwLo4qqjo6MTAXRx1dHR0YkAurjqhB2Hw8HMzAx2ux29jZDOmcqZbpatE0aklNjtdqxWK9PT0wghEEIQHR2t/jMYDAhxxvon65xB6OKqExaklMzOzuJwOBgaGqK/v5/09HRSU1MBsFqtAAghmJmZITU1VRdbndOaM7376xl98uHC4XAwOzuL3W6ntbWVsbEx8vPzGRsbY2RkBIfDQWpqqiq2hw4dYsuWLQD6yPb054z9MHVx1QkZKSU2m436+nqWLVtGc3MzWVlZFBUVYbVaVZG02+2MjIyo/ywWCwUFBaSnp5OSkoLB8EHoXxHbmJgYoqKidLE99TljPzxdXHVCQhsG2L9/P1arlbKyMtLT09XXvIni/v37KS4uVsVWCEFaWpoqtkIIpJRqzDYmJkYd2Sp/0zllOGM/LD3mqhM0NpsNq9WK3W6nqamJ6elpKisrSU9PD2h9IQRZWVlkZWUBznjsyMgIQ0NDtLa2EhUVpYqt0WhkdnaWmZkZAAwGAzExMerIVhdbnaWKLq46AaOEAWw2GxMTE9TV1VFQUIDD4SAqKirk7cbExJCdnU12djYAs7OzjIyMYDKZaGlpITo6mrS0NNLS0khJSWF2dpbZ2VkANWwQFxenxmx1dJYCurjqBITD4VBHq93d3XR3d1NWVobRaGR0dDSs+ayxsbHk5OSQk5MDOMV2eHiY/v5+WlpaiImJUcXWaDRSV1dHSUkJ8fHxGAyGeTFbHZ2TgS6uOj7R5q5arVYaGhqIjY1l27Zt6mhViZFGitjYWJYtW8ayZcsAmJmZYXh4mN7eXo4dO8bMzAx9fX1kZ2eTnJw8b2Sri63OyUAXVx2vSCnV0ero6CgNDQ2sXr2a3NzcOctFWlzdiYuLIzc3Vz2OQ4cOERMTQ1dXF+Pj48THx5Oenk5aWhqJiYnzxFaZINPFVieS6OKq4xEld9XhcHD8+HEGBwfZtGkTiYmJ85YNVlzDLcRRUVHk5OSwcuVKpJRMTU0xMjJCR0cHExMTJCQkqGKbkJDAzMzMnAmy2NhYoqKi1GwEHZ1woIurzhy0k1azs7McOXKE1NRUqqurvY7yFnvk6gshBImJiSQmJpKXl4eUksnJSUZGRjh+/DgTExMkJSWpYhsfH8/09LS6flRU1JyRrS62OqGii6uOipKfevDgQQoKCmhubmb9+vVkZmb6XC9YcV1MwRJCkJSURFJSEvn5+UgpmZiYYHh4mLa2NqampuaJ7fj4OK2traxbt04XW52Q0cVVB0AdrSrVVA6Hg61btxIXF+d33aUwcg1U9IQQJCcnk5yczIoVK5BSMj4+zvDwMC0tLUxPT5OUlITFYmF2dpa4uDh9ZKsTErq4nuFowwBTU1McOXIEIQRbtmwJSrBOprguZN9CCIxGI0ajUY3Zms1mxsbG1EwEo9GojmyFEExNTanvjS62Ot7QxfUMRslddTgc9PX1cfz4cTZu3Eh9fX1QInGyxTWcKCPbxMREKisrcTgcWCwWhoeHaWxsxGq1kpqaqubZ6mKr4w1dXM9AtLmrNpuNpqYmHA4H27ZtIzo6+EvidBJXBUUUDQYDqamppKamsmrVKhwOB2NjYwwPD9Pd3Y3dbp/j+OUutlrHL11szyx0cT3D0OauWiwW6uvrWblyJfn5+SHf+KebuPo6F4PBoI5ai4qKsNvtqth2dnbOs1e02+1qLHtgYICCggJ1ZKs7fp3e6OJ6BqHNXe3s7KS3t5eKigqSk5MXtF13cbXZbPT19ZGSkkJiYuIpJyCKI1cgREVFkZ6erprWaO0VT5w4gZRSLdM1mUwsX74cm80GON83bRhBF9vTC11czwC0k1ZWq5X6+nri4+PnlLAuBK24WiwWjhw5Qnp6OoODg0xNTZGcnKwKUHx8/IL3F2mCEVd3oqKiyMzMVNPXbDYbIyMjDA4OYrFYOHTo0Bx7RWUZ0I3DTzd0cT3N0fqujoyM0NjYSElJiVqnHw6EEDgcDrq6uujs7KS8vJzY2FhVdMfHxzGbzRw9elS1DzSZTKSnpxMTExO24wgXCxFXd6Kjo8nKyiIpKQmr1cr69etVsfVkr6h8AYIutqc6uriextjtdtra2sjJyaGzs5Ph4WE2b95MQkJCWPcjpaSrq4vk5GS2bduGwWBQa/m1qU6FhYU4HA7effddxsfH6erqwuFwqOKSlpYWlpH0UsThcKi+BoHYKypiqxjmgN6l4VRDF9fTEG0YoL+/n56eHjIzM9m6dWvARiWBjt4sFgvHjx8nIyOD8vJydV1vKC5VRUVFCCHUx2az2Ux7ezsGg0ENIbi3gFkswjly9bdNb/aKfX19NDc3q/aK6enpJCcnzxNbgPj4eH1kuwTRxfU0Q5u7OjAwwOjoKKWlpeTn5we8DYPBEJABdnd3NydOnFBnwLX4usm1bVyUx2alK4FWXI4dO0ZcXJwqtsnJyV63G05RWUxxdcebvWJPTw8Wi0V9P9LS0khOTqampob169cTExOjd2lYYujiepqgzV11OBwcO3aMqakpMjMz1YmTQPGXWmW322loaEBKybZt2+jr61MnZYI5Xk+4i8vU1BTDw8N0dHQwPj6u+gCkp6eTkJBwyoiHEhYIFnd7RcXxS7FXnJ6epr+/n6ysLNVecWZmBiHEPC9bXWwXF11cTwO0YYDJyUmOHDnC8uXLKS0tpb6+HofDEdT2lJGrJ8bHxzly5AgrVqxQc2ODFY1gbvCEhAQSEhJUhyvFdEXxATAajUxPTzM7O0tsbGxQx+GNkzly9YfyfixfvhwpJQcOHCA6OpoTJ07Ms1f05GWrG4cvHrq4nuIouatSSnp6eujo6KCsrEwdrfoSSm94G7n29PRw/PhxysvLMRqNc14LtogglKIDd9MVpTR1ZGRErTJTEvjT09NDqjZTji3c4hrqyNUXyki0oKCAgoICv/aKCQkJutguIrq4nqJowwB2u53GxkaAeSWs4RBXZft2u91jiazBYJgnlr4EKlzCpZSmJiQksGHDBqKiohgdHVXDCIA6GZSamhpwJsJSHrn6IhR7Rb1LQ+TQxfUURJu7OjY2RkNDA6tWrSIvL2/esqGIq3adiYkJamtr1dGRN4EI1s81EuWyUVFRZGRkkJGRAXzQslvJKY2OjlZHtUaj0at4RGrkutjxzkDsFZUCj7S0NOLi4tQuDR0dHaxatUoX2wWgi+sphraEtaOjg/7+fiorK0lKSvK4vKdRpT8U8evt7aW9vX1OmCFc+1gM3HNK3WfelV5b6enpJCUlRVT8pJQnXZw82Ssqjl/u9oomk4nCwkJVbKWUc0IIeksc/+jieoqgNVyxWq3U1dWRmJioJu17I5SRK0BzczMwP8zg6/gC5WQZvWhn3pVeW8PDw2p8UhnFRSJf9GSMXP0hhCAlJYWUlBS1wMNisTA0NMTMzAzvvffeHHtFg8GgGodfcsklvPPOO8Hu73+ATwImKWWZh9cF8ABwKTAJXCelPLjQ8zxZ6OJ6CqCEAf7xj3+wdu1ampqaWLNmjZp47gulNDVQJiYmGBoaoqCggLVr1wYkCO5iqUy0LEacMVS0vbaU+KTyyNzf38/09DQ2m00d2S40E2EpjFz9oY1hj42NUVFRMc9eMSUlhYMHD6px2iD5LfAg8JiX1y8B1rj+nQU85Pr/lEQX1yWOUmvucDiYmZmhra2NLVu2BGyAEszIta+vj7a2NtLT08nNzY1YJ4KlaFGofWROTExkZGSErKwsVVi0mQhpaWlBZyJEapIsEtjtdjXG6m6vODAwwIEDB+jo6OCss87i8ssv5zvf+U6gx/umEGKVj0V2AI9J54ntE0KkCSGWSyl7w3Bai44urksUbe7q9PQ0R44cwWAwUFlZGVBfK4VAxNXhcHD06FFmZmaorq6mubk5qNHuyRbLSOzbk7CMjIyoYQQhhDqqTU1N9TsqjURYIBLpXeAUV0/bjYqKIjc3l+985zs0Nzfz0ksv0draGs5d5wOdmt+7XH/TxVUnPGhLWE0mE62trWzYsIG2traghcRgMKi16J6YnJyktraW3NxcSktL5zzSB8pSGLlGuvzV3UrQarUyPDysGq7ExMSQnp5ORkaGxzJdKWXYTWkiJa7+Sp+VSrmkpCQqKirCuWtPH+LSesQJAl1clxDuuavNzc3MzMywbds2tXY8lLQqb0LW399PS0sLGzduJC0tTf17sHFad7F0OBx0d3cTHx9/SjpdBfIIHxMTM8dwZXp6Wu1GMD4+rlZKZWRkkJCQcMoUJsAHYQFvKJN/EaALWKH5vQDoicSOFgNdXJcI2myAiYkJ6urqyM/PV0eT4Bw92e32oLbrSZAdDgdNTU1MTU2pwu2+Tqgj0enpaWpqakhNTWViYoK2tjaio6PV/FNlVLfUYq4LJT4+nuXLl6tlqZOTkwwPD9Pa2srU1BSAGrMNJqzji5MlrpOTk15T/xbI88AtQogncE5kjZ6q8VbQxXVJoM1d7enpUQ2n3UtMQ6220q6jDQOsX7/e42gq1Mf8oaEhjh49SmlpKcnJyepobWZmBrPZrJqvKFaIubm5S7IzwUJHmdpKKaUstampSTW8sdlspKSkkJGRQVpaWsiG4Ysdc1VQwgLBIoT4I3AhkCWE6AK+C8QASCkfBnbjTMNqwZmKdX3QO1lC6OJ6EtFOWtlsNhoaGoiOjvaaW7rQaiuTyURzc/O8MEA49jM2Nsbk5CRbt24lLi5uTqpOXFzcnFFdTU0Ndrtd7UyQmppKRkbGgvwAwkm4H+GFEMTExJCamkpWVhYOh2NOma7SZysjIyOoMt2TGXMNJSwgpbzGz+sS+HrQG16inPwr+QxFW8I6OjpKY2MjRUVFLF++3Os6UVFRIYmrImQTExNUV1f7zdkMZuRqtVrVUdn27dv93uyKh2tubi6rV69Wz99sNnPixAl1Fj4jI+O0M8tWzkVrCA4f9NkaGhpSW78o74GvMt1APHdD4STGXE8rdHE9CUxPT2O1WomKiuLEiRMMDAxQVVVFYmKiz/UUoQwGq9VKf38/q1atYt26dQEXBQQi4mNjY9TV1ZGXl8fo6GjAQqg9BnehUWbhtWbZSrx2sTrJRiIe7EuwvRmGuxtkZ2RkzCnTjWRYwFeoYnJyUhfXANDFdRFRwgD9/f2MjIwwPj6O0Wikuro6oJsk2Md1k8lEU1MTRqORoqKigNcLZEJLaUZYWVkJwMjISMDbB+8C5j4LPzU1hdlspq2tjcnJSYxGoyq24fJv9cTJnNn3ZhjubiMYipduIPgbuY6PjwdUHXimo4vrIqHNXZ2YmKCnp4eKigp1tBIIgYqrw+GgubmZ8fFxysvL6ezs9LuOFl8jV639YHV1NdHR0UxOTs4TS1+hhWDNsvPz89USVYvFgtlspq6uDrvdTlpampplEa5OskvNctCbYXhfX59qrKKMbMPxhePvi0APCwSGLq4Rxr39SmtrK2azmZycnKCEFQJLxZqenqa2tpasrCw2b97M1NRUyHFad5RMg7y8PFasWKGKRSipVaGaZStGI6tWrVKrpkwmE7W1tXMS+Y1GY8hiFgkfgHBVaGltBGNiYpiZmSEtLW3OF44yQRhKmS7oMddwoYtrBNHmrk5NTVFXV0d2djalpaV0dXUFvT1/MdeBgQGOHTtGaWmp6mkabEGAso67+Cnb9pRpEErqVjhQqqYSEhIoLy/H4XAwPDxMV1cXFotFfXxWEvlPJpES7KioKFJTU0lNTVXLdJVMBKVMV8lESElJCWgCLJA8V11c/aOLa4TQtl/p7++nra1NFabx8fGgJ6bAeymrw+GgpaWFsbExNRVKu85CxFVKSUtLCyMjI14zDZZK+au7peDk5CRms1n1KtW2gPEVQlhqYQFveHp892YYPjAwQEtLi2oY7mt0H8jINUJFBKcVuriGGW3uqmKIYrPZ5lRChZJSpaznLspKGCAzM5MtW7bMu1kWkhs7OztLbW0tqampbN261WfblqVWcaVN5Ff6bSkjOiUGrc0t1YrUqVKq6nA4/D72ezMMV7rHxsfHqznGSjaGvxQvPSwQGLq4hhFt7ur4+Dh1dXWsWLFiXnuUUMpYYb5QDg4O0tTUxPr161VDEX/rBIIQgqmpKQ4cOBCQb+xSGLn6w1PKlxKvbW5uVlO+0tPTI+JgtVQE25thuJKNkZyczOTkpM9uuhMTE/OqB3Xmo4trmNBOWnV3d9Pd3U1FRYXHb/iFiqv2Ud09DOBtnUBRylgHBgY466yz/ObewtIcufrDfUSnTXcaGRkhISEBh8NBRkZGWLwAIjVyXcg2vRmG19bW0tLSwuzsrFqmqw2lhDJyFUJcjLPLQBTwiJTyPrfXU4HfAytx6tJ/SikfDfnklgC6uC4Q9xLW+vp6YmNj2bZtm9dHq1CKAcApylarlffee4/09HSfj+oKwYyW7HY79fX1zM7OkpubG5CwKvtY6iNXf2jTndrb2zEYDMzMzKheAEoX2fT09JCqopbKyNUXimF4bGwsVVVVagNMs9lMV1cXDoeDV155RR1EBLHdKOC/gY/jdL46IIR4XkrZoFns60CDlPJyIUQ20CSEeFxKGVLLg6WALq4LQJu7OjIyQmNjI6tXryY3N9fneqHeZBaLBZPJRFVVldcwQKgoXV5XrFhBQkICJpMp4HWXolgulMTERLKzs+ekfCkjW4PBoE4aBZryFSlxjaSdo9YwHJxlul1dXTz11FNcfPHFJCcn88orrwRyDNuAFillG4DL9WoHoBVXCRhdfbSSATNgC+8ZLS66uIaAkrva19eH2WwmOjqawcFBNm/eHJGUHyklra2tDAwMkJ6eHnZhVXxdy8rKSE1NZXh4+JQzyw4n7kLobpQ9OzuL2Wymu7sbi8VCQkLCnEkhTyzFsECwREdHs2PHDn7yk5+wb98+xsbGAhV3Tx0G3HtjPYjTcrAHMAKfk1IGP+u7hNDFNUi0YQC73U5vby95eXkBl7AGy8zMDLW1taSlpVFRUcGxY8fCtm1tJZc2zWqheaszMzPU1dUhhCAjI4PMzMyweZguBWJjYz2mfLW0tDA9Pe0xTgknt6Q2nChfPqmpqYGuEkiHgU8Ah4GPAKuBvwgh3pJSjoV8oCcZXVyDQJu7Ojg4yLFjx0hOTmbt2rUR2Z/ZbKaxsZF169aRlZXFzMxMSLFaTyiinZGRwebNm+eZqYSSKgYwOjpKXV0dJSUlREdHMzw8TGNjI1arVU19klKeUiNXX3hK+VLilJ2dnUgpSU9PV1PzwimGkRBXf59LiJ9bIB0Grgfuc9kOtggh2oH1wP5QdrgU0MU1QOx2uyqszc3NTExMUF5ezokTJ7yu8/3vR1FUJPnSlzwLlbebWEpJW1sbQ0NDczq9LkT0tAwPD9PQ0KCKtjuhPrYrZi6bNm0iLi4Om82G0Whk5cqVatzSbDZjMpkYGRkhNzd3ntNTKIRbqBcSH/UUp1Qcrt5//31iYmLUeO1CzzsS4uqvgGB2djYUg/MDwBohRBHQDVwNfN5tmQ7go8BbQohlwDqgLdgdLSV0cQ0Qg8HA5OQkdXV15Obmsm7dOr8jyT/9ycCWLZ7FVXGecr+5lMR9o9HI1q1b59w8oRYfKEgpOXHiBP39/T7jw8GWzDocDqamphgaGlKzJGy2uXMR2rilwWAgMTERKaXq9GQ0GsnMzCQ9PT0k85FINygMlejoaLKzszl+/DjV1dVMT09jNpvV805OTlbFNtjQib+OAaHgT7DHx8cDziJRkFLahBC3AHtxpmL9j5SyXgjxz67XHwZ+APxWCHEEZxjhm1LKwRBPY0mgi2uAWK1Wjhw5woYNG9RYU3R09DwR0ZKUBJOTnl9Tcl21F7Iyoly7dq2af6klFJ8A7fEraWL+4sPB9NBSemYZDAYqKioCFqWYmBiysrJUpyeLxcLQ0BDd3d04HA61RDOQttWnEvHx8eTl5annPT4+roZ/rFZrUKYrJ2PkuoAuBLtxtnHR/u1hzc89wEVBb3gJo4trgMTGxrJ9+/agKq0SE2Fy0rPYKOvGxMQgpaS9vZ2BgQG/I8pQcDgcHDhwgFWrVpGXl+d3+UDDAsqXwfr162lqagoqTun+u+J2VVRUpD5KK9VTSommYpgdaSKRNuUJJa/UaDRSWFiomq4oI1ulqsxbR4KTIa566Wvg6OIaBJ5EwRdJSTAx4fk1pZBgdnaWI0eOkJycHJGMg97eXqampqiurvbZN8v92HyNkKWUdHZ20tPT4/PLwJdI+RJv5VE6OztbLdHUzsZHuufWYomrO+6mK+4dCbRfMsp7Hu7jPImdX087dHENgmAnehITJSaT95HryMgIx48fD6h+P1iU9tlKalAwIz5f5+lwOGhoaMDhcFBdXR1SEnswgqAt0SwoKJjXc8tgMDA9PY3FYiE1NTVsYnMyxNUdbUcCT18yMzMz9Pf3hxyn9kSkOr+eiejiGkGcYYH5f1fc5MfGxiJSeKDEQXNycli/fj0HDx4MKlbrbeSqbDc3N5eVK1eGLEALKSJwN2CZnZ3l/fffp7u7m6amJpKSkhacW7sU08Q8fcns37+fyclJNU4dSgdZdyLV+fVMRBfXCOKc0JorQMrEmMPhYO3atWEX1qGhIY4ePTrHKSvYFC5P4qdMumiNuJcCsbGxxMbGsn79egwGAxMTEwwNDc3xBFAmiAIVnHCHBSIh1gaDgaioKIqKitQ49cjICIODg7S2ts7pypCcnBzw+egx1/Chi2sQeBIdZQbf06NUYqKcE3MdGRmhvr6ekpISxsbGQr7pPN382kkxbW6st+P2hbtZdkdHB319ffO2GyqRKn/VtkBRJoiU3Nq2traAc0wjIa6RDjO4d5CdmZnBbDbT0dGhPsor5+7rM9TFNXzo4rpAlHQsTzGvxETnhJbDIenoOEFfXx+bNm0iMTGRiYmJBdkOam8Aq9VKXV0d8fHxHifFQh25Ki5ZBoMhYuW9kcTdE8A9x1QpU83IyJhTphpuMYyEP6w/4uLiWL58OcuXL1fDUGazmaNHj6opX0p4RTspGIi4BpJxoqOL64LxlY6VlAR2u+C992owGuPYtm2bKlAL9XRVbgCLxcKRI0coLi726sYVirgq6Vv5+fmsWLHC/0oetuHrtZMR13TPMdXa6UkpVaGNRMVXJJL9A0U7ol+5cuW8SUEhhBpCsNvtPifH9P5ZgaOLaxB4EgxfImkwTAFGUlPzWLcuZ956obZ6Udbr6enh+PHjXk25PziO4MR1aGiIyclJqqur1YmjpUyonWS1jf2sVqvartpsNmO1WsnOzlYbIC70+CLR2SBUwfbUlUE598HBQeLj47FarXNavyhEwijbtcyFwM+AGGBQSnlBSCe3hNDFdYF4ElclTjk6KoF1JCXNT7NSjK+DRWlS2NLSgtVqZdu2bX5zPQMVV6UcdWBggMTExKCFNVARiVSDwoUQExNDTk4OOTk5zM7OUlBQwOTk5JzmhqHm1i51u0HtuTc3N2M0GrHb7bS1tTE1NYXRaCQ9PZ2UlJSgmxMGYpQthEgDfglcLKXsEEKENy/xJKGL6wJxL4HVlplu2LABUDIG5opJqGEBKSW1tbXk5+cHnA4ViLgqXRRiYmLYunUr+/btC+q4FMFcCvmhC0UIoXq0esutDcYs+1Qyyrbb7SQmJpKSkqK2frFYLJjNZu69915effVVZmZmmJ6e5oILLggk3S0Qo+zPA89IKTsApJSBO7UvYXRxDQJ/YYGxsTHq6uooKipi+fLltLY6Rxaecl1DafUyODjI8PAw69ato6CgIOD1/HkFTE5OUlNTozZTDAejo6NMTEyQmZk5b6R3qplle8qtVWK1FovFr/lKJCa0ImHaAvNHxNrS5AceeICrr76a888/nz179lBeXs7y5cv9bTIQo+y1QIwQ4g2cRtkPSCkfW/DJnGR0cV0gigNUR0cH3d3dVFZWqo9NCQlOAfFUAhvMyFXpRDA8PExOTk7QFTK+DF8GBgY4duyY2oUgVJTRscFgoKOjg56eHlJSUujs7MRgMKiz9sqxn0ri6o67WbZivqLNrc3MzFST+SM1oRUJcfWXLTAzM8MnP/lJvvrVrwa6yUCMsqOBLTgtBxOAd4QQ+6SU4XOGPwno4rpAhBAcP36c1NTUeU0JFQ30NHINdEJL6z2wZcsWjh07FvREmKewgOIZazab/XaQDRSHw0F9fT12u53NmzerAqDkXCopUNHR0SQlJWG1WuekQC0VgjXLdjdfGR4enpPMn5SUhM1mC2t44GSJawgTWoEYZXfhnMSaACaEEG8ClYAurmcK7jfG2NgYJ06cID09nY0bN85b/gNxnR9zDSQsoHX1X7ZsmbreQsXVZrNx5MgR4uPj2bJlS9hu0kOHDpGTk8OqVatwOBzqPt1zLltbWxkfH6empgYhBJmZmUE1+4s0CxHBqKioOcn809PTdHd3Mz4+zv79+73m1gbLyRLXqampYLMnAjHK3gU8KISIBmJxhg1+GsxOliK6uIaAlJKuri66urrU7qCeSEoKPSyguPpXVVXNCQOEksKlFdeJiQlqamoCth8MhLGxMcbGxigrK/Mbg1Nq5OPj4ykoKFDjl52dnYyPj2M0GlXxCZcZyckkPj6ezMxM7HY7JSUl6uRQV1cXgJpfmpKSEpRYRkpcA9luMBNpgRhlSykbhRAvA7WAA2e6Vl2o57BU0MU1SGw2Gw0NDRgMBrZt24bZbMZsNntcVvmCD0Zc7XY7DQ0NSCnnhRkgtIkwJX1L8UctLy8nJSXF73qBjOD6+vpoa2tTR2SBosRc3eOXimn2kSNHANRutykpKYs2qo1UhZbBYPCaW9vU1ERCQoI6ivc3Ooxkc0JfZcGhxMr9GWW7fv8x8OOgN76E0cU1CCYnJzlw4ACFhYXk5+cDzlQsXxVazvU8Zxm4j0CVWfuCggIKCgo8XuShhAWEEAwMDCCEmNPl1d86/vxYW1paGBsbo7q6mtra2jk3nr8KLW9/15pmW61WzGYzPT09HD16VJ2Vz8zMjOioNhLeAp6EUJtfqu0iq82tzczM9NiV4GR1foWlYcd4KqCLaxDEx8fPyQYA34/3ioVqIKlYyqjS36x9sGEBq9Wq5mZu27Yt4BtDO/vvjs1mo7a2lqSkJLVzbDCtYSCwbIGYmJg5fqbj4+MMDQ1RV1entoLJzMwMS9PGSBKIWHvqIquYzihdCbS5tSdDXJdyhsdSRBfXIIiKipqXBuVLXGNiICZGegwLKDeb0k1WGQH6G5Epj/iBMD4+Tm1trdoUMFiTak83kzK6LiwsnBezDfTmC2Xko52VX7VqldoKpq+vj8nJSerr69V0r4VmPkQiLBCsEGrFFJwpUMPDw2psOioqioSEBGZmZsKS6REI09PTi9Jm53RBF9cF4issAL6bFDocDt577z3S0tLYsmVL2Kqt4INYaHl5OTMzMwwNDfldR4sncVW8Yj2NroMpDAhHEYG2FczY2BhFRUVzck2VUW0oDQ6XouVgXFzcnNh0a2srU1NTc853oUbZ/j6TUDq/nsno4hoE3iq0/HWAnZiYv97IyAiTk5OsW7fOY6dXbwTS36q5uRmLxUJ1dTUxMTFYrdYFZRhoPV295cSezKor5ZFacX1SjKO1DQ4DnSjSbjNchPsRXghBbGwsycnJ5ObmzsmtbWlpITY2Vj1fd+MVX0Sq8+uZii6uQeIuIv47wMo5I1dFqHp7e0lMTAxKWP3tz2q1UltbS0pKihoLhdAnwaSUas8sKaVPT9fFHrn6QmscrfSeGhoa4tixY8zOzobUnWAhRNpbwD23Vum11dbWxuTkJCkpKWRmZpKenu4ztzaQAgK9f1bg6OK6QPzdNNo+Woo5SlRUFNXV1bz77rtB33jehFLxdV29erVacOBvHX/7mZ6epra2NqCeWUvVL0Dbe2rFihVqd4KhoSFaW1uJi4ubM8qLBJESV29fdAkJCeTn55Ofn4/D4VDT2zo7nSX+3nJrA+n8qo9cA0cX1wjj7EYgmJiYoLa2do45iqeuAv7wJJS9vb20t7d79XUNdiYfnF8EdXV1bNiwQXXy98VSGrn6wr07gTKqVTqqpqWlYbPZ/ApNMDgcjrC3AA/UuEWbWwsfeLf29vbS1NREYmKiOnHmL3yhd34NDl1cI0xSEgwOWqmpqaGsrGxO8r7yiB/MTaxNxXI4HDQ3NzMxMeHT19WXcYsnenp6GB0dpby8PCBhVfYR7lSsxSAhIUHNK1bSn/r7+zl48CDR0dGqEAcTu3RnsUeuvvCVWzs5OYnBYGBwcNBjbm2kjLJdy1UD+4DPSSmfCvrEliC6uAaJNxHxdAM5HA7sdgsWS7Q6uaQl1FJWu93O7OwsNTU1pKens2nTJp83bzBm2ceOHWNiYoLs7OygUnyCHbkuRZT0J6UXmdJzSzGN1sYugxmJLlWzbPfc2sHBQfr6+hgZGZmTW5uZmUlycnLQ4uqaG/BplO06jijgfpwlsqcNuriGAU/VTEq8MimpAocjgZiY+bmpoRhmKy5TBw4cYM2aNeTk+DdtD0RctZNhmzZtUiexAsVdXK1WKxaLxWvZ6lLrROAJbc8td8PsqKgoVXh8dZKFpTVy9YWUksTERIqLi4G5HWT37NnDX/7yF1avXk1fX5/Xfm1a9u/fD/6NsgFuBZ4GqsN4OicdXVzDgHsHWLPZTGNjI+vXr2fZsmQOHPB8Y4UiriaTCYvFwjnnnBNw/MufuCpmLtomh6F2jAVUx6vExEQ1fUd5vI6JiVmyI1dfuBtmu9so+pqRP1XE1T1EpXUzW79+PaOjo3R3d/OFL3yBO+64g0suucTn9rq7u8GPUbYQIh/4FPARdHE9s/HVjUDpQWUymdiyZQvx8fEkJUmvRQTBmLA4HA6ampqYnp4mOTk5qIkFX0KpmGW7m7kEG0NVlh8cHKSpqYmNGzcSHR2NEEKN69XW1gKouaanclsYdxvFsbExhoaG6OjomPc4vVTDAp626S3+r1SEfe5zn+Oqq64KaHterh/3P/4M+KaU0n6qXgve0MU1DERFRTEzM8PRo0eJi4ubkw+akOB0xZIS3K+dQEeuMzMz1NTUkJWVxbp164Lub+XNLFtpRuip7DbYSTAhhDqq3rRpEzExMeoXjpIKpaQGdXZ2MjQ0pPqbKqlQ4Z5RXyy0nWSLi4tVG8WOjg7Gx8fVLxGj0Rg2c/BIjVx9fQYTExMYjcaAt+fKivFnlL0VeMIlrFnApUIIm5TyuYB3tEQ5Na/mJYaUUs0xda+3T0oCh0MwMwPx8XPXC2RCSzHMXrdunZokHizuIwK73a7m227dutXjTRpM+pbD4WBwcJCoqCh1ck0IoQqJw+FASondblfzTh0OB8XFxYyPjzM8PExHR8ecNKlgZueXSuaBgruNYl1dnRrTBtTUp4XYKEZKXH1NYgZbRFBdXQ1+jLKllEXKz0KI3wIvng7CCrq4Bo37zdDT04PZbGbt2rUezaeVa3FiwrO4BmKYvWnTprAluE9PT3P48GHy8vJYuXKl1+UCDQtYrc40s+joaPLz81Vh1b5PQ9NDGISBzIRMent76e3tZcOGDQghSE5OJjk5mRUrVqg5mMrsvGK5l56eviiVVJFA+ZLJy8sjJSUlbDaKkRJXX9sMNlvANQr2aZS9oANe4ujiGiIOh4OjR48yMzNDfn6+1xsjMdEpUJOT4J4y6i3m6nA4aGxsxGazeTTMDpWRkRHq6+spLS31a2wdyISWdiJMCEFzczP9/f1kZ2er7lSNg41c9cxVlKSX8F+b/4uxsTG2bNmiPn4q7WAcDsecHExwdjgYHh6mvb2dmJgYsrKyyMzMnOcPsNRjddrYcqA2ioF0JohkSa0nJicngwoLQGBG2Zq/XxfUxpc4uriGwPT0NDU1NSxbtozS0lKOHz/u1bzlA0/X+X20PI1ctdsuLCwM2w00OztLY2MjmzdvDsi8xN/IVXHI2rBhA8nJyWovrImJCQYHBzly5Aj7B/fz/aPfJz46nquXXY3NZqOqqmpe22pFRBQxt9vtOBwO1Ti7sLCQ2dlZhoeHVSNpRYTS0tIW9sYsAt5Gmb5sFJXqqXDZKAZCBJoTntHo4hoko6OjHD58eM7oz9fjvTYs4E5UVNQcb9bh4WEaGhpYv359wJVR/lCyDGw2G+edd17Ak0a+JrQ6Ozvp7u5WJ660YQDlMf+vI3/lzro7WZ2ymm+t+ha55GKz2TCZTF6b8ykCpBVbZyGGndjY2DkjvtHRUdUFSmkCmJmZSbx77GUJEGhWhNZGUUrJxMSEVxvFSKC7YoUXXVyDJDExcZ7tnrtIzl3ef5NCKSWdnZ309PQEPLIM5IZV4qHp6ekkJCQENRvvaUJLSqmGQhTXLff4qt1h51tvfIv/fv+/+Vjhx7gt7zYqSyvJyMhgbGyMgYEBNQk/MzOT7Oxsr5NXyqg2OjpanRRTxDYtLU0dtb7//vvY7XY1lKLEMINt+hcplB5awaCNR3uyUVS+UIKxUfSHv5jr9PT0kvzyWqro4hoksbGx826UqKgopqenPS6vhAWmpua/pohyfX09DoeD6urqgOKr/vpbwQddCBSXrP7+fr/b9bQPBZvNRk1NDampqZSUlKjLaI/BMmvhuheu4+W2l/nKxq/wufTPUVleqY52tAYi09PTqmHK1NQUaWlpZGVleZ28Um76qKgoYmJiVJEdGxsDYPny5ers/MjIiPponZSUpD5an6xust56aAWDu43i/v371acSq9VKWlqaOqoNNUbvL+aqtPPRCQxdXIMkWMPsD8IC82OuVquVvr4+iouL/Vr6ue/P12yxyWSipaWF8vLyORMQwSTta9vJKK1dVq1apT6yugtr51gnn3760xwdOsrdW+/mwuQLqays9Cpo8fHxc2zxtGbPcXFxZGdnk5WV5XWkZDAYMJvNtLS0UFVVpebVOhyOOe1RFG9TZcJIEVqj0bhoE2GhjFx9oYjcihUr5tgoat+/UGwUfYUFQu38eiaji2sY8NXqxVtYQHEiMhqNFBYWBrU/JcvA/TFfSklbWxtms5mtW7fOEbZARrtalOWVUt6NGzeqOY7uon6g5wCfffazzNhn+MVZv2Bz2mY2bNgQ8AjKYDDMsQGcmJhgaGiIhoYGrFYrGRkZZGdnk5qaqh5/V1cXfX19bN68WT1PZX/KqFZKOcfbVBGhrq4uLBYLRqNRFSElBhwJAQnHyFWLu1gHYqOYkZHhN6UtkPSupZ6ZsZTQxTUMBDKhpYQFtC1TysrK6OjoCHp/ntKk7HY7R44cIS4uji1btsy7SXx1c/W2j+HhYfr7+6mqqlLDIe4319NHn+af9vwTyxKX8ZOKn1CZV6mmZoWK4tSkxBqVvNDGxkaSk5NV4dy0aZPXEIJ2UkwpYABUMTUYDIyPj2M2m+ns7FQFXvEOCCeL3fDQk42i2WxWU9qUmLSnWLe34wz36PtMQBfXMBBIe+2JCaFWRhkMBqqrq5mdnQ2pLbS7uE5NTXH48OE5Rtz+1vGFlJKenh4mJiY466yzPE5cSSn5j33/wfff/j7bl2/nm6u+SUVJRUBuScEQHR2t5r7a7XZqa2uxWq0IITh48KAah1TSwdzRxmrhg1Gtw+GYY7enpEGdOHGCiYkJjh49GpK9oCfCLUzBfklqwyRKrFtbqKGMan0xNTWlG2UHiS6uQeLpJlFcsTyhXI+jo1YOHDhAfn4+K1Y4y62DMW7Roi2bVdK3NmzY4PMGCVRcbTYbtbW16kjOk7DO2Gb42t6v8UTDE3y65NN8JecrVG6sjGjOqVI+mp2drVaWzc7OMjg4SHt7OxMTE6SmpqqFBt4ef91Htco/ZcIoIyOD6elpcnJyVNcr5TVltBcKJ0tc3XGPdWttFCcnJzlx4oRHG8VQOr++/PLLXHLJJU14McoWQnwB+KayC+BmKWVNSCe2BNHFNQTcZ9J9jVyjoiAuTtLWZmLdunVzBDAUy0H4QJSVfFPFgcvfOv7EVRkBFxYWkpKSQmNjI/v37yczM5OsrCzS0tIYnBrkmueu4Z3ud/jXqn/lUuOlVFZWhi0dyBNKYUVRUdEc/9rY2Nh5fquK2CoVXVlZWV5FwZPQNjU1kZaWhtFoJDk5mVWrVqkFDNoYplLAEGhceamIqxZ3G8V3332X2NhYjzaKoRhlf/3rXwe4BO9G2e3ABVLKYSHEJcCvcbMkPJXRxTUMeBNJxXkqNnY1qanLSU+fO1kSSuNAcN6obW1taqPDQNO3fO1LOwJWbqLNmzdjt9sZGhqit7eXvQf38u+N/87Q7BD/Uf0fnJ1yNuXl5WFzevKExWKhrq6O0tJSnyNjd6GYmppS7Q9nZmbIyMhQvyA8CZOUkvr6epKTkykuLp6TU+telquM9lpbW9W23YtZwBBJu0HFRlHb2HDfvn3cd999JCQkcPDgQaqqqvzuf//+/ZSUlNDa2urVKFtK+Q/NKvtwumadNujiGgY8jUqUBn+xsbGkpEQxPS0Bm9/1/DE7O8vAwAAZGRmUlZUFlVrlbSa8u7ubjo4OjxNXUVFR5OTkcGTyCP9S9y/ERcXxk7KfUCyKcTgc9PT0kJWVFZF4nJJaVFlZGfQjaUJCwpxUJbPZTH9/v1pWqvU/UEIOOTk5ashGSXfyVMCg5OsWFRUxPT3N8PAwR48exWq1zqmiilRO6GL4w2obGxYXFxMdHc2vf/1rfvrTn1JZWckdd9zhc3vd3d3qe+linlG2G18B9izgFJYcuriGgL+6eyUvVJlgSkrCq2F2MFgsFmpra0lLSyMnJycocfbm6drc3Mz4+DhbtmzxGF8F+H+H/h//+uq/sj5zPXevvZuyFWWsWLGCmZkZBgcH1Yohf6PDYOju7lYr1haa/B8VFTWvrFTxP7DZbMzOzvqdDFS2oy1gkFISFxenluUCahXVsWPHVG+AUJ5OfLEYXQjciYmJoaqqip/85CcBbS9Ao2wAhBAfximu5wW08VMEXVzDjOLsX1ZWplYjJSV5Ln8Nhr6+Ptra2qisrMRkMoXU2FC7js1m48iRIyQlJVFRUeGxMMDusHPnG3fyy/d/yccLP8438r5B+dpysrOzAefkiJLy4z46TEpKUgsBggkbKLm64+PjbN68OexWg9qy0pycHGpqasjLy2NiYoJ9+/aRkpKiTl55yxLwluqluFopoQmlweH09DTvvffenLLchcRgT4a4BhtzLSgooLOzc86fmG+UjRCiAngEuERKORTwDk4BdHENE1JKWltbMZvN85z9ExIUV6zQttvS0sLY2JjaQXZoaGhB4qpMXK1cuVI1QvFUyvrlF77M3ra93LjxRj6T9hkqyyu9Ws65jw7Hx8cZGBjg0KFDGAwGdXLJVzM/h8NBQ0MDMTExVFRURDSvcmxsjPr6esrKytRz0hrCBON/oJy/cg6K0MbHx7N8+XL6+vooLy9neHh4jo+rtq9YMJwK4lpdXU1zczO+jLKFECuBZ4AvSSmPhXTgSxhdXMOAzWZjenqamZkZjwn8SUkwMuJ9fW9J5sroMjExUTVKgdBSuBRx1Xq6Go1Gj9VDHaMdXPXMVRwdOsr3qr/HBUkXUFFREbDtndZKr7i4WA0ftLa2MjU1RXp6uuojoOxbSQHLzMwMumItWIaGhmhubqaqqmpOloMQYo4hjNb/YHJycs5xB5LqZbPZaGxsJCsra04VmhBCdbxSuhMor3nL19VyMnpyBduFIDo6mgcffJDLLrvMl1H2XUAm8EvXOduklFsXcBpLCl1cQ0B78SuG0XFxcZSUlHi8QJOSJN3d3jvAejLMUOK2hYWF8zocaOv+gznmoaEhhoeHqaysJC4uzmN8VVvK+uDZD1KZXElZWdmCHs/j4uLm+Qgo4ZPExETS0tLo7e2lqKhIjV1Gir6+Pjo6OgKK5fryP4iPj1dH456yBBwOB/X19aSkpFBUVKT+TRnVKn3FtCEVpeeWEprwVsDgb5QZCoGMXANp467l0ksvRUq5Vvs3rVG2lPJG4MYgD/WUQRfXBaDYv5WXl9PS0uJ1NOkrLKCMQrUXtmJErY3balEaIgaKlJLh4WFsNps6svYkrE8dfYp/2v1PLE9ezs+qfsbGZRtZvXp1WB/PtSM4KSUmk4mjR48SGxtLR0cHk5OTPiuuFkJHRweDg4Ns3rw56KorT/4Hg4OD1NfXqzaHiv+Bw+GgpqaG7OzsOTPm3goYADWkIoTAYrGoif2e+oqdrLCAXqEVHLq4hoASBx0ZGVHjq/6qtLxlC2hzZLW+A74KA4IJCyieA1JKVqxY4VFYpZTc/879/ODvP2D78u3cWXQnG4s2euwJFk6UevetW7eSlJQ0r+IqkMfwQFDi4VNTUwHlaAaCUjpbWFio+h90d3erZjM5OTk+S4E9TYopYqtMuBUWFnrsKyalDLrdij/8iWsoLV7OdHRxDYGOjg7sdruavgT+zFuk12wBZT1lMkdKOac1tycCLT5QmhEWFBQQExPDsWPHGBkZITs7m4yMDKcPrW2ar738Nf7U+CeuWnMVN2TfQMWGiogYmGjp6emhu7t7zuO5e8XVyMgIAwMDNDc3k5CQoGYfBNPyROlHFhUVFVRecDAo/gdpaWkcOnRIPf7Dhw8DLMj/QNudAJwTccePH8dsNjM4OKiOahdaIefPKHt8fFwfuQaJLq4hUFhYOE9I/Zm3TE6ClOB+bylG2/X19QH3zQqmJbcycQWwfft2RkdHGRgYcI7kDFPc1XAXBwcO8v9t/v/4ROInqKioiOhNpKRaWSwWn6lWWsMRKSWTk5MMDAxw5MgRHA4HWVlZZGdn+wwfKKP21NRUVq1aFdHsA+WLbO3atapJSnFxsUf/A+2Xm7dz95bqpfQVUyYMvfUVC3Z0rgi5N/T+WcGji2sIeDPM9iWuUgqmp53xVy02m00tO83Kygpo//7CAn19fbS3t1NRUUF8fPycMICSh9k42Mj1T19P/3g/3177bc42nE1OTo5qihwJIdKOIisrK4PyllUew1etWoXVamVwcJDjx48zPj7uUbCUFjfLly8nPz8/7OeiZWJigiNHjlBaWjovRu7N/6C1tZXY2NiA/A/gg1Ht4OAgIyMjrFixgqioKJYtWzanLFeZcEtISAiquaHdbve5nB4WCB5dXEMgeGesDwyzteLa29uL2Wxm7dq1AQsreA8LKLHF0dFRNm/e7HXi6q/tf+VLz3+J+Oh4fnX2r9iQuoHVq1er9fJKulR2dnZYqq3gg1SrjIyMBXe1jYmJmVMDrx2Nx8XFkZaWRl9fHyUlJerjdKRQvA+0+bLe0PofrFmzJmj/g8HBQdra2tRQivuoNi0tTTUUV8pyGxoasNvtflt26xNa4UcX1zARmKcrZGXNLTtV4qHB4M0sW/EyqKysVD1EfZWyfm/d9yjNK1XFThlh2e121Si7qakJo9Go1uKH4m06PT1NbW0tK1euDLvfq7thy+DgIA0NDcTFxdHe3s7Y2BhZWVkLroryxMjICEePHg3J+wCC8z8YGBigvb1d9X9Qzh18FzAoX0AjIyP09vbS1NSkFjBkZGSo2/IXc52YmNBHrkGii2uY8N0B1vn/5KTAZnM+rqakpLBp0yaOHz8edEGAe8xVseRbvnw5eXl5amGAeynrN1//Jg8dfIiLVl3ErctvpaykzGPuYlRU1JxmeBaLhYGBAY4fP05MTIw6sRTIJMr4+DhHjhxh/fr1EZ8kGx0dpaWlhc2bN5OcnIzValU7DVgsFlJSUlTBWmieqFJcUFVVFRY3LF/+BzMzMzgcDjZs2ODzi9hXBwYlfi2EYHJyErPZzJEjR9TXpqenfYrrzMzMSWvweKqii2sIeIu5eusAqzxNDQ1NMTx8kKKiIpYvX66uF0opq3LTjI2NqeKlpOm4j1jHZsa47sXr2Nu2l6+WfZWrUq+ioqyClJSUgM5VmURZvXq1+ijb2NiI1WpVJ5Y8NfxT+oRFepIMUOOYWm/ZmJgY1VRFKW1VRoDBfkloMZlMHD9+nE2bNkVEcLT+B3FxcXR2dpKfn09vby/Hjh0L2P8A5l5fyqhW21dMO6ptaGhQDcczMjLmbDvYzq8vv/wyt912G8eOHWvBs1G2AB4ALgUmgeuklAeDeZ+WOrq4hglfMVelSWFtbStf/GLZHFFbSClrf38/ra2t6sSV8pqWjtEOPv3Mp2kaauIH237AeQnnzVk+WLSPsjabjaGhITo6OrBYLHMmlkwmE52dnWzatCmo1KlQUNK6fImde2nr1NQUAwMD6peEYgiubYLobV89PT1s2rQpoj62yr56e3vVoof8/PwF+R+4FzAooYCMjAx6e3spKSlRc3Y7OjrUwonJycmgBgCKUfZf/vIXVq9evQHPRtmXAGtc/84CHuI0MsoGXVzDhi/D7NHRXmAVK1asJyUlat56s7OzQe1LebTr7OxU05k8xVf39+zns89+lln7LL8855eUJZZRVla24J5QCtHR0erIUJuX2tDgvIfCXd3liePHjzM8PBy0g1ZCQgIrV66c0wSxu7ubxsZGryPDzs5OBgYGvDZGDCddXV2YTCaqqqrm7MuT/4Fi++ju2+Av1Uvxqu3p6UFKSWxsLHFxcWpzSKWA4a677qKrq4uvfe1rfOpTn+Kiiy7yeeyKUbbLdHzWk1G26/fHpNObcJ8QIk0IsVxK2buQ920poYtrCASaiqU0JATnI7HVGgM4/K7nC2WbUkoqKyvV4/FVyvrzTT9nfdZ61qxZEzGxMxgM6iy9UvY5NDRETY2zJZIST/Q2ugoWZVJwdnaWysrKBWU0aJsgSikZGxtTY8xK/6yZmRmmp6fDVuHli87OTgYHB6msrPQr4lrbx2D9DwD6+/vp7+9XvzC0XrVKAcNjjz3Ghz/8Ya6++mpOnDjh9/gDNMrOBzrdlskHdHHVmYt7WEA7yVRR4Wyo56lKK5iY68zMDIcPHyY3N5eRkRGGh4fnpe1IKbnvnfv44d9/yPa87dy56k5KV5Z6NYIOF4qDV3p6upp9YDQa1R5Uyg0/NTU1pw4/FKFSqtliY2PZuHFjWL8whBCqA39JSQlTU1M0NDQwMTFBbGwsra2tavggEiJ74sQJ1Vwn2O378z9wD3309fXR09MzZ3TsaVLsgQceIDo6mgsuuCCg4wjQKNvTh+bdgf4URBfXEPHVpFBr65eRkUF/v3MZT/4CgY5clS4E69atIy0tjYSEBNXxXkmVSk5L5ta/3MqfGv/EZ9d+luuzr2fjuo3qzRYpZmZm1M4LykSdFm0ivZJy1Nvby9GjR4NO81Laayv5spFESkl7ezvJycls3rwZh8Oh9hNTjl0JH4Qj/tre3o7FYqGioiKi/geNjY1ER0djtVp9hjiEEDz00EMcPnyYAwcOBLzfAI2yu4AVfpY5pdHFNUwoIqn0o9q0aZOa+6hMlE9MBFfZpWAymWhpaVEnoqSUc5ylxsbGONp1lFufu5X6sXq+vv7r7MjYQVVlVcRn6cfHx6mrq5tT9ukL95Qj7SO4MoOfnZ3t8TF2dnaWmpoaCgoKPIp4OHE4HNTV1ZGUlERxcTFCCLWfmBI+UFLUOjo65qSvhfKet7W1MTExQVlZWURGxNrQR19fH8ePHycnJ4e6ujpgvv+BlJJHHnmEN954g2eeeSaorAjFKLu9vZ3i4uJYPBhlA88Dt7jisWcBo6dTvBV0cQ0bigGyyWSiurp6zihMyfTxNHL1lS2gdI9VbPI8TVwJIeix9nD936+nf7Kf/9z2n5RHOTuyNjQ0qGIVCZFVUq3KyspCqjv39AiuTIgpj7FKmpcSZikpKQmqmi0UAhkdu6eoaSeWguknplTVzczMRMxYRouSxbF161b1GnX3P9izZw8zMzM0NjbywgsvBJ3toRhlf+ITnwBoxLNR9m6caVgtOFOxrg/fWS4NhK9Ge2cAIZ+81WpVY6VKHbvFYuHCCy/0eIOkpcXyta/ZueeeuUKqmLZs2bJlzt8Vs2WDwcDatU6/YU8TV39p/wvXPn8tCTEJ/Hjzj1mXvI4NGzZgMBjUTrEmk4mZmRkyMzPJyckJS7VSb28vnZ2dqvF2uFHSvAYGBhgdHWV2dpbi4mIKCgoiOlNvs9moqakhNzc3ZE8CJfSh+AB46yemTMrZbDZKS0sXRVhPnDhBVVWV1zCGw+Hg3nvvZdeuXcTFxVFUVMSTTz65kGOL7EktYfSR6wIZHx+ntraW1atX09ra6vUidDYp9BwWcJ/Qmp2d5fDhw+Tk5FBQUOCxMADg14d+zR2v3sGGrA3cve5u1uaspaioSF0uNjZWTRa32+0MDQ2p1UppaWlqTmowj6HKaHpkZCQk0+lAUdK8YmNjsVgsrFu3DovFwoEDB+bYD4YziV953wsLCxfUESHQfmJdXV0AiyKsSujFX37uU089xdtvv82+fftITk5mcHAw4sd2uqKL6wJQOhFUVFRgNBppbW31umxCgvdsAW1YQJm4Wrt2Lenp6R6F1eawcefrd6qlrN/I+walRaU+6/a18UJ3r9RAO7U6HA6amprUNLBIpySZTCba29vZtGmTGoNVSkMHBgbmpXktJPShhB1Wr14d1rCDt35ihw8fxmazkZubi9lsntNPLNwoVWn+hPW5557jkUce4aWXXlLDPJEOwZzO6OIaIsePH6e/v39ep1dvdn1JSZKpqfnbUdp2AKrYaSeu3G+4sZkxvvzCl3ml/RX+qfyf+HTKpynfWO6xHYw33L1Sx8fHMZlMHDp0aI45s3ZSSUm1SktLi7g3KjhzJZXqJK0gaEtDi4qK1NCHNtYZbJrX5OSkmokRaf+D2NhYRkZGWLZsGcXFxeqXnNJPLNwjciWW6k9YX3rpJX7xi1+we/fuoK4lHe/oMdcQGR4eJj4+fs4NvG/fPqqrqz3GBM85J4acHMlzz80vkf373/9Ofn4+JpOJ8vJyoqOjPYYBtKWsd2+7m3PjzqWiomLBLvRalEmlgYEB7Ha7OinT3NzMypUrIz5Lr4QdRkdHKS8vDyq+qsQ6lThtIGleirHMxo0bA/JaWAhKfm5CQoKagaCgHZEPDg4C/jsY+GNoaIjW1tY5TlqeeOWVV7j33nvZvXt3JNL2ztiYgj5yDRGj0ei1G4EnQUhM9BxzdTgcTE9PY7FY2LRpE+B54kpbyvrQuQ+xIX6DKsThRFsWarVa6e7upqamhpiYGMbGxoiPjyctLS0iI1cpJU1NTTgcjpByPYNN8xodHaWhoWFRjGWU1C6j0ah2g9XiaUS+kH5iimuXP3OZ119/nR/96Ee89NJLEc+HPtPQxTWMKFVani7mxEQYGpr7NyVv02AwUFpa6nXi6snGJ7lpz03kJedx35b7KEkrYe3atRGPeY6Pj9Pb20t1dTUJCQlzkv/Dad8HH4hPYmJiWDwJvKV51dfXY7fbSUxMZHR0dE4+cqRwOBxz2s0Egrd+YkpZq69+YmazOSBhfeutt7jrrrt46aWXgm6breMfPSwQIna7fZ4LVm1tLUVFRR5Nha+5JpqjRwWHDjk9X5Usg5KSErWkMicnZ06+qLaU9ey8s7mz6E7WFqxl5cqVoR52wPT19dHR0eHRQUtr3zc0NER8fDw5OTkhxwqVLgVZWVmLdm4tLS0YjUYmJydDzpwIBCVnNjMzMyznpu0nNjg4qPYTy8rKUntqNTc3U1VV5TNF7p133uGOO+7ghRdeiHRp9BkbFtDFNUQ8iWt9fT35+fmqY5GWr3wlmrfeMnDs2Kza2qO8vJyEhASsVitDQ0OYTCamp6edN0q6kTvfuZM/N/6Zz677LNdnXc+GtYH32QoVKSUnTpzAbDZTUVERUNhBKZ5Q0na0Ji3+UNKfItGlwBN9fX10dnaquZ7aUaHZbA7rpJLdbqempkZNqYsESj+xwcFBRkdHsdlsrF27lmXLlnl9ojhw4AC33XYbu3btingJMbq4nrGEVVyPHj2qPiq7c+ut0Tz7rIG//72Zvr4+VbjcwwB2u52jnUe58S83UjtSy00lN/GprE9RVVUV8TYbSszTbrdTWloa0ihuZmZGnRCbnZ1VCxc8mWkrs/Rr1qxZlHhfV1cX/f39VFZWevzScJ9UEkKoZuDBxmTtdrtqshPpBong9LNobGykuLiY0dFRzGYzcXFx6heF8vRx6NAhbr75Zp599llWr14d8ePiDBZXPeYaIsF3gHUwPo7aPFDZhvt2moab+Mwrn8E0YeJn5/6MjThnsevq6khJSSEnJ8dnW+ZQUdpQp6SkzClECJa4uDjVAs/dTFtpepiens7ExAR1dXWLMksPqIUP7v6oWtwnlZScVPc0L38TejabjcOHD5Ofnx/x7Ar4oJeXkg+sFEAo4YP6+nrefPNNGhsbeffdd3nhhRcWS1jPaPSRa4g4HI55PbPa29uJi4sjLy9vzt+tViu33mrmt79dydjYBEJ4nrjSlrL+ZMtPWJ2wmo0bNxIVFaXGOU0m07zH14U6MikTa/n5+fOOPVwoXqPKqNBqtbJ69Wry8vIiVuUFH9TuT01NsXHjxpBjqoGmeVmtVjXMsZAqr0AZHR2lsbHRby+vffv28a1vfYvc3Fza29v51a9+xdlnnx3x40MfueqEA08j14mJCZevaxUAk5OS5OT5wvqrg7/ijtfuYGPWRr63/nuszlw9Z9Zc60DvKfE/JyeH7OzsoOv8JyYmOHLkSMQfzRWvUZvNxujoKKWlpYyMjPD+++8TExMT8vH7QglzSCkXbIoSSJpXWloajY2NrFq1alFm3xVhrays9Cmsx44d4/bbb+fxxx+nvLwcu90edGshneDRR66hrijlvPYs3d3dWK1WNd1maGiIo0ePUl5ezmOPpXDHHXG0t0+ive+0payfWPUJbiu4jTWFa4IaQU5NTWEymRgYGACc5aA5OTl+iwuUOF1ZWdmitE3u7OzEZDLNi3lqj9/hcMwpZw1VEJWE/bi4OEpKSiJaUTY1NaXa+MXGxpKbm+u1aWO4GBsbo6GhYU5DRk+0t7dzzTXX8L//+79qHvVCuOGGG3jxxRfn2BVqkVJy2223sXv3bhITEzly5MiW063xYKDoI9cwEh0dzZSrxrWjo4Pe3l617FCZD5maEiiari1lvaniJq40XsnG9RuDLsFMSEigsLCQwsJCdUJJ251VSfHS3uj9/f2qkUc4WkP7Qkqp+pVu2rRp3qO59viV5HnlUT4jI4OcnBy/jQO1KHmlSvw40hgMBkwmExUVzo66Q0NDnDhxgvHx8YikeQUqrB0dHXz+85/nN7/5TViEFeC6667jlltu4dprr/X4+p49e2hubqa5uZl3332Xs88++7RrPBgouriGkaioKGw2G42NjczOzrJp0yY1tvqBYbbz/xOjJ7jqmas4Zj7Gj87+Edujt1NRUbHghHbthJKSptPe3s7k5KQ6IaPMJm/ZsiWi8U5wCmtjYyNCCMrLy/0KpKeuBYp7vra7rLdJKSX9SenjFWkUw5c1a9aoZuG5ubnk5ubOM8gJR5qXxWJRq8p8CWt3dzdXX301Dz30ENXV1SHtyxMf+tCHOH78uNfXd+3axbXXXosQgu3btwOcdo0HA0UX1xDxJBJSSvr6+lixYgVr1qzB4XCo4qq0156cFLzbvY/PPfc5Zu2z/Oq8X7Emeg0VFRVhb9UcExPD8uXLWb58uWo52NjYyMzMDMuWLWN0dDSibkx2u31OyWewj8jucU5lQq+1tVW1HczOzlbfN2UyaTE6FYBTWA8fPuzV8MXdIEfr5hVKmpfFYqG+vt7vl3BfXx+f/exneeCBBzjnnHNCPr9Q8NKc8LRqPBgouriGicnJSZqamtRHXMXRShEU5V7Y3fg3fvbmDvKN+dxXcR9FxiLWr18f8VJWcBpcL1u2jFWrVqlCpe3BlZWVFbYUL6vVSm1tLcuWLQtLAr37hJ4iVIcOHSIqKor09HRMJhOrV68mOzs7DGfgm6mpKWpqali/fr3HohFPx+8rzUvpuuAt/KG00/EnrCaTiauuuoof//jHATcUDCcBNic8I9DFdQEovYbMZjONjY2sWbOG5uZmtSur9iZJjHfOzt7/5s85+8Ob+VbRtyjJLWHlypURt+9TUq3y8vLUhHbtiEqZ+W5vb1dLWbUjwmBRGhZGatbcXahGRkaora0lLi6OtrY2xsbGPMaZw4VS/LBhw4aQc3Tj4uLmGJmbzWZ6enpobGwkJSVFbXwYHR2tOnf5M5gZHBzkqquu4p577uFjH/tYqKe3IAJsTnhGoIvrAunq6qK7u1uduFqzZg29vb00NTWRlpZGTk4O6enptNbtAr7ABc1JPNHYScyOgyR+oRQZYWH1l2rlbnCilLIqI0JFaAOd9FL2txjeqMr+jh49SmVlJampqWopsdZNKicnx28vq2D2V1tbG9YMC19pXgaDQc3R9SWsw8PDfOYzn+G73/0uF198cViOKxSuuOIKHnzwQa6++mreffddOA0bDwaKnooV6opScuTIEaamptiwYYMaW1VGSspkRn9/PyMjI8TFpvCPNw9yUe0BNh5+i5hjxwCwb9uGfccO7Dt3IgN0TAqUhaZaTU9Pz0mRUjIPvN3kioXfYqV2WSwW6urqvO5PKVwwmUyMjIwE3cbbHWUEWV5eHlJDxmCZmJhQ2/2MjY1ht9vVcmLtqHx0dJRPf/rT3HHHHVx55ZURPaZrrrmGN954g8HBQZYtW8b3vvc9tZjmn//5n5FScsstt/Dyyy+TmJhIXV1dtZTyvYge1BJFF9cF0NfXh9Fo9GoVqO7EVSXU3d1NTEwMRqORfIuF7DffJOaFFzAcPgyAo7IS244d2HfsQK5fv5BDU1uk+EswDxQlRcq92aGSyzk0NERzc7Pf9KBwoZR8BpphoR0RDg0NERsbq44WAylcUIR8Mbxf4YMRslbIlVH5wMAA4+PjDA0NYbFYePTRR7n11lv53Oc+F/HjCoEztkJLF9cQsVqt/OIXv+Cyyy4jPz/f6yOnIqzanvRjY2OYTCaGhoZISEggb2aGZX//OzEvvECU81EKx/r12HfswHbFFcjKSggifNDR0cHAwEBEMhDgg86sJpOJ8fFx4uLimJ6eZsuWLRHpBOuOYgS9kC8Ope5+YGAAKaXPPlxKXmk4UuUCPbaamhqfTwAOh4O3336b7373u/T19bFt2za+/e1vU1VVFfHjCxJdXM9QQj75qakpfv3rX/Pss88yOzvL5Zdfzo4dOygsLFRHsHa7nfr6ehISEjxWCWnLWAcHB505nkDuO+8Q+8ILGN5+G+Fw4Cgqwn7FFdh37sSxdSv4EPJjx45htVrV9tqR5sSJE/T29mI0GhkbG1PNZTIzMyOyf5PJxPHjx/22LgkGpQ/XwMCAOnOvtCAfHR1VY7qLMSIPRFjBef199rOf5fOf/zw33HADDQ0NpKamRtqbNRR0cT1DWfDJSynp7+/nmWee4ZlnnmFsbIzLLruM7du38/zzz3P77bcHbDmnTCYNDAwQHR3N8uhocvfvJ/6llzC8/jrCasWRl+cU2h07cJx7LrhSp5Sc0qSkpLA4+ftDSklLSwvT09OqIYo2F3VoaIikpCTVRDscxQo9PT309PRQWVkZkRE5oOYDDwwMMDw8rPqj5ubmRvzLSknv8peFMD09zec//3l27tzJTTfdFJbP+uWXX+a2227Dbrdz4403cuedd855fXR0lC9+8Yt0dHRgs9m44447uP766wPZtC6uZyhhP/nBwUEefPBBfv7zn1NWVsYFF1zAjh07gu5Nr9Tbm0wmhBDkxseTd/Ag8Xv2EPXKK4jpaWRWFvbLL2f60ks5lJ7OshUrFmXk4nA4aGxsJDo6mrVr13otqHAflSueB6GMODs6OhgcHKSysjLsdoueMJvNNDU1qR1azWZzwC3IQyFQYZ2ZmeFLX/oSF110EbfeemtYhNVut7N27Vr+8pe/UFBQQHV1NX/84x/ZsGGDusw999zD6Ogo999/PwMDA6xbt46+vr5APsszVlz1VKwwk5GRQW1tLfv27SM7O5vnn3+e73//+3R2dnLRRRfxqU99So29+sLdL8BkMnFo40YcpaXk/Pu/k19bS+KePRj+/GeSH32Uc1NScFx2GfadO7F/9KMQoUdYpW1Jenr6nBCIO0IIjEYjRqOR1atXMzk5iclkUquTAjWXkVLS3t6OxWKhqqpqUUIdSkx38+bNxMXFsWzZMvXLQlu4EGyamjcUYS0tLfUprFarlRtuuIELL7wwbMIKsH//fkpKSiguLgbg6quvZteuXXPEVQiBxWJR34eMjIyIl06f6ugj10VibGyMF198kaeffpqWlhY+9rGPsXPnTo9GJr5Q4oP9/f1MT09jn5igcmCArL/9jaiXXkIMDyOTkrB/4hPOFK9PfALClBZltVrVYoSF+L4q5jImkwmbzTYnxcu93XRzc7MaQ450qANQiyn8xXSnp6fVc1BakGdnZwdduKCU0JaWlpKamup1OZvNxle+8hWqqqr4t3/7t7C+F0899RQvv/wyjzzyCAC/+93vePfdd3nwwQfVZSwWC1dccQVHjx7FYrHwpz/9icsuuyyQzesjV53IkpKSwuc//3k+//nPMz4+zp49e/jv//5v6uvr+fCHP8yOHTvYtm2b30fe2NhY8vPziYmJoa2tjfy1a2nLzqahuJis22+noKUF41/+QvQLLxD9zDPIuDjsH/2oU2gvuwxCTOxXDEqKi4sXXF7qyVxGccHSpngdPXqUqKioRRNWZbJMKQjxRXx8PCtWrGDFihVzDHImJibmdCzw9cWpCOv69et9Cqvdbufmm29mw4YNYRdW8Fyy6r6PvXv3UlVVxWuvvUZraysf//jHOf/88xeli8Spij5yPclMT0/zyiuv8OSTT3Lo0CHOP/98duzYwTnnnOP1sctTqpWSHtXf38/ExASZaWkUdHSQ+te/EvX88xi6upDR0TguuMCZ4vXJT0KATvlKXXugdfShokwmKZN6SpvtSHRldae/v5+Ojg61cWGoOBwOzGYzJpOJ0dFRry3IZ2ZmOHTokN9KNrvdzq233kpubi733ntvRL5k3nnnHe6++2727t0LwL333gvAt771LXWZyy67jDvvvJPzzz8fgI985CPcd999bNu2zd/mz9iRqy6uS4jZ2VleffVVnnzySfbv38/27dvZuXMn559/vtqptKWlhZmZGZ8tS7QiZbFYSE9LI7+3l/TXXyd61y4Mra1IIXCce65zRHvFFUgvE2FKsn5ZWdmiVCUpMd2MjAyMRiMmk4nh4WGMRqOa4hXuCa3e3l66u7upqqoKaxxRKVxQsifi4+PJzs4mJSWF+vp6v8LqcDj4l3/5F4xGI//5n/8ZsS8YJSPi1VdfJT8/n+rqav7whz+wceNGdZmbb76ZZcuWcffdd9Pf38/mzZupqakJpBuxLq5nKEv25K1WK3/729946qmnePvtt6mqquL48ePcdNNNXHnllUEZR2tHUqkpKRSMjJChCG1DAwD2rVs/KMN1TWwMDAzQ1tYWtiovf9hsNmpqauZ1THUXKU92g6HS09NDb2+vz8aF4WJiYoLe3l46OjqcxSN5eV5bkDscDr75zW8ihODnP/95xEfuu3fv5vbbb8dut3PDDTfw7W9/m4cffhhwlrX29PRw3XXX0dvbi5SSO++8ky9+8YuBbFoX1zOUU+Lk+/v7ufjii8nOzqanp4eKigp27tzJRz/60aAS26WUaq29MhrMn5gg6803iX7+eaIOOrtxOMrKGP7IR2jbtIm1O3cSE6ZkfV/Mzs5y+PBhCgsLfTb2U+wGlRSv6OhoNfMg2Oqwrq4ute3MYqR3zc7OcujQIdasWUNSUtK8FuTKqFZKyV133cX4+DgPP/zwomRIRBBdXM9QTomTf++99+jv7+eyyy7Dbrezb98+nn76af7617+ybt06du7cyUUXXRRUzbunhP98m43st97C8eSTJBw+jJASx9q12K+4AtvOnciqqqDKcANFmSxbvXp1II+Zc5iamlJn7ZUy1pycHL9lqh0dHQwNDVFRUbGowlpSUjLPnUyJlw8MDHDPPfcwODhIQkICzz///KJUhUUYXVzPUE7pk3c4HLz//vs89dRT7N27l6KiIq644gouueSSoGZxpZRYLBb6+/vp7u4mKiqKtUYjufv2Octw33oLYbfjKCz8oDrsrLO8luEGg+KNGg6LQiVNzWQyMTs767V/2IkTJxgZGaG8vHxRRoXKqHz16tU+O+xKKbnnnns4fPgwxcXFvPXWW7zxxhsRnURcBHRxPUM5bU5eacr35JNPsmfPHpYvX84VV1zBZZddFpBoORwO6uvriY+PJzc3V31kjY2NJTc6muXvvUfciy9ieO01xOwsjtzcD4T2vPNgARZ+GzduDHtKj81mU128lPSonJwchoeHGR8fD6iQIxxYrVYOHTpEcXGxz1G5lJKf/exnHD58mD/84Q/ExMSobmuh4q+kFeCNN97g9ttvV5tZ/u1vfwt5f17QxfUM5bQ8eSklDQ0NPPXUU7z00kukp6ezY8cOPvnJT3q8wW02G7W1tWRmZlJYWDjnNaWyymQyERUVRW5CAssPHSJ+926i9u5FTE0hMzOxX3YZth07cHz4wxBA7FPxfl0MCz+Hw8HQ0JCaS5uTk8OyZcsinuKlCGtRUZHP3GApJb/85S/5+9//zp///OewGNIEUtI6MjLCOeecw8svv8zKlSsxmUyR6Byhi+sZyml/8kqV01NPPcULL7xAQkICO3bs4PLLL2fZsmUMDw/T0tLCihUr/Db1U8yzlfjmMqOR/CNHSNi9m6iXX0aMjSFTUrBffLGzDPfjH/+geZgGs9nMsWPHFs1pSjGZmZ2dpbS0VI01m81mkpOT1RSvcKZhBSOsjzzyCK+88grPPPNM2CwbA8ld/eUvf0lPTw8//OEPw7JPL+jieoZyRp28lJLjx4/z9NNP89xzz6kG2A888AAf+chHgnoEVfwOlPLPnNRU8o8eJenll51luENDyMRE7Bdd5AwfXHIJpKSo6V1VVVWL4v2q2DDa7fZ55jlKrFmZ1IuNjVX9AhYyelS60BYWFvodCf72t79l165d7Nq1K6zpboGUtCrhgPr6eiwWC7fddhvXXntt2I7BhS6uZyhn7MnX1tZyzTXXcPHFF/P+++9jtVq5/PLLueKKK3wasnhi3kRSWhoFra0kv/IK0c8/j+jvR8bGMnnuuXRu20beTTcRHWB12EKQUtLU1ATAunXr/J6T0lF2YGAAg8Gg5tIGM7q22WwcOnSIlStX+kwpA/j973/Pn/70J1544YWwm3A/+eST7N27d4647t+/n1/84hfqMrfccgvvvfcer776KlNTU5x99tm89NJLrF27NpyHcsaKq+4tcIYSFRXFs88+y9q1a+d40n7jG9/AYrFw2WWXsWPHDo8m3+4ofgf5+fnYbDYGBgY4mp/P1NVXk3nzzRR0diKeeYbEvXtZ//rryP/8Txznn499505sl18OublhPz8ppWqLuGbNmoC+LJKSkkhKSmLVqlWqMUtDQwN2u10VWl9VajabjcOHDwckrH/+8595/PHHeemllyLS3cC9C2tXV9c8s52CggKysrLU8/7Qhz5ETU1NuMX1jEUfuQaAv1lXKSW33XYbu3fvJjExkd/+9rds3rw5Ige8GAwMDPDcc8/x9NNPMzQ0xCWXXMKOHTtYv359UCNau92uGppMTk6St3w5KwYGSH31VWd1WHOzswz37LPVzAO5cuWCj1+Z0IuLiwuLcbjValVHtFpzmZSUFHXbirAWFBSQ6+fL4tlnn+Xhhx/mxRdf9GnYshACKWltbGzklltuYe/evczOzrJt2zaeeOIJysrKwnkoZ+zIVRdXPwQy67p7925+8YtfsHv3bt59911uu+02pa3wKY/ZbOb555/nmWeeoaurS/Wk9eVtoKD0D5uamqK0tFStDhsbGyM9LY284WHS33iDmF27MNTVAWDfvNlZhrtjB3LNmqCPV0kpS0xMpLi4OOxGJ/N8G9LTyczM5MSJE+Tn5/udFHzppZf46U9/qmZxRBJ/Ja0AP/7xj3n00UcxGAzceOON3H777eE+DF1cz1D8nnwgs6433XQTF154Iddccw3gjO+98cYbfm+0U43R0VFefPFFnnnmGVpbW/nYxz7Gjh07PHrSKvFOKeW8Ea97y+uUlBTyp6bIfOMNol94gaj3nJ2YHRs2OEMHO3YgN270Wx3mcDioq6vDaDRSVFQU/jfAw/6GhoZobGwEUHNpvZnLvPLKK9x7773s3r3bZzHBacYZK656zNUP3d3drFixQv29oKBg3qjU0zLd3d2nnbimpqbyhS98gS984QuqJ+2DDz5IY2Oj6klbXV2N3W7n4MGDZGZmeozZGgwGMjMzyczMRErJyMgI/SYTjeeei/Hii1nuKsONefFFou+9l5h77sGxerVqLOPYvHme0DocjjkdEhYDKSUdHR2UlJSwfPlyRkdH1WyIhIQEtX9YTEwMr732Gj/60Y946aWXziRhPaPRxdUPgRgJB7LM6UZycjKf+cxn+MxnPsPU1BSvvPIKjz76KLfeeisGg4GLLrqI7373u37fByEE6enppKenz3G/aq6uJvGCC1huMJDzj38Q+8ILRD/wADH/9V84Cgo+6Ia7fTt2nNkPWVlZc77kIondblcdvJSJorS0NNLS0igpKVHNZX72s5+xa9cuzGYzu3btikSSvs4S5ZS221kMAp119bfM6YxSmPDwww+Tl5fHhz70IYaHhznnnHO47bbbeP3117FarX63I4QgNTWVNWvWcNZZZ1FUVMRYQgL7Kit59wc/oH3fPib/+79xVFQQ/ZvfEH/RRcSXlDD+pS9ReOwYKyKQdeAJRViXLVvmsbOvEILk5GSKi4s577zziI6O5vrrr+f222/n97///aIco87JR4+5+iGQWdeXXnqJBx98UJ3Q+sY3vsH+/fsjeuBLkZmZGd58800+/vGPA85Z9jfeeIOnn36at99+m61bt7Jz504uvPDCoJP03duOL0tMZNn77zPzhz+Q+e67GKamkBkZ2C+91Bk++MhHIAIetA6Hg5qaGrKzs/122j1w4AC33XYbzz//PCtdWRCL4Reg7Hv79u386U9/4qqrrgp5f2Hg9H6E84EurgHgb9ZVSsktt9zCyy+/TGJiIo8++ihbt271uj1/N8jjjz/O/fffDzgfvx966CEqKytDPcclgc1m4+233+app57ijTfeoLKyUvWkDbYyaWpqit7eXk6cOEFcXBwrsrLIU8pw9+xBjI4ijUZnk8adO7FfdBGEwb9AEdZAwg+HDh3i5ptv5rnnnlO7qi6UQDJXlOU+/vGPEx8fzw033KCL60lCF9dFJpAb5B//+AelpaWkp6ezZ88e7r777tMmtQtQPWmfeuopXn31VdavX8+OHTsC9qRVyktXrlxJamqqOqJ1OBxkp6ZScOwYiXv2EPXii4jBQWR8PPaPf9w5or3kEgjBwk+ZMMvMzPQrrEeOHOGrX/0qTz31VFgT8gPJXAH42c9+RkxMDAcOHOCTn/ykLq4nCX1Ca5EJpEf8Oeeco/68fft2urq6Fv04I0lUVBTnnnsu5557rupJ++STT/LjH/+Y4uJi1ZPW6KEluGKIsmrVKnVyaOXKlaxcuZLZ2VlMJhO1eXnYvvxlsm67jYL2dpJeeYWoXbuIfuEFZEwMjgsvxLZzp7MbbgCdbBU7x4yMDL/C2tDQwFe/+lWeeOKJsFc6BZq58uyzz/Laa69x4MCBsO5fJzj0Ca1Fxlvaljd+85vfcMkllyzGoZ0UDAYD1dXV/Md//AcHDx7krrvuorm5mUsvvZTPfe5zPP7444yMjABOizzFacrTrHtsbCwFBQVs3ryZTZs2EZ+cTMOyZbz1mc9Qt3s35hdfxHbzzYhjx4j7+tdJKC4m7pJLiP7VrxC9vR6PTxHWtLQ0NW7qjaamJq6//noef/zxeY/q4SCQrJTbb7+d+++/f1G6K+j4Rh+5LjLBpG29/vrr/OY3v+Htt9+O9GEtCQwGA1VVVVRVVfHDH/5Q9aTdsWMHSUlJdHV18dBDD/m08FOIiYkhLy+PvLw81Tj72PQ0kzt2kHn99eSZTKT+9a9EP/88sf/n/8D/+T/Yzzrrg+qwVavUooTU1FS/ubNtbW18+ctf5rHHHqO8vDxcb8kcAslKee+997j66qsBGBwcZPfu3URHR7Nz586IHJOOd/SY6yITaNystraWT33qU+zZs+eMN9Lo6uri4osv5pxzzuHIkSMkJiaqnrQ5OTlB+x1oy1czMjLIGxn5oO14bS0A1q98hfdvvDGgaq8TJ05w9dVX88gjj1BdXb2gc/VFIJkrWq677jo95noS0Ueui0x1dTXNzc20t7eTn5/PE088wR/+8Ic5y3R0dHDllVfyu9/97owXVoChoSEefvhhzjvvPKSUtLe38/TTT/PFL36R6OhorrjiCnbs2MHy5cv9Cm1UVBQ5OTnk5OSobcc77XbqLriAtB07WD45Scabb9IRExOQsHZ3d3PNNdfw0EMPRVRYAaKjo3nwwQf5xCc+oWaubNy4cZ5fgM7SQB+5ngT8pXbdeOONPP300+qjaHR0NO+56u29cQrmPy4YKSVdXV2q+bfNZuPyyy9nx44drFixIqgRrcPhcJbh9vfT19dHQkKC2lDQm0FNb28vV111FQ888AAf+tCHwnVapxtn7MhVF9fTgFM0/zGsSCnp6+vjmWee4dlnn2V8fFz1pA3UdlBrVZiVlUV/f/+cVjBZWVnqRFF/fz+f/vSn+fGPf8xHP/rRSJ/eqYwurmcop8XJn6L5jxHF3ZP20ksvZceOHV47EnjzgNW2HR8aGqKpqYmRkRGefPJJ7rnnHi6++OLFPrVTjTNWXPVUrNOAQNK7lPzHMyUul52dzVe/+lVefvll9u7dS1FREd/97nc5//zz+cEPfkBdXR0OhwP4oGtBbGzsvFGuEIKUlBTV72Dt2rU899xzjI+P8+CDD7Jv376TdYo6SxxdXE8D9PxH32RkZHDdddfxwgsv8Prrr1NWVsZ9993Heeedx1133cXVV19NU1OT35Y2Y2NjfOc73+Gb3/wmx44d46c//anfrgP+ePnll1m3bh0lJSXcd999815//PHHqaiooKKignPOOYeampoF7U9nEZFSnsn/Tgv+8Y9/yIsuukj9/Z577pH33HPPnGVWrVolCwsLZWFhoUxKSpLZ2dny2WefXeQjXVqMjY3Jiy66SFZUVMiysjJ52223yddee01aLBY5MTEx519fX58877zz5BNPPBG2/dtsNllcXCxbW1vlzMyMrKiokPX19XOW+fvf/y7NZrOUUsrdu3fLbdu2hW3/i8TJvsdP2r+TfgAn+d9pgdVqlUVFRbKtrU29Sevq6rwu/+Uvf1k++eSTi3iES5PGxkZ55513SofDIScnJ+Wzzz4rv/jFL8qysjJ58803y71798qxsTFpMpnkBRdcIB977LGw7j+QL0UtZrNZ5uXlhfUYFoGTfY+ftH96WOA0QJv/WFpaymc/+1k1/1FJ8QoVf4+tAG+88QZVVVVs3LiRCy64YEH7W0zWr1/PvffeixCChIQEdu7cye9+9zvee+89LrvsMv74xz+yfft2tm7dyrXXXsuXvvSlsO5fL4U+zTnZ6n6S/+n4IJDH1uHhYVlaWipPnDghpZSyv7//ZBxqxJidnZVPP/20dDgcYd/2n//8Z/mVr3xF/f2xxx6Tt9xyi8dlX3vtNbl+/Xo5ODgY9uOIMCf7HtdHrjpLD62DV2xsrOrgpeUPf/gDV155pWpqcrq1MYmJieHKK6+MSNueQDtY1NbWcuONN7Jr1y69/9YphC6uOl4J5LH12LFjDA8Pc+GFF7JlyxYee+yxxT7MUxZtKfTs7CxPPPEEV1xxxZxl9FLoUxfdW0DHK1L6T/Gy2Wy8//77vPrqq0xNTXH22Wezfft2XQgCIBCvgO9///sMDQ3xta99TV3HXym0ztJAF1cdrwTanDErK4ukpCSSkpL40Ic+RE1NjS6uAXLppZdy6aWXzvmbttDjkUce4ZFHHlnsw9IJA3pYQMcrgTy27tixg7feegubzcbk5CTvvvsupaWlJ+mIdXSWDrq46nglkBSv0tJSLr74YioqKti2bRs33ngjZWVlPrfrL71rdHSUyy+/nMrKSjZu3Mijjz4akfPT0YkkunGLzqISiIPXPffcw+joKPfffz8DAwOsW7eOvr6+oNtx6ywJdOMWHZ3FIJD0LiEEFosFKSXj4+NkZGQQHa1PD+icWujiqrOoBJLedcstt9DY2EheXh7l5eU88MADXg2rlxr+Qh5SSr7xjW9QUlJCRUUFBw8ePAlHqbMYnBpXrM5pQyDpXXv37qWqqoqenh4OHz7MLbfcwtjY2GIdYsjY7Xa+/vWvs2fPHhoaGvjjH/9IQ0PDnGX27NlDc3Mzzc3N/PrXv+bmm28+SUerE2l0cdVZVAJJ73r00UfVqqiSkhKKioo4evToYh9q0AQS8ti1axfXXnstQgi2b9/OyMgIvV7aeuuc2ujiqrOoBJLetXLlSl599VXA2U6lqamJ4uLik3G4QRGoaXkwZi06py76LIHOohJIVdK///u/c91111FeXo6Ukvvvv5+srKyTfOT+CSTkEcgyOqcHZ3oqls5phBDif4BPAiYp5bxkW+FUsQeAS4FJ4DopZdhmlIQQZwN3Syk/4fr9WwBSyns1y/wKeENK+UfX703AhVJKPTZwmqGHBXROJ34L+OoYeAmwxvXvn4CHwrz/A8AaIUSRECIWuBp43m2Z54FrhZPtwKgurKcnelhA57RBSvmmEGKVj0V2AI9J5+PaPiFEmhBiebjETUppE0LcAuwFooD/kVLWCyH+2fX6w8BunCPnFpyj5+vDsW+dpYcurjpnEvlAp+b3LtffwjZylFLuximg2r89rPlZAl8P1/50li56WEDnTMLTzJE+6aATEXRx1TmT6AJWaH4vAHpO0rHonObo4qpzJqFPJuksGnrMVee0QQjxR+BCIEsI0QV8F4gBfTJJZ/HR81x1dHR0IoAeFtDR0dGJALq46ujo6EQAXVx1dHR0IoAurjo6OjoRQBdXHR0dnQigi6uOjo5OBNDFVUdHRycC6OKqo6OjEwH+f7h3bDH3T9MfAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vector_a = np.array([1,0,0])\n",
    "vector_b = np.array([0,1,0])\n",
    "\n",
    "cross_product = np.cross(vector_a, vector_b)\n",
    "\n",
    "def plot_vector(vector, axes, color):\n",
    "    x_vals = np.array([0, vector[0]])\n",
    "    y_vals = np.array([0, vector[1]])\n",
    "    z_vals = np.array([0, vector[2]])\n",
    "    axes.plot(x_vals, y_vals, z_vals, color)\n",
    "\n",
    "axes = plt.axes(projection='3d')\n",
    "\n",
    "plot_vector(vector_a, axes, \"r\")\n",
    "plot_vector(vector_b, axes, \"g\")\n",
    "plot_vector(cross_product, axes, \"b\")\n",
    "plt.title(\"Two vectors (red, green) and their cross product (blue)\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Cross Product and Parallelograms\n",
    "\n",
    "Another property of the cross product of vectors a and b is that its magnitude (i.e. length) is equal to that of a parallelogram with the vectors a and b as two of the sides, as shown in this diagram (from [wikipedia](https://en.wikipedia.org/wiki/Cross_product#/media/File:Cross_product_parallelogram.svg))\n",
    "\n",
    "![Diagram](https://upload.wikimedia.org/wikipedia/commons/4/4e/Cross_product_parallelogram.svg)\n",
    "\n",
    "I'd like you to verify this for the vectors a and b below. Do this by dividing the parallelogram into two triangles, and working out the area of each. The length of the sides of these triangles will be ||a||, ||b|| and ||b - a||. To work out their area you may want to make use of [Heron's formula](https://en.wikipedia.org/wiki/Heron%27s_formula)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Area estimated with cross product 13.856406460551018\n",
      "Area estimated by finding triangle areas 0.0\n"
     ]
    }
   ],
   "source": [
    "vector_a = np.array([1,2,3])\n",
    "vector_b = np.array([-3,2,-1])\n",
    "\n",
    "cross_product = np.cross(vector_a, vector_b)\n",
    "cross_area = np.linalg.norm(cross_product) # This function works out the norm, i.e. the length of the vector in this case.\n",
    "\n",
    "print(\"Area estimated with cross product\", cross_area)\n",
    "\n",
    "# Work out the lengths of the 3 sides of the triangle\n",
    "# Find the area of each triangle using heron's formula\n",
    "# Double to get the parallelogram area.\n",
    "# Your code here (replace this line)...\n",
    "triangle_area = 0.0\n",
    "\n",
    "## This should be the same as cross_area above.\n",
    "print(\"Area estimated by finding triangle areas\", triangle_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mathematical Operators and Matrix Multiplication\n",
    "\n",
    "As a reminder, all the standard mathematical operators in numpy `+ - * /` work **element-wise**. This means that `+` is vector addition, but `*` is **not** matrix multiplication.\n",
    "\n",
    "To do matrix multiplication you must use the `np.matmul()` function or the `@` operator (these do the same thing)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sum of the vectors: [5 7 9]\n",
      "Broadcasted element-wise multiply:\n",
      " [[1 0 0]\n",
      " [0 2 0]\n",
      " [0 0 3]]\n",
      "Proper matrix multiplication, giving a 3-vector: [1 2 3]\n",
      "Proper matrix multiplication, giving a 3-vector: [1 2 3]\n"
     ]
    }
   ],
   "source": [
    "another_vector = np.array([4, 5, 6])\n",
    "vector_sum = vector + another_vector # This works fine, and does vector addition as you would expect\n",
    "print(\"Sum of the vectors:\", vector_sum)\n",
    "\n",
    "multiple = matrix * vector # This does not do matrix multiplication\n",
    "# This will use numpy's \"broadcasting\" feature, and make copies of the vector to make it a 3x3 matrix.\n",
    "# Then it will multiply both together, giving a 3x3 result\n",
    "print(\"Broadcasted element-wise multiply:\\n\", multiple)\n",
    "\n",
    "# To do real matrix multiplication, use np.matmul()\n",
    "matrix_multiply = np.matmul(matrix, vector)\n",
    "print(\"Proper matrix multiplication, giving a 3-vector:\", matrix_multiply)\n",
    "# You can also use the @ operator, which does the same thing\n",
    "matrix_multiply = matrix @ vector\n",
    "print(\"Proper matrix multiplication, giving a 3-vector:\", matrix_multiply)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy can also invert a matrix for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multiplication of matrix with its inverse:\n",
      " [[ 1.00000000e+00 -2.07850389e-16  4.33108961e-16]\n",
      " [ 1.42311088e-16  1.00000000e+00  3.78593168e-16]\n",
      " [-2.35269174e-16 -6.72015753e-16  1.00000000e+00]]\n",
      "Bad test for equality: False\n",
      "Good test for equality: True\n"
     ]
    }
   ],
   "source": [
    "matrix = np.random.rand(3,3)\n",
    "\n",
    "matrix_inverse = np.linalg.inv(matrix)\n",
    "\n",
    "# When multiplying the matrix by its inverse, the result should be close to the identity matrix\n",
    "# It won't be exactly the same due to floating point errors\n",
    "multiple = np.matmul(matrix_inverse, matrix)\n",
    "print(\"Multiplication of matrix with its inverse:\\n\", multiple)\n",
    "# It's really really close though, those non-diagonal entries are in the order of 10**-17.\n",
    "\n",
    "# Because it's not exactly the identity, this test won't work:\n",
    "print(\"Bad test for equality:\", np.all(multiple == np.eye(3)))\n",
    "\n",
    "# If you want to check if one numpy array is \"close enough\" to another, you can use np.allclose():\n",
    "print(\"Good test for equality:\", np.allclose(np.matmul(matrix, matrix_inverse), np.eye(3)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can only invert square matrices, so the following won't work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "ename": "LinAlgError",
     "evalue": "Last 2 dimensions of the array must be square",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-44-8c3fac153a39>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mmatrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mmatrix_inverse\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmatrix\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# This throws LinAlgError as the input isn't square.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m    538\u001b[0m     \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwrap\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_makearray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    539\u001b[0m     \u001b[0m_assert_stacked_2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 540\u001b[0;31m     \u001b[0m_assert_stacked_square\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    541\u001b[0m     \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresult_t\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_commonType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    542\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_assert_stacked_square\u001b[0;34m(*arrays)\u001b[0m\n\u001b[1;32m    201\u001b[0m         \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    202\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mm\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 203\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Last 2 dimensions of the array must be square'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    204\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    205\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_assert_finite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0marrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mLinAlgError\u001b[0m: Last 2 dimensions of the array must be square"
     ]
    }
   ],
   "source": [
    "matrix = np.random.rand(2,3)\n",
    "matrix_inverse = np.linalg.inv(matrix) # This throws LinAlgError as the input isn't square."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, as covered in the course, not all square matrices are invertible. The matrix below is not invertible because its rows are not linearly independent (row 1 is just two times row 0).\n",
    "\n",
    "As such the matrix is *singular* and can't be inverted. Trying to invert it will give a `LinAlgError` that will inform you the matrix is singular:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "LinAlgError",
     "evalue": "Singular matrix",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-33edef8d9a2c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mmatrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mmatrix_inv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmatrix\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmatrix_inv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m    543\u001b[0m     \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    544\u001b[0m     \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 545\u001b[0;31m     \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    546\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    547\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m     86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     87\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     90\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix"
     ]
    }
   ],
   "source": [
    "matrix = np.array([[1, 0], [2, 0]])\n",
    "matrix_inv = np.linalg.inv(matrix)\n",
    "print(matrix_inv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the course you discussed decomposition of a matrix into eigenvalues and eigenvectors, the *eigendecomposition*.\n",
    "\n",
    "Numpy can simply compute the eigenvectors and eigenvalues of any *diagonalisable* square matrix for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: [1. 2. 3.]\n",
      "Eigenvectors:\n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "# Let's start with a simple example:\n",
    "matrix = np.array([\n",
    "    [1, 0, 0], \n",
    "    [0, 2, 0], \n",
    "    [0, 0, 3]])\n",
    "\n",
    "# The eig function returns first an array containing the eigenvalues\n",
    "# then a second array containing the eigenvectors\n",
    "# For this matrix as you would expect the eigenvectors are \n",
    "eigen_decomp = np.linalg.eig(matrix)\n",
    "\n",
    "print(\"Eigenvalues:\", eigen_decomp[0])\n",
    "print(\"Eigenvectors:\\n\", eigen_decomp[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, you can compute the eigenvalues of any *diagonalisable* square matrix.\n",
    "\n",
    "This is a related concept to being *singular*, but not the same thing. For example, the following matrix is singular but is diagonalisable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([1., 0.]), array([[1., 0.],\n",
      "       [0., 1.]]))\n"
     ]
    },
    {
     "ename": "LinAlgError",
     "evalue": "Singular matrix",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-16-700386618963>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meig\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msingular_matrix\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msingular_matrix\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m    543\u001b[0m     \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    544\u001b[0m     \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 545\u001b[0;31m     \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    546\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    547\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m     86\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     87\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     89\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     90\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix"
     ]
    }
   ],
   "source": [
    "singular_matrix = np.array([\n",
    "    [1, 0],\n",
    "    [0, 0]\n",
    "])\n",
    "\n",
    "# The eigendecomposition works fine:\n",
    "print(np.linalg.eig(singular_matrix))\n",
    "\n",
    "# But this will throw an error as the matrix is singular:\n",
    "print(np.linalg.inv(singular_matrix))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we covered in the course, Singular Value Decomposition (SVD) is a more general matrix decomposition that works for any n by m matrix. Numpy can compute SVD of a matrix too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U:\n",
      " [[-0.19972938  0.89075334  0.40824829]\n",
      " [-0.51626998  0.25845432 -0.81649658]\n",
      " [-0.83281059 -0.37384469  0.40824829]]\n",
      "Sigma: [4.3218954  0.56676285]\n",
      "Vt:\n",
      " [[-0.5048459  -0.86320949]\n",
      " [-0.86320949  0.5048459 ]]\n",
      "Sigma (matrix):\n",
      " [[4.3218954  0.        ]\n",
      " [0.         0.56676285]\n",
      " [0.         0.        ]]\n",
      "Reconstruction:\n",
      " [[-5.55111512e-17  1.00000000e+00]\n",
      " [ 1.00000000e+00  2.00000000e+00]\n",
      " [ 2.00000000e+00  3.00000000e+00]]\n",
      "Reconstruction is close enough to original matrix.\n"
     ]
    }
   ],
   "source": [
    "matrix = np.array([\n",
    "    [0, 1],\n",
    "    [1, 2],\n",
    "    [2, 3]\n",
    "])\n",
    "\n",
    "u, s, vt = np.linalg.svd(matrix)\n",
    "print(\"U:\\n\", u)\n",
    "print(\"Sigma:\", s)\n",
    "print(\"Vt:\\n\", vt)\n",
    "\n",
    "# Note that the \"Sigma\" is returned as a vector rather than a matrix, so we can't just multiply\n",
    "# U * Sigma * Vt to get the original matrix back.\n",
    "\n",
    "# s is a vector containing the \"diagonal entries\" of the Sigma matrix.\n",
    "\n",
    "# We need to convert Sigma into a 3x2 matrix first:\n",
    "sigma = np.zeros((u.shape[1], s.shape[0]))\n",
    "for i in range(s.shape[0]):\n",
    "    sigma[i, i] = s[i]\n",
    "\n",
    "print(\"Sigma (matrix):\\n\", sigma)\n",
    "\n",
    "# Now we can reconstruct the original matrix by multiplying (matrix multiplication)\n",
    "reconstruction = u @ sigma @ vt\n",
    "\n",
    "print(\"Reconstruction:\\n\", reconstruction)\n",
    "\n",
    "# We can check it's \"close enough\" to the original matrix with np.allclose():\n",
    "if np.allclose(reconstruction, matrix):\n",
    "    print(\"Reconstruction is close enough to original matrix.\")\n",
    "else:\n",
    "    print(\"Uh oh - reconstruction is not close enough to original matrix.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise - Pseudoinverse\n",
    "\n",
    "Non-square matrices don't have true inverses, but the SVD can be used to find a pseudoinverse (also called the [Moore-Penrose Inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse)) of any matrix. This has important applications in solving least-squares problems, amongst other things.\n",
    "\n",
    "For a matrix `A`, the pseudoinverse `pinv(A)` doesn't behave exactly like a true inverse, but does satisfy some conditions:\n",
    "\n",
    "1. `A @ pinv(A) @ A == A`\n",
    "2. `pinv(A) @ A @ pinv(A) == pinv(A)`\n",
    "3. `pinv(A) @ A == I` *only* if the columns of A are linearly independent\n",
    "4. `A @ pinv(A) == I` *only* if the rows of A are linearly independent\n",
    "\n",
    "Look up the formula to compute the pseudoinverse of the matrix here, and check these conditions hold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "matrix = np.array([\n",
    "    [0, 1],\n",
    "    [1, 2],\n",
    "    [2, 3]\n",
    "])\n",
    "\n",
    "u, s, vt = np.linalg.svd(matrix)\n",
    "\n",
    "# Construct the matrix sigma\n",
    "sigma = np.zeros((u.shape[1], s.shape[0]))\n",
    "for i in range(s.shape[0]):\n",
    "    sigma[i, i] = s[i]\n",
    "\n",
    "# Find the pseudoinverse of matrix\n",
    "# Your code here...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test the 4 conditions\n",
    "# Your code here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Reading\n",
    "\n",
    "You will cover more on linear algebra in your computer graphics and geometry courses, and SVD in particular will be very important in machine learning. For more information on numpy's linear algebra functionality, check the [documentation here](https://numpy.org/doc/stable/reference/routines.linalg.html).\n",
    "\n",
    "If you're looking for more resources on Linear Algebra (or anything else) during your degree, do check the [UCL Library](https://www.ucl.ac.uk/library). It includes lectures and other resources as well as books, and many resources are available online. If you only want to see online resources, there's a check box for that in the search engine.\n",
    "\n",
    "If you find that numpy does not have a function you require, scipy is a more complete library with many more linear algebra functions (see [documentation here](https://docs.scipy.org/doc/scipy/reference/linalg.html))."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.5 ('base')",
   "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.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c921a71fa7d5301530abc76a09084737c3fca4b7328b11bbb5dba4e94f78eaed"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}