pyqm-week6/notebooks/pyqm-6.1-scipy-linear-algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc8fd83d",
   "metadata": {},
   "source": [
    "<center>\n",
    "    \n",
    "# Python For Quantum Mechanics\n",
    "# Week 6: Linear Algebra\n",
    "\n",
    "</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "4bcff732",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEAAUDBAcHBgcHBwcHBwcHBgYHBwgIBwcHCAcHBgcHBgYHBwcHChALBwgOCQYGDRUNDhERExMTBwsWGBYSJBASExIBBQUFCAcIDQcIDRIICAgSEhISEhIeEhISEhISExISEhISHhISEhISEhISEhISEhISHhISEhITEhISEh4SEhYeEv/AABEIAWgB4AMBIgACEQEDEQH/xAAcAAADAQADAQEAAAAAAAAAAAAAAQIHBAYIBQP/xABGEAEBAAIBAgEHBwoFAgMJAAABAAIDBAUREhMXIWWkpeQGFBgxZLPjCCIkQWZxc3SBsgcVIzJhNFEzkbQlQkNSVGJydqH/xAAaAQEBAQEBAQEAAAAAAAAAAAAAAQIEAwUG/8QAMREBAAEBBQcEAQMEAwAAAAAAAAERAgNSsdEEFBUhMVGhEhNBcWGBkfAiQsHhBTKi/9oADAMBAAIRAxEAPwDxlERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBF6dPyTf2k9z/HR9E79pPc/x10bredstXl71ju8xRenH8k/9o/c/x0n8lD9o/c/x03W97ZanvWO7zJF6ZfyUv2j90fHS+in+0fuj46u6XvbLU9+x3zeZ4vS/0VP2i90fHSfyVv2i90fGzc73Dlqnv2O+bzTF6V+iv+0Xuj42X0WP2i90/Gzc73Dlqe/Y75vNcXpN/JZ/aH3T8bT9Fv8AaH3T8bXcr7DlqbxY75vN0XpD6Ln7Q+6fjZfRd/aH3T8bNyvsOWpvFjvm84RejvovftB7p+Nl9GD9oPdXxtdxvsOWpvF33zecovRj+TD+0Hur42X0YvX/ALq+Mm432HzGqbxd983nSL0U/kyev/dXxlP0ZfX/ALq+Mm4X+HzGpvN33zed4vQ7+TN6/wDdXxkvoz+v/dfxldwv8PmNTebvvm88xehfo0+vvdfxkvo1evvdfxk3C/w+Y1N5u++bz3F6D+jX6+91/GSfybPX3uv4ybhf4fMam83ffN59i9A/Rt9e+6/jJfRu9e+6/jJuF/h8xqbzd983n+L0B9G/177s+Ll9G/177s+LnD7/AA+Y1N5u++bAIt++jh6992fFyfycvXvuz4ucPv8AD5jU3m775sCi3z6Ofrz3Z8XJ/J09ee7Pi5w+/wAPmNU3q7756MEi3r6Ovrz3Z8XL6O3rz3b8XOH3+HzGpvV33z0YNFvH0d/Xfu34uX0ePXfu34uvD7/D5jU3q7756MIi3f6PPrv3b8XL6PXrv3b8XOH3+HzGpvV33z0YTFuv0e/Xfu34uX0fPXXu34qcPv8AD5jU3q7756MLi3P6Pvrr3d8VL6P3rr3d8VOHbRh8xqb1d989GGxbj9H/ANde7vio8wHrr3d8VOHbRh8xqb1d989GHRbh5gfXPu74qXmC9c+7vipw7aMPmNTerrvnoxCLbvMH6593fFS8wnrn3d8VOHbRh8xqb1dd89GJRbZ5hfXPu/4qXmG9ce7/AIqcO2jD5jU3q6756MUi2rzD+uPd/wAVHmI9ce7/AImcO2jD5jU3u6756MVi2nzE+uPd/wATLzFeuPd/xNeHbRh8xqb3dd89GLxbP5i/W/u/4mXmM9b+wfEzh20YfMam93XfPRjMWy+Y31v7B8Tftwf8B/Lb9On/ADjw+V3atXi/y/v4fK544eLt869PbxfV/wAUn/jtojn6fMam9XXfPR7N7+j+knKjv6JLfUiHGblStK0uVqiKcqFktC1oilpWlaVtURS0OUloWtEUtK0rStqgpaFktC1oinKhZLQ5WqJVS0uVK0raiEUtK0rStaJVS0LJaVqilpWlaVgpaVktPeBrJZLStQ1pWS0rVKmtKyWlYGslpWlYKWlZLStaBrSslpWClpWS0rA1pWS0rWgpaVktK1SprSslpWIpaVktKwNaVkslga0rJaVrRDWlZLSsRS0rJaVga0rJZLUNbldDf07ifznF++wuEtyuiP6dxP5zi/fYUtRylY6th7+gpWnv6P6SW+bEOs1pWlaVrRDWlaXKlbVBS0rStK2qIa0rJaHKtEUtDlJaVtUSprS5UrStaIpaVpWlbVEqa0rJaVqiloWS0rBS0rSslgay70rJahrSslpWpU1pWS0rENaVktKwUtK0rJZQNaVktK1DWSyWlYGtKyWla0FLStKyWqGtKyWlYilpWlZLA1pWS0rBS0rSslrRKmtKyWlYilpWS0rA1pWS0rBS0rJaVrQNaVktK1FLcnor+m8T+b4332Fw1uT0V/TeL/N8b77ClrpKx1bB39H9KVpMvR/Slb50Q6qqWlaVpW1EIpaFktC2qJVS0rStK1iENaVpcqVtURTlStK0raoiloWS0rVDWS0rSsFLQslpWClpWSyWtA1pWS0rUNaVktKxDWlZLSsFLStK0rBS0rJaVqGslpWSwNaVktK1oGslpWlaopaVktKwNZLSsliGtKyWlYGslpWSyiVNaVktK1DWS0rSsRS0rJaVgayWlZLWga0rJaVqKWlaVksDWlZLSsDW5XR39M4v83xvvcLhLcno7+mcX+a433uEtdJI6tf7+j+lK09/QfuktwRDpqa0uVLlStqIQ3KlZOVC2qJVTlStK0rWgpaFktK2kNaVktKxDWlZLSsFLSslpWtA1kslpWBrSslpWqVNaVktKwNaVktKwUtKyWlahrSslpWClpWS0rUNZLStKwUtKyWlapU1ktK0rEUtKyWlYGtKyWlYKWlZLStUNaVktKxFLSslpWBrJaVpWopaVktK1DWS0rSsFLSslpWBrJaVktaJU1pWS0rUNbk9If0zjfzXH+9wuGtyOkv6Xxv5rj/e4UnoQ17xeg/cUrSPoP3FLlcUQ6aqWlaVpW1RDWlZLStaJU1pWS0rVDWlZLSsDWSyWlYGslpWS1DWlZLStUqa0rJaVga0rJaVgpaVpWS1oGtKyWlYGslktK1DWlZLSsDWS0rJapU1pWS0rEUtK0rJYGtKyWlYKWlaVktaIa0rJaViKWlaVksDWlZLStRS0rJaVqGtKyWlYKWlZLSsDWlZLStUqpaVktK1DWlZLSsDWlZLSsDW5HSn9L438zx/vcLiLcjpT+lcb+Z0fe40noQ10fQfuKVpH0H7pLcsQ9zWlZLStUNaVktKwNaVkslga095LStRS0rJaVqVNaVktKxDWlZLSsDWS0rSsoKWlZLStQ1ktKyWBrSslpWtA1ktK0rVFLSslpWIayWlaVgpaVktKwNZLSslrRKmtKyWlYhrJaVksDWlZLSsDWS0rJa0DWlZLStRS0rSslga0rJaVrQUtK0rJaoa0rJaVga0rJaVga0rJaVgayWS0rKBrcjpb+lcf+Z0feY3FW/bpj+lcf8AmNH3mMnojXB9H9CS0j6D9xJble5rSslpWClpWS0rUUtKyWlahrSslpWIa0rJaVga0rJaVlBS0rJaVqGtKyWlYKWlZLStQ1pWS0rVFLSslpWFTWlZLSsRS0rJaVga0rJaViVUtKyWlahrSslpWIpaVktKwNaVktK1oKWlZLStQ1pWS0rBS0rJaVqGslpWlYilpWS0rUNaVktKwNaVktKwNZLT3ksD7yWlZLVKmtKyWlYhrft01/SdH8xo+8xuMt+3Tn9J0fx9P3mMmFhro+g/cSWkfQfuJLctHuayWlaVqKWlZLSsSprSslpWBrSslpWBrJaVktQ1vy5O/HVrz2Zvhw14Z7M3svhw14ueb2xO72B+qpb5/wAo3/2fzf5Hl/cZ0tTSK9iOr53TPlr0rlZ6MNXJzMuT4fmzu4nN4mHIcjxYY6NvL04YbskRDFV/VffWzPoXT+p9T6D0PiZcficbha9HRORnyvnee/kZauFhp36zRxzRiad2brwO7m+Eyy/3X3/k91rk8jj9e2bcxy4XV+r8bjdsMMTXp4enVloxQPz0yclcu6+J/wCA57q+maeqOvT4+Kz1ny9bd3EdPh93rnWeLwNJu5e7HThlsx14dzLPPbtz7+DVp1axz3bHs9scRfQ3G6L8ouFztmzVo25+X0hlt0btG/i78MMv9ubo5WGObre5+cCfq73Sfk1zdvP670PdzU2ZHyOw5+hcMccXn8zdq1c3bgYhjjn5LwnYPQZ31/lwGHWvk5v1+jk5c/lcV7KOzh7eLs2crHMP9+GLhqy9P1L3/XIvrUx64/6RMRT550/P56ak3cR/T/dSXZemdX4vLz5OHH3Y7cuHyc+LyQMh1b9YOet8QeLt4vrO53H0+ho5fWuJp5vG4OzcY8vmY7s+Nq8OzJ2Y8fB2bsvFji44Bjjl/uTv4Xt37WdfJbP/AC7r2fK9Jx+udX6707kKnhw5/C6ly8+m5Kvfvnqy3agP1k8X5z8qejdUfq5XM67xeIv/ANB03p27Rqyx/wDt2b8+btH9eO3CzG0zMRy/qm1ET9TMc/2mP1WbqKz2pP79vEtP27McMcs88jDHHFyyyyTHHHHE75ZZZPoxAF7/APF83ovW9HOMs+PjvdIYuG/PRs06d49/zuPltB3Yeg/ODwvc7LfE/wAW89n+R79evt35HJ6fxkyyywxyw5HN0atmGeeIuOGWOTioL2yZ59S6nwedwNPNy4PI43P3Z8XHLjcfdxc+JyMdOzfpxMdm7Yb9Kac8e/5qeh/4va1e0t+n+2KV/WZiP5zYixWzX5mvh2pZLda/xK53M4nR+XyuDv18fbx9WW1yz0G9yxPzfDrMszHXl3yHxZGZ6PqvqZZcz55r7fNvmPzXLynfynzn5148fJ+Dt+Z5HweLv39Pe36/6ps0msU810Y9PKvevhyOoc3TxtOzfv2YadOrFz2bM3w44Yn61/qH/PcvmdH+U3B5u3yOjZt8o63bhhu4vK4rs1DiOzV8514+VwHPDu49+3iO/wBZR8t+i5dT6bu4mvcaNmWWjbq2ZYeUwNnH3Yb9Zswf92C6wf39+z27Xyem/KfqGjncbgdZ4WvTs5XlMePzOJtdvF37dWDs2a8teZ5Tj5OOPc79+/8A29HcxbvJs24ieViac6V59uvL46tWbMTZrHO1/P3dvW4XTeq8flvIOPs8p815W3h7/wAzZh5PkaPC7df5+J4u3jx9OPce/ob4HyO5HO6jo09U28zLVp5Gezbq4OGjjeSx4vjzw0Y7duWDuy3uOOOTljmHd7GPouByPlJzceH1bZjnrdmj5R4dN4znqHHTo27OBoHwYOPlcsfnO3Lvk/Wnf0Hae/FItUmLM1n7ile/I9rnNnrMfzs7wtK3U8eby+D1rgcHfzM+do6np53gdunjatujkcDDHfkmfFwwxy1Za8sjwuKjiem43H69y8vkzzue7D51px6zlrz8nr7Y/NeTydfHPJ+HwpjjrwPSPft6e9r345xMTE2a16fFJ7/mE9qeXaaea6O6LcPgdQ08k2ujPxmjkbuNtfDlj4d+hMd2H5weLsvbuej/AJr4mxz06ssnvllq15ZP1d3LAX0H1elut/4dv+l1X/8AY+t/+ptTb/qsx8WqsxZ5TPZ9/idS0bt3I0atni28TPXhyMfDmeTy26zdrPFkBn3wyH81bkLdK678puVxdXyi24OvN6dt6fhxcM8PzcfnOjjZbHPwdstn5/IyfS/qC5G/lc3g9X6fx9vNz5nG5+vqBsx26ONrz0bOFoOSZ6cuNrx745C4+HLxdux6WxF/HSefP9q2psx894/y37U/z6iZ+Py7YtK2f9F+UXUeTjwuZicvYcvlazdwzp2ZxdHC37XXjt1c3yJlnt14eDNyc3HL847Hov0+WHyg2afn23h8/lZbODiuWjX0z5xwcNmnDHZt0cvlmnuZovdNmPg8R6PQ95vdj0+vnT9O1e/jr+F9i1X066f6d6WlbrHI5nL5fUs+Hp5TwtOrp3E5eWWrTp27tmfK2bcPB4+TjlhhrxNX6sO6v1n1XZMsg+t7B9a/UB9a3vYtxarTpDytWafq4j1XjfPfmPlT5383eT5Ltl38gZmvx+Lt4f8Ac/V37/W9qdfU9eXN28Mxz8pq42nkZZIeBw37NurDHF793Lvpz7+j9Z9dnW/qvHM+N1zx5eXy6vls2+LVvwMekcrE6drxc3DwuJq18bd2Xt4s87s3J5mPG6z1Tk5i48foHE35B9bjo28/bkH/AD2xbmsbT6u1K/8AmYmY/Xk9rVzT9vNYdqWlbpOvqXUtXTeH1PZzMdry9vT8tnFeNqx0a9PUd2rAw0Z4dt2OeOO7H87LLMXF9Hpo6r13nbOR1L5pnycXgbXRxtOnp+XK1crfr04btmPK3mtcDLLPHAMMsEPzle96TtVmI5xMTPxy+Yme/wCPtj2ZnlEw7UdU1vPy4Phz8rjw8OW5dsfB5PZuz4+OI9+/j8WvJ+rt27em5i3TM8N/J6/3w37eDln8neFs2eHVpz24uXN5GTq/SMMscEXs98V9H6qH5Q83Hg8jW56nmcfrWnpXzh1d8M8d+3TjhynRjkHlPJ7jvgId8f6UjaaV9UTSs0/T4+/sm6rSn48u6LSt17qu7Lh8fXq3dT5WW/fyE17MOHo38raY4eLPTo4+jQ49jw+JycMu3d9PpO0/I7qW7ec3Vvzz2vF5Zq17NvH+a7s9WenVvw8tp7BjmeUyO5jiID2vWL6PVFilLU/X38Sz7c+n1dYj70dhWlZLSt7PM1pWSyWBrSslktUqa0rJaViGtKyWlahrSslpWoa37dOf0jR/H0/eY3HW/Xp7+kaf4+n+/GkjXh9B+4ktI+g/cUrcr3UtKyWlYGtKyWlYGtKyWla0FLSslpWBrJaVksDW4/UNBu0btKuJu07dSnpcTbhlrUP1p4r9VpWtPgcPoXAx4XB4vDxzdmPE4vH42OaA548fXjqMsg9Apj3/AK3TvlLwtfS8epbMOsfM9XVM+TyHhZaePu37uZu1GvcdPc3xuexNZ4PBs7Kdu13taMscXIyccXLHuY5djxYmX+4MvrO/Y/8AK87d1FqzSOVOnXl8fEw1Zt0mvWrqfTfku7ekdEx2bNvB6j03gcTHVv1mvLZx9jxterk6M8Mxw26svD4csH0PhPTfh8kOn8zLqW3ldXw5O3n8bXt08bkeHjYdOw42zanfga9WbmbtmBg5uw8R28PcPQ9zWlbPsWYmJw/tNOkz+fye7POO7rXV/khp5PTeTwMt+3Dy/Ud/UMN+AG3j79/Ny5/fV/2cXPPAfr7ZNyc/k5oN/SNuvLLVh0fXydfH1AOOeHI4hwgzy+s8OJ39H1t9paVtezY6058vE1jynuWu/fy4Xyh6Xp6hwt3D3+I17sAcsHw7NeWGWOzVt15dvzdmOeGGQ+n04l83j9A3Zcri8nn87Lm5cHyjxsDj6+Nh5XZrdOXJ344ZPlt/gyzBPDieNTEvvLStqbuzM1nrp0+6JFuYikPh9Q6Ds5PTOZwOTzdu753lyO251a8c9OrdsM9enHDHsZGAdu79f/HoDn5cTZ8717zk7TVhxs9LxQw8lszyyMjfk9vF4wO319v/AO9+WslrF3Zjn88u/wAdM09cuJ1nh5cnQ6sOTv4mfiwzx3cfLDHZjlryMw7bMcsc8Ht2cU7IpfI4nyc2vM4/M5/P3c/ZwzZ80w8jp4unVntx8nnuz1aT/V3eDuCvY8T2D9XYFpWWruzams/H38eCLcxFIdf6N8nt3B7aeNz88OBjvy3a+M8bRnnqxz2u7Pja+Tl9XHXLM7ODkGT2yPQnxflr0g4vS+b5Pdn4ud8oOBzHLw4d9Oe/m8LX2wMhMvD5HFPEP1+ku8LTn2Tsgn/ZO5/3+psWrizNn0xy5TEdeVeXdqL2Yms83xOJ0HZ/mGvn8zlvM28bVu08TE4+vj69ByPCb9rji5OzdljgYuXcO3f826r8r+nnB6dz+Dp6o9ubjy8+J07yGjby9mzm555ZaNGR3zy47tzy7vgXEcvzjt3tEWjw4+LxdjxdvD4ux4vD37+Hv9fbuvo/5pb2eLUUjlM158/nr88/1LN7MTWecR9J4uDhq14P14a9eL+/HExe3/lfC6d0HkcXZyHRz3HVyOfyebnqy4urPtnytvldmBs8Rl4f1f8Ae+8slvabuzNK/wBv3GTEWpj9XXur/JbXytXVdbvzwOq7OJnsTDFdPzTDRrxMO/8Au7nHPr/+Zuf1LpWO/ncHmOeWOXBeY44AOO355p8hkZr6Tsen0X0FpWntWO3X/Ezazmq+uf59UyfC6T0Hdwg0cfn7MODhudmvjvH057NWGWzyuXGw5Of/AMDu5HZwcgy7GR2E4vI+SinO0a+bt1cLqO7dv5PHNWrPM2cnt868hyNnfyWGfb0jjl27vbtdlWS09ixSlOUfmf269Px0X3LXXR8/jdKw187bzMc3vs4fG4vk+3oxx42e3PHLxfWr5VP6X69Y4jyeLv4/lMtXl9Wep2YdnPDHYeDNw7/Vl4XLs/q73JWlb0i7inp+J/yx6prX5hxep9P1cjhbeFkeHTt4+fH7Y/8AuYZYeTx8J+px9Cf/AIl8zT0Bw369ryXaPTtXT+Zhs1Y5HL1aTZ4M++OQ6dq7su7+cIp2PrvtrJaWrqzM1mOcEW5iKOg9Y4Lp18TpWrqTyXVzuneR4Rq0O/Xx9PI17nLl7MO+Rp16sFMvDh/tx7r9T2Tb0bbhyt/I4nLeN87ywz5OvLRhyMMtmGBrN2nxZHktjjjgP+4fCei+t2BUAcu3iQBy7HY7v6/QELeVnZ7MTMz+KdeVK/mvz9fhub2Z/kc/D5HL6Pm8jHl6OVnq5Rw8eJnns06t2G7Vhm7cctmo8HbPx5ZvfBxPzu3b6ri6vkxrx4mPHd+zPN6jr6jv35Y4ePkcnDdjvzyyxx7Y4C4Y49j6jEvvrSt6Tc2J5zHX7+erPuWv2fM670n5zt4vI17suPyeHnty07DDHbj4d+Br369urJDZhliY/rE8J2auj9N+bPIzy257t3K3+W3bMscMByNeGnDHXrwO2GBhrwOz3fR6VuetK19uz6vVTn/qn105J65pT4NZLJaVvSjBrJZLStUNaVktK1DWlZLSsDWS0rJahrSslpWClv14D+kaf4+r+/G4y368F/19P8bV/fjEa8PoP3ElpH0H7iS3I6DWlZLSsoKWlaVktQ1pWS0rBS0rJaVrQNaVktK1RS0rJaViGtKyWlYKWlZLSsDWlZLStaIpaVktKxDWS0rSsFLSslpWBrJaVpWtBS0rJaVqGtKyWlYKWlZLStaFTWS0rJYhrSslpWoa0rJaVga0rJaVgayWlZLKIayWlaVqVUtKyWhYilpWS0rUNaVkslqGtKyWlYKWlaVpWIpaVpWSwNb9OE/6+r+Nq/vxvwW/Thv+tq/i6v78ajXsX0H7iFoH0H7iFuR0GtKyWlYKWlaVktQ1pWS0rVFLStKyWFTWlZLSsQ1ktKyWBrSslpWUSqlpWlZLVDWlZLSsFLStKyWBrSslpWtBS0rJaVqGtKyWlYKWlZLStQ1pWS0rEUtKyWlahrSslpWBrSslpWBrSslksDWlZLJapU1pWlaViKWlZLStaBrSslpWoayWS0rA1ktK0rEqpaVpWlYKWlZLStQ1pWS0rBS18R/1tX8XX/eX4rXxH/V1/wAXX/cQa+PoP3ELQPoP3ELcroNaVktK1DWS0rJYlTWlZLSsQ1ktKyWBrSslpWJU1ktKyWoa0rJaViGslpWSwNaVktK1oGslpWS1DWlZLSsDWS0rJYGtKyWlaoayWlZLUNaVpWSwNaVktKwNaVktKwUtKyWlapVS0rSsliGtK0rJa0DWlZLStQ1kslpWBrJaVpWJVS0rStKwUtK0rJahrSslpWBrJaVksDWlZLKBrXxf/F1/xNf9xfnfpxv/ABNf8TD+4qjXB9B+4ktI+g/dJbldFVLSslpWIa0rJaVgpaVktKwNaVktK1RS0rJaViGtKyWlYKWlZLStQ1pWSyWoa0rJaVga0rJaVgpaVktK1KmslpWlaopaVktCwUtKyWlYGtKyWlYGslktKyiVNZLStK1FLStK0rVFLSslpWoa0rJZLA1pWS0rBS0rStKxFLStKyWoa0rJaVgayWlZLA1pWS094GsoiAiIgL9OKf6uv+Jr/uL879uCd9+n+Nq/vxg1Tv6JLLvStzPY1ktKyWBrSslpWtEUtK0rJYhrSslpWClpWlZLUNaVktKygpaVktK1DWlZLSsFLSslpWtCprSslpWqKWlaVpWClpWlZLA1pWS0rA1pWSyWUDWS0rStUUtK0rSsRS0rStK1FLSslpWoayWS0rA1ktK0rEUtK0rSsFLSslpWoa0rJaVgpaVktKwNZLKICIiAiIgIiIC5HTDvydH8xo+8xuPcrpB35fGPtPH+9wkjS1pWS0rc1HqayWlZLVDWlZLSsDWS0rJYGtKyWla0DWS0rJahrSslpWBrJaVktSprSslpWIayWlaVqKWlaVksDWlaVksDWlZLJYH3pWS0rVKqWlaVpWIpaVpWS1oGtKyWlahrLvJaVgayWlaViVUtK0rSsFLStKyWoa0rJaVgpaVpWSwNaVl3iAiIgIiICIiAiIgIiIC5fRP+t4n83xvvsLiXM6H/ANdxP5zi/fYUnoQ0NaVktK3g9TWS0rSsRS0rJaVrQNZLStK1FLSslp7wNaVkslga0rJaVqhrJaVpWopaVpWlYKWlaVpWClpWS0rA1kslpWqGslpWlYilpWlaVrQUtK0rJahrSslksDWlZLSsSqlpWlaVgpaVpWlailpWS0rA1ktKyWBrSsllAd4iICIiAiIgIiICIiAiIgIiIC5nQv8AruH/ADnF++wuHczoP/X8P+d4n3+ulroQ76tK2NP+OPqj2/4aXnw9U+3/AA18viOz4vE6O3dL3tlq2RaVsd89/qn2/wCGl57vVPt/w1eI7Pi8Tobpe9stWxLStj3nt9U+3/DS89nqn274acS2fF4nQ3S97ZatgWlbIPPX6q9u+Gl56/VXt3w84ls+LxOhul72y1a+tK2Q+er1V7d8PLz0+qvbvh68S2fF4taJul72y1a6tK2Seef1V7d8PLzzeq/bvh5xLZ8XidDdL3tlq1taVslf8ZfVft3w8vPJ6r9t+HnEtnxeJ0N0ve2WrWVpWyfzx+q/bfh5P+MXqz234ecS2fF4nQ3S97ZatXWlbKfPD6s9t+Hl54PVntvw84ls+Lxa0N0ve2WrVVl3sq87/qz238CPO96s9t/ArxLZsfi1om6XuHLVqayWyzzu+rPbPwJedz1b7Z+BOJ7Ni8WtDc73Dlq1JaVsu87fq32z8CXna9W+2fgV4ns2Pxa0Tc73Dlq1BaVsw87Pq32z8CXnY9Xe2fgTiezY/FrQ3O9w5atOWlbMvOv6u9s/Al51vV3tn4FeJ7Nj8WtDc73Dlq01aVsz86vq72v8CXnU9Xe1/gTiezY/FrQ3O9w5atLWS2a+dP1d7X+BJ/xS9X+1/gziezY/FrQ3O+w5atJWlbN/Oj6v9r/Bl50PV/tf4M4ns2Pxa0Tc77Dlq0daVs585/q/2r8GXnO+we1fg14ps2Pxa0NzvsOWrRVpWzvzm/YPavwZecz7B7V+DOKbNj8WtDc77Dlq0RaVs885f2D2r8GXnK+we1fgzimzY/FrQ3O+w5atCWVn3nJ+we1fgx5yfsHtX4M4ps2Pxa0NzvsOWrQYs+85P2D2r8GPOT9g9q/BnFNmx+LWhud9hy1aDFn3nJ+we1fgx5yfsHtX4M4ps2Pxa0NzvsOWrQYs+85P2D2r8GPOT9g9q/BnFNmx+LWhud9hy1aDFn3nJ+we1fgx5yfsHtX4M4ps2Pxa0NzvsOWrQYs+85P2D2r8GPOT9g9q/BnFNmx+LWhud9hy1aDFn3nJ+we1fgx5yfsHtX4M4ps2Pxa0NzvsOWrQYs+85P2D2r8GPOT9g9q/BnFNmx+LWhud9hy1aDFn3nJ+we1fgx5yfsHtX4M4ps2Pxa0NzvsOWrQbm9A/6/h/zvE+/wBdmPnJ+we1fg379P8A8UfI8jTu/wAv8Xkd+nb4fnfbxeSzx2eHv5F7d/D27/8ANJ/5TZqf9/FrRY2O9w5as3iIvyL7giIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiIP/Z",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"700\"\n",
       "            height=\"400\"\n",
       "            src=\"https://www.youtube.com/embed/EQWi8rliCqY\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x15fdb2610>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('EQWi8rliCqY', width=700, height=400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "c6b583ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "\n",
    "import scipy\n",
    "import scipy.linalg as linalg\n",
    "from scipy import stats\n",
    "# Here we imported scipy as well as linalg from scipy. We should normally avoid such redundancy.\n",
    "# SciPy's recommended way of importing is following -> from scipy import linalg, stats, ...\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.linalg is linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ab35b4f",
   "metadata": {},
   "source": [
    "You may be wondering why we are diving into another linear algebra library considering we have already studied `numpy.linalg`. Simply put, `scipy.linalg` contains everything that `numpy.linalg` does but offers additional functionality. It is used in conjuction with numpy ndarrays. \n",
    "\n",
    "> Q: Try finding out what are key differences and similarities between numpy.linalg and scipy.linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f6a755e",
   "metadata": {},
   "source": [
    "## Linear Systems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b2a0419",
   "metadata": {},
   "source": [
    "We can use `scipy.linalg.solve()` to solve, for $\\bf{x}$, linear systems of the form\n",
    "$$ A x = b $$\n",
    "where $A$ is a matrix, $x$ and $b$ are vectors.\n",
    "\n",
    "> Q: Mathematically, will you always find a solution of $A x = b$? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "93e64a95",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [[1 0 2]\n",
      " [0 1 2]\n",
      " [2 1 0]] \n",
      "\n",
      "b =  [[ 6]\n",
      " [12]\n",
      " [18]] \n",
      "\n",
      "x =  [[ 4.]\n",
      " [10.]\n",
      " [ 1.]]\n",
      "Ax - b =  [[0.]\n",
      " [0.]\n",
      " [0.]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,0,2],[0,1,2],[2,1,0]])\n",
    "b = np.array([[6],[12],[18]])\n",
    "\n",
    "x = linalg.solve(A, b)\n",
    "\n",
    "print('A = ',A,'\\n')\n",
    "print('b = ',b,'\\n')\n",
    "print('x = ', x)\n",
    "print('Ax - b = ', A@x - b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider following code, check out the output of this -\n",
    "\n",
    "``` python\n",
    "aa = np.array([[1,1], [1, 1]])\n",
    "bb = np.array([1,2])\n",
    "print(aa, bb)\n",
    "x = linalg.solve(aa, bb)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24d3844d",
   "metadata": {},
   "source": [
    "We can also solve systems of the form \n",
    "\n",
    "$$ A X = B $$\n",
    "\n",
    "where $A$, $X$, and $B$ are all matrices. Typically, this is often simple generalization of previous $Ax = b$ problem, for a solving equation of given $A$ but for many copies of $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "195e4a48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [[1 0 2]\n",
      " [0 1 2]\n",
      " [2 1 0]] \n",
      "\n",
      "b =  [[6, 3, -6], [12, 6, -12], [18, 36, -18]] \n",
      "\n",
      "x =  [[  4.  11.  -4.]\n",
      " [ 10.  14. -10.]\n",
      " [  1.  -4.  -1.]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAJSCAYAAAAS48v9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzPklEQVR4nO3df3AUdZ7/8dcYYBI1GY0hyWRNQtxVFonrYuJCOH5KbTTccqKshecWwp2ySwl4kKLcDe5+Dd4XU7os35SHguwqiKx71FXWXwWF5FuaoEWiwCaHtSKFddFk+WZEWJwIJxMS+vsHMudsfpDuSX9mMnk+qrrK7umefndvr3n7+vR0eyzLsgQAAADXXRbrAgAAAIYLGi8AAABDaLwAAAAMofECAAAwhMYLAADAEBovAAAAQ2i8AAAADKHxAgAAMITGCwAAwBAaLwAAAENovAAAQMKpqqrSrbfeqtTUVGVmZmru3Lk6cuRIxDqWZamyslI5OTlKSUnRjBkz9Oc//9nVumi8AABAwqmvr9fSpUvV2Nio2tpadXV1qbS0VGfOnAmv89RTT2n9+vXasGGD9u/fr+zsbP3whz/Ul19+6VpdHl6SDQAAEt3nn3+uzMxM1dfXa9q0abIsSzk5OVqxYoV+/vOfS5JCoZCysrL05JNP6mc/+5krdZB4AQCAISEUCqmjoyNiCoVCA9o2GAxKktLT0yVJLS0tCgQCKi0tDa/j9Xo1ffp07du3b/CL/9oI174ZAAAktLNnz6qzs9PY/p566imtXbs2Ytljjz2mysrKfrezLEvl5eWaMmWKCgsLJUmBQECSlJWVFbFuVlaWPv3008Er+m/QeAEAANvOnj2r0SkpOm1wn9nZ2frss8+UnJwcXub1ei+53bJly3To0CG9++67PT7zeDwR85Zl9Vg2mGi8AACAbZ2dnTotaaWkS7c+0QtJ+j+BgJKTk5WWljbg7ZYvX67XX39de/fu1bXXXhtenp2dLelC8uX3+8PLjx8/3iMFG0zc4wUAABxLkXS5gSnFZl2WZWnZsmX64x//qLfeeksFBQURnxcUFCg7O1u1tbXhZZ2dnaqvr9fkyZNt7m3gSLwAAEDCWbp0qV5++WW99tprSk1NDd/T5fP5lJKSIo/HoxUrVuiJJ57Q9ddfr+uvv15PPPGELr/8ct13332u1UXjBQAAHLtMZobP7O5j48aNkqQZM2ZELN+yZYsWLVokSXrkkUf01Vdf6aGHHtKpU6c0ceJE7dmzR6mpqdEX3Aee4wUAAGzr6OiQz+fTLyUlX3Lt6J2V9L914bEQdu7xijckXgAAwLF4TbziVaIcBwAAQNwj8QIAAI6ReNmTKMcBAAAQ90i8AACAYyRe9iTKcQAAAMQ9Gi8AAABDGGoEAACOMdRoT6IcBwAAQNwj8QIAAI6ReNmTKMcBAAAQ90i8AACAYyRe9iTKcQAAAMQ9Ei8AAOCYR2ZSHI+BfZhA4gUAAGAIiRcAAHDMIzNpFIkXAAAAbKHxAgAAMIShRgAA4FjS15OJ/SQCEi8AAABDSLwAAIBjPEDVnkQ5DgAAgLhH4gUAABwj8bInUY4DAAAg7pF4AQAAx0i87EmU4wAAAIh7JF4AAMAxEi97EuU4AAAA4h6JFwAAcIzEy55EOQ4AAIC4R+MFAABgCEONAADAMYYa7UmU4wAAAIh7JF4AAMAxz9eTif0kAhIvAAAAQ0i8AACAY5dJSjK0n0SQKMcBAAAQ90i8AACAY/yq0Z5EOQ4AAIC4R+IFAAAcI/GyJ1GOAwAAIO7ReAEAABjCUCMAAHCMoUZ7EuU4AAAA4h6JFwAAcIzEyx5Xj+PUqVNasGCBfD6ffD6fFixYoC+++KLfbRYtWiSPxxMxTZo0yc0yAQAAjHA18brvvvv0l7/8Rbt375Yk/fSnP9WCBQv0xhtv9LvdHXfcoS1btoTnR40a5WaZAADAIRIve1xrvA4fPqzdu3ersbFREydOlCT99re/VUlJiY4cOaKxY8f2ua3X61V2drZbpQEAAMSEa41XQ0ODfD5fuOmSpEmTJsnn82nfvn39Nl51dXXKzMzUVVddpenTp2vt2rXKzMwc0H4ty9KXX34Zdf0AAAxFqamp8ng8xvZH4mWPa41XIBDotVnKzMxUIBDoc7uysjLdc889ys/PV0tLi371q1/ptttu08GDB+X1enusHwqFFAqFwvPBYFB5eXmDcxAAAAwxX3zxhXw+X6zLQB9sN16VlZVas2ZNv+vs379fknrtuC3L6rcTnz9/fvifCwsLVVxcrPz8fO3cuVN33313j/Wrqqp6refdd9t05ZVp/dYJDFTBgimxLgEJxvfBi7EuAQnntKRpxvfq+XoysZ9EYLvxWrZsme69995+1xkzZowOHTqkzz77rMdnn3/+ubKysga8P7/fr/z8fB09erTXzysqKlReXh6e7+joUG5urq68Mk2pqTReGBxpSUmxLgEJJzXWBSBBmRxmhH22G6+MjAxlZGRccr2SkhIFg0G9//77+sEPfiBJeu+99xQMBjV58uQB7+/kyZNqa2uT3+/v9XOv19vrECQAAEC8ce1etXHjxumOO+7Q4sWL1djYqMbGRi1evFg/+tGPIm6s/+53v6tXXnlFknT69GmtWrVKDQ0N+uSTT1RXV6c5c+YoIyNDd911l1ulAgAAh5IMTonA1R8J/P73v9dNN92k0tJSlZaW6nvf+55eeumliHWOHDmiYDAoSUpKStIHH3ygO++8UzfccIMWLlyoG264QQ0NDUpNJZYHAABDm6sPUE1PT9f27dv7XceyrPA/p6Sk6M0333SzJAAAMIg8MvOoh0S5cy1RHosBAAAQ93hJNgAAcIwHqNqTKMcBAAAQ90i8AACAYyRe9iTKcQAAAMQ9Ei8AAOAYiZc9iXIcAAAAcY/ECwAAOEbiZU+iHAcAAEDco/ECAAAwhKFGAADgGEON9iTKcQAAAMQ9Ei8AAOCYR2ZeYM1LsgEAAGALiRcAAHAs6evJxH4SAYkXAACAISReAADAMX7VaE+iHAcAAEDcI/ECAACOeWQmxeFXjQAAALCFxgsAAMAQhhoBAIBj3FxvT6IcBwAAQNwj8QIAAI6ReNmTKMcBAAAQ90i8AACAYyRe9iTKcQAAAETYu3ev5syZo5ycHHk8Hr366qsRny9atEgejydimjRpkqs10XgBAADHLjM42XXmzBndfPPN2rBhQ5/r3HHHHWpvbw9Pu3btcrCngWOoEQAAJKSysjKVlZX1u47X61V2drahiki8AABAFDwGJzfU1dUpMzNTN9xwgxYvXqzjx4+7tKcLSLwAAMCQ0dHRETHv9Xrl9XodfVdZWZnuuece5efnq6WlRb/61a9022236eDBg46/81JovAAAwJCRm5sbMf/YY4+psrLS0XfNnz8//M+FhYUqLi5Wfn6+du7cqbvvvjuaMvtE4wUAABxL+noysR9JamtrU1paWnj5YCZTfr9f+fn5Onr06KB959+i8QIAAENGWlpaROM1mE6ePKm2tjb5/X5Xvl+i8QIAAFGI5weonj59Wh9//HF4vqWlRc3NzUpPT1d6eroqKys1b948+f1+ffLJJ1q9erUyMjJ01113DV7hf4PGCwAAJKQDBw5o5syZ4fny8nJJ0sKFC7Vx40Z98MEH2rZtm7744gv5/X7NnDlTO3bsUGpqqms10XgBAADH4jnxmjFjhizL6vPzN99803lBDvEcLwAAAENIvAAAgGMemUlx3HqAqmkkXgAAAIaQeAEAAMfi+R6veJQoxwEAABD3SLwAAIBjJF72JMpxAAAAxD0aLwAAAEMYagQAAI4x1GhPohwHAABA3CPxAgAAjnlk5uGmPEAVAAAAtpB4AQAAx5K+nkzsJxGQeAEAABhC4gUAABzjV432JMpxAAAAxD0SLwAA4BiJlz2JchwAAABxj8YLAADAEIYaAQCAYww12pMoxwEAABD3SLwAAIBjHplJcXhlEAAAAGwh8QIAAI5xj5c9iXIcAAAAcY/ECwAAOEbiZU+iHAcAAEDcc73xevbZZ1VQUKDk5GQVFRXpnXfe6Xf9+vp6FRUVKTk5Wdddd502bdrkdokAAMChywxOicDV49ixY4dWrFihRx99VE1NTZo6darKysrU2tra6/otLS2aPXu2pk6dqqamJq1evVoPP/ywampq3CwTAADACFcbr/Xr1+uBBx7Qgw8+qHHjxqm6ulq5ubnauHFjr+tv2rRJeXl5qq6u1rhx4/Tggw/qn//5n7Vu3To3ywQAADDCtcars7NTBw8eVGlpacTy0tJS7du3r9dtGhoaeqx/++2368CBAzp37lyv24RCIXV0dERMAADADIYa7XHtOE6cOKHu7m5lZWVFLM/KylIgEOh1m0Ag0Ov6XV1dOnHiRK/bVFVVyefzhafc3NzBOQAAAIBB5noD6fFEPuTfsqweyy61fm/LL6qoqFAwGAxPbW1tUVYMAAAGisTLHtee45WRkaGkpKQe6dbx48d7pFoXZWdn97r+iBEjdM011/S6jdfrldfrHZyiAQAAXORaAzlq1CgVFRWptrY2Ynltba0mT57c6zYlJSU91t+zZ4+Ki4s1cuRIt0oFAAAOkXjZ4+pxlJeX63e/+51eeOEFHT58WCtXrlRra6uWLFki6cIw4f333x9ef8mSJfr0009VXl6uw4cP64UXXtDzzz+vVatWuVkmAACAEa6+Mmj+/Pk6efKkHn/8cbW3t6uwsFC7du1Sfn6+JKm9vT3imV4FBQXatWuXVq5cqWeeeUY5OTl6+umnNW/ePDfLBAAADvHKIHs81sW71xNER0eHfD6fmpuDSk1Ni3U5SBDXzZsQ6xKQYDzN/xHrEpBwvpR0i4LBoNLS3P/7d/Hv7V8lmfhr2yEpXTJ2fG7hJdkAAMAxEi97EuU4AAAA4h6JFwAAcMyjvp+1Oaj7SZA7o0i8AAAADKHxAgAAMIShRgAA4NyIEZKBoUZZltTV5f5+XEbiBQAAYAiJFwAAcI7EyxYSLwAAAENIvAAAgHMmE68EQOIFAABgCIkXAABwjsTLFhIvAAAAQ0i8AACAc0lJ0mUGcpzz593fhwEkXgAAAIbQeAEAABjCUCMAAHBuxAiGGm0g8QIAADCExAsAADhH4mULiRcAAIAhJF4AAMA5Ei9bSLwAAAAMIfECAADOJSVdmNzW3e3+Pgwg8QIAADCExAsAADg3YoSZxMvEi7gNIPECAAAwhMYLAADAEIYaAQCAcww12kLiBQAAYAiJFwAAcI7EyxYSLwAAAENIvAAAgHMkXraQeAEAABhC4gUAAJxLSrqQemFASLwAAAAMoUUFAADOjRhB4mUDiRcAAIAhtKgAAMA5Ei9bSLwAAAAMofECAAAwhGwQAAA4x1CjLSReAAAAhtCiAgAA50w9QNWy3N+HASReAAAAhpB4AQAA50zd40XiBQAAEL/27t2rOXPmKCcnRx6PR6+++mrE55ZlqbKyUjk5OUpJSdGMGTP05z//2dWaaLwAAIBzFxMvE5NNZ86c0c0336wNGzb0+vlTTz2l9evXa8OGDdq/f7+ys7P1wx/+UF9++WW0Z6VPDDUCAICEVFZWprKysl4/syxL1dXVevTRR3X33XdLkl588UVlZWXp5Zdf1s9+9jNXaiLxAgAAzhlOvDo6OiKmUCjkqOyWlhYFAgGVlpaGl3m9Xk2fPl379u0blFPTGxovAAAwZOTm5srn84WnqqoqR98TCAQkSVlZWRHLs7Kywp+5gaFGAAAwZLS1tSktLS087/V6o/o+j8cTMW9ZVo9lg4nGCwAAOGf4cRJpaWkRjZdT2dnZki4kX36/P7z8+PHjPVKwwcRQIwAAGHYKCgqUnZ2t2tra8LLOzk7V19dr8uTJru2XxAsAADhn6pVB58/b3uT06dP6+OOPw/MtLS1qbm5Wenq68vLytGLFCj3xxBO6/vrrdf311+uJJ57Q5Zdfrvvuu28wK49A4wUAABLSgQMHNHPmzPB8eXm5JGnhwoXaunWrHnnkEX311Vd66KGHdOrUKU2cOFF79uxRamqqazXReAEAAOdM3ePlIPGaMWOGrH5eNeTxeFRZWanKysooCrOHe7wAAAAMIfECAADOxXHiFY9IvAAAAAwh8QIAAM6ReNlC4gUAAGAIjRcAAIAhDDUCAADnGGq0xfXE69lnn1VBQYGSk5NVVFSkd955p8916+rq5PF4ekwfffSR22UCAAC4ztUWdceOHVqxYoWeffZZ/d3f/Z2ee+45lZWV6cMPP1ReXl6f2x05ciTiBZijR492s0wAAOCUqVcGdXe7vw8DXE281q9frwceeEAPPvigxo0bp+rqauXm5mrjxo39bpeZmans7OzwlJSU5GaZAAAARrjWeHV2durgwYMqLS2NWF5aWqp9+/b1u+2ECRPk9/s1a9Ysvf32226VCAAAonXxHi8TUwJw7ShOnDih7u5uZWVlRSzPyspSIBDodRu/36/NmzerqKhIoVBIL730kmbNmqW6ujpNmzat121CoZBCoVB4vqOjQ5JUULlQaSNHDtLRYLj7f83NsS4BCWdlrAtAwjkX6wIwAK63jx6PJ2Lesqweyy4aO3asxo4dG54vKSlRW1ub1q1b12fjVVVVpTVr1gxewQAAYOBMpVHc49W/jIwMJSUl9Ui3jh8/3iMF68+kSZN09OjRPj+vqKhQMBgMT21tbY5rBgAAcJNrLeqoUaNUVFSk2tpa3XXXXeHltbW1uvPOOwf8PU1NTfL7/X1+7vV65fV6o6oVAAA4ROJli6tnqry8XAsWLFBxcbFKSkq0efNmtba2asmSJZIupFXHjh3Ttm3bJEnV1dUaM2aMxo8fr87OTm3fvl01NTWqqalxs0wAAAAjXG285s+fr5MnT+rxxx9Xe3u7CgsLtWvXLuXn50uS2tvb1draGl6/s7NTq1at0rFjx5SSkqLx48dr586dmj17tptlAgAAp0w9x6ury/19GOCxLMuKdRGDqaOjQz6fT8G5c/lVIwbN//uP/4h1CUgw39KPYl0CEs45SW8qGAxGPITcLeG/txUVSktOdn9/Z8/KV1Vl7PjcwkuyAQAADEmMp5EBAIDYMHVzfYI8QJXECwAAwJDEaB8BAEBskHjZQuIFAABgSGK0jwAAIDZIvGwh8QIAADAkMdpHAAAQGyRetpB4AQAAGJIY7SMAAIgNU68MSkpyfx8GkHgBAAAYQuMFAABgCEONAADAOW6ut4XECwAAwJDEaB8BAEBskHjZQuIFAABgSGK0jwAAIDZIvGwh8QIAADAkMdpHAAAQGzxA1RYSLwAAAENIvAAAgHPc42ULiRcAAIAhNF4AAACGJEZuBwAAYoOhRltIvAAAAAxJjPYRAADEBomXLSReAAAAhiRG+wgAAGKDxMsWEi8AAABDEqN9BAAAscErg2wh8QIAADCExAsAADjHPV62kHgBAAAYkhjtIwAAiA0SL1tIvAAAAAyh8QIAADAkMXI7AAAQGww12kLiBQAAYEhitI8AACA2eICqLSReAAAAhpB4AQAA57jHyxYSLwAAAEMSo30EAACxQeJlC4kXAACAIYnRPgIAgNgg8bKFxAsAAMAQGi8AAABDEiO3AwAAscFQoy0kXgAAAIYkRvsIAABig1cG2ULiBQAAYAiJFwAAcI57vGwh8QIAADAkMdpHAAAQGyRetpB4AQAAGJIY7SMAAIgNEi9bSLwAAAAMofECAAAwJDFyOwAAEBs8QNUWEi8AAABDSLwAAIBz3FxvC4kXAACAIYnRPgIAgNgg8bKFxAsAAMAQVxuvvXv3as6cOcrJyZHH49Grr756yW3q6+tVVFSk5ORkXXfdddq0aZObJQIAgGhcTLxMTDZUVlbK4/FETNnZ2S6dhIFztfE6c+aMbr75Zm3YsGFA67e0tGj27NmaOnWqmpqatHr1aj388MOqqalxs0wAAJCAxo8fr/b29vD0wQcfxLokd+/xKisrU1lZ2YDX37Rpk/Ly8lRdXS1JGjdunA4cOKB169Zp3rx5LlUJAAAci+N7vEaMGBEXKdc3xdU9Xg0NDSotLY1Ydvvtt+vAgQM6d+5cjKoCAADxoqOjI2IKhUJ9rnv06FHl5OSooKBA9957r/7rv/7LYKW9i6vGKxAIKCsrK2JZVlaWurq6dOLEiV63CYVCPf5HAAAAZpzXZcYmScrNzZXP5wtPVVVVvdY1ceJEbdu2TW+++aZ++9vfKhAIaPLkyTp58qTJ09ND3P020+PxRMxbltXr8ouqqqq0Zs0a1+sCAACx19bWprS0tPC81+vtdb1v3up00003qaSkRN/+9rf14osvqry83PU6+xJXiVd2drYCgUDEsuPHj2vEiBG65ppret2moqJCwWAwPLW1tZkoFQAAxEBaWlrE1Ffj9beuuOIK3XTTTTp69KjLFfYvrhKvkpISvfHGGxHL9uzZo+LiYo0cObLXbbxe74BPOgAAGFxdXRcmE/uJRigU0uHDhzV16tTBKcghVxOv06dPq7m5Wc3NzZIuPC6iublZra2tki6kVffff394/SVLlujTTz9VeXm5Dh8+rBdeeEHPP/+8Vq1a5WaZAAAgwaxatUr19fVqaWnRe++9px//+Mfq6OjQwoULY1qXq4nXgQMHNHPmzPD8xTHVhQsXauvWrWpvbw83YZJUUFCgXbt2aeXKlXrmmWeUk5Ojp59+mkdJAAAQp+I18frLX/6if/zHf9SJEyc0evRoTZo0SY2NjcrPz3enwAFytfGaMWNG+Ob43mzdurXHsunTp+tPf/qTi1UBAIBE9+///u+xLqFXcXWPFwAAGFriNfGKV3H1q0YAAIBERuIFAAAcI/Gyh8QLAADAEBIvAADgWHe3mTSqu9v9fZhA4gUAAGAIjRcAAIAhDDUCAADHuLneHhIvAAAAQ0i8AACAYyRe9pB4AQAAGELiBQAAHCPxsofECwAAwBASLwAA4BgPULWHxAsAAMAQEi8AAOAY93jZQ+IFAABgCI0XAACAIQw1AgAAxxhqtIfECwAAwBASLwAA4BiJlz0kXgAAAIaQeAEAAMdIvOwh8QIAADCExAsAADjGK4PsIfECAAAwhMQLAAA4xj1e9pB4AQAAGELiBQAAHCPxsofECwAAwBAaLwAAAEMYagQAAI4x1GgPiRcAAIAhJF4AAMAxHqBqD4kXAACAISReAADAMe7xsofECwAAwBASLwAA4BiJlz0kXgAAAIaQeAEAAMdIvOwh8QIAADCExgsAAMAQhhoBAIBjDDXaQ+IFAABgCIkXAABwjFcG2UPiBQAAYAiJFwAAcIx7vOwh8QIAADCExAsAADhG4mUPiRcAAIAhJF4AAMAxEi97SLwAAAAMofECAAAwhKFGAADgGA9QtYfECwAAwBASLwAA4Bg319tD4gUAAGAIiRcAAHCMxMseEi8AAABDSLwAAIBjJF72kHgBAAAYQuIFAAAcI/Gyh8QLAADAEBIvAADgGE+ut8fVxGvv3r2aM2eOcnJy5PF49Oqrr/a7fl1dnTweT4/po48+crNMAAAAI1xNvM6cOaObb75Z//RP/6R58+YNeLsjR44oLS0tPD969Gg3ygMAADDK1carrKxMZWVltrfLzMzUVVddNfgFAQCAQdXVJSUlmdlPIojLe7wmTJigs2fP6sYbb9Qvf/lLzZw5s891Q6GQQqFQeL6jo0OS9Ku8F+X1pvW1GWDLr/VxrEtAgpkw4TuxLgEJpru7Q4cO+WJdBi4hrhovv9+vzZs3q6ioSKFQSC+99JJmzZqluro6TZs2rddtqqqqtGbNGsOVAgAAicTLrrhqvMaOHauxY8eG50tKStTW1qZ169b12XhVVFSovLw8PN/R0aHc3FzXawUAALArrhqv3kyaNEnbt2/v83Ov1yuv12uwIgAAcBGJlz1x/wDVpqYm+f3+WJcBAAAQNVcTr9OnT+vjj//npuSWlhY1NzcrPT1deXl5qqio0LFjx7Rt2zZJUnV1tcaMGaPx48ers7NT27dvV01NjWpqatwsEwAAOMQDVO1xtfE6cOBAxC8SL96LtXDhQm3dulXt7e1qbW0Nf97Z2alVq1bp2LFjSklJ0fjx47Vz507Nnj3bzTIBAACMcLXxmjFjhizL6vPzrVu3Rsw/8sgjeuSRR9wsCQAADKKuLukyAzcucY8XAAAAbKHxAgAAMCTuHycBAADiF0ON9pB4AQAAGELiBQAAHCPxsofECwAAwBASLwAA4BiJlz0kXgAAIGE9++yzKigoUHJysoqKivTOO+/EtB4aLwAA4NjFVwa5PTl5ZdCOHTu0YsUKPfroo2pqatLUqVNVVlYW8dYc02i8AABAQlq/fr0eeOABPfjggxo3bpyqq6uVm5urjRs3xqwm7vECAACOdXVJHo+Z/djR2dmpgwcP6he/+EXE8tLSUu3bt28QK7OHxgsAAAwZHR0dEfNer1der7fHeidOnFB3d7eysrIilmdlZSkQCLhaY38YagQAAENGbm6ufD5feKqqqup3fc/fxHGWZfVYZhKJFwAAcMz0UGNbW5vS0tLCy3tLuyQpIyNDSUlJPdKt48eP90jBTCLxAgAAQ0ZaWlrE1FfjNWrUKBUVFam2tjZieW1trSZPnmyi1F6ReAEAAMfi9eZ6SSovL9eCBQtUXFyskpISbd68Wa2trVqyZMngFzhANF4AACAhzZ8/XydPntTjjz+u9vZ2FRYWateuXcrPz49ZTTReAADAsXhOvCTpoYce0kMPPTS4xUSBe7wAAAAMIfECAACOdXebSbycvDIoHpF4AQAAGELiBQAAHHN671W87sdtJF4AAACGkHgBAADHSLzsIfECAAAwhMYLAADAEIYaAQCAYww12kPiBQAAYAiJFwAAcMzUg015gCoAAABsIfECAACOdXVJluX+fki8AAAAYAuJFwAAcIzEyx4SLwAAAENIvAAAgGMkXvaQeAEAABhC4wUAAGAIQ40AAMAxhhrtIfECAAAwhMQLAAA41t1tJvE6f979fZhA4gUAAGAIiRcAAHCsq0u6zECMQ+IFAAAAW0i8AACAYyRe9pB4AQAAGELiBQAAHCPxsofECwAAwBAaLwAAAEMYagQAAI51d5sZBjTxkFYTSLwAAAAMIfECAACOdXVJHo/7+yHxAgAAgC0kXgAAwDESL3tIvAAAAAwh8QIAAI6ReNlD4gUAAGAIiRcAAHCMxMseEi8AAABDSLwAAEAUzsuyTLzBOjHekk3iBQAAYAiNFwAAgCEMNQIAgCh0fz2Z2M/QR+IFAABgiKuNV1VVlW699ValpqYqMzNTc+fO1ZEjRy65XX19vYqKipScnKzrrrtOmzZtcrNMAADgWLfBaehztfGqr6/X0qVL1djYqNraWnV1dam0tFRnzpzpc5uWlhbNnj1bU6dOVVNTk1avXq2HH35YNTU1bpYKAADgOlfv8dq9e3fE/JYtW5SZmamDBw9q2rRpvW6zadMm5eXlqbq6WpI0btw4HThwQOvWrdO8efPcLBcAANjGPV52GL3HKxgMSpLS09P7XKehoUGlpaURy26//XYdOHBA586dc7U+AAAANxn7VaNlWSovL9eUKVNUWFjY53qBQEBZWVkRy7KystTV1aUTJ07I7/dHfBYKhRQKhcLzHR0dg1s4AADox3mZebgpD1C1ZdmyZTp06JD+8Ic/XHJdz9+89Mn6+gVNf7tcunADv8/nC0+5ubmDUzAAAMAgM9J4LV++XK+//rrefvttXXvttf2um52drUAgELHs+PHjGjFihK655poe61dUVCgYDIantra2Qa0dAAD0h1812uHqUKNlWVq+fLleeeUV1dXVqaCg4JLblJSU6I033ohYtmfPHhUXF2vkyJE91vd6vfJ6vYNWMwAAgFtcTbyWLl2q7du36+WXX1ZqaqoCgYACgYC++uqr8DoVFRW6//77w/NLlizRp59+qvLych0+fFgvvPCCnn/+ea1atcrNUgEAAFznauO1ceNGBYNBzZgxQ36/Pzzt2LEjvE57e7taW1vD8wUFBdq1a5fq6ur0/e9/X//6r/+qp59+mkdJAAAQl87LzDBjYtxc7/pQ46Vs3bq1x7Lp06frT3/6kwsVAQAAxA4vyQYAAFHgAap28JJsAAAAQ0i8AABAFEi87CDxAgAAMITECwAARIFXBtlB4gUAAGAIiRcAAIgC93jZQeIFAABgCI0XAACAIQw1AgCAKDDUaAeJFwAAgCEkXgAAIAokXnaQeAEAABhC4gUAAKJwXmbSKB6gCgAAABtIvAAAQBR4ZZAdJF4AAACGkHgBAIAo8KtGO0i8AAAADCHxAgAAUSDxsoPECwAAwBAaLwAAAEMYagQAAFFgqNEOEi8AAABDSLwAAEAUSLzsIPECAAAwhMQLAABEwZKZ1/lYBvbhPhIvAAAw7I0ZM0Yejydi+sUvfjHo+yHxAgAAUUice7wef/xxLV68ODx/5ZVXDvo+aLwAAAAkpaamKjs729V9MNQIAACi0G1wkjo6OiKmUCg0aEfy5JNP6pprrtH3v/99rV27Vp2dnYP23ReReAEAgCEjNzc3Yv6xxx5TZWVl1N/7L//yL7rlllt09dVX6/3331dFRYVaWlr0u9/9Lurv/iYaLwAAMGS0tbUpLS0tPO/1evtct7KyUmvWrOn3+/bv36/i4mKtXLkyvOx73/uerr76av34xz8Op2CDhcYLAABEwezN9WlpaRGNV3+WLVume++9t991xowZ0+vySZMmSZI+/vhjGi8AAIBLycjIUEZGhqNtm5qaJEl+v38wS6LxAgAA0Rj6j5NoaGhQY2OjZs6cKZ/Pp/3792vlypX6h3/4B+Xl5Q3qvmi8AADAsOb1erVjxw6tWbNGoVBI+fn5Wrx4sR555JFB3xeNFwAAiMJ5mXllkHv7uOWWW9TY2Oja938Tz/ECAAAwhMQLAABEYejf42USiRcAAIAhJF4AACAK52UmjTJxH5n7SLwAAAAMofECAAAwhKFGAAAQBW6ut4PECwAAwBASLwAAEIWh/wBVk0i8AAAADCHxAgAAUeAeLztIvAAAAAwh8QIAAFEg8bKDxAsAAMAQEi8AABAFEi87SLwAAAAMIfECAABR4CXZdpB4AQAAGELjBQAAYAhDjQAAIAq8MsgOEi8AAABDSLwAAEAUeJyEHSReAAAAhpB4AQCAKJB42UHiBQAAYAiJFwAAiAKJlx2uJl5VVVW69dZblZqaqszMTM2dO1dHjhzpd5u6ujp5PJ4e00cffeRmqQAAAK5zNfGqr6/X0qVLdeutt6qrq0uPPvqoSktL9eGHH+qKK67od9sjR44oLS0tPD969Gg3SwUAAI6QeNnhauO1e/fuiPktW7YoMzNTBw8e1LRp0/rdNjMzU1dddZWL1QEAAJhl9B6vYDAoSUpPT7/kuhMmTNDZs2d144036pe//KVmzpw5oH1YliVJCoU6nBcK9PBlrAtAgunu5t9RGFwXr6mLfwcRn4w1XpZlqby8XFOmTFFhYWGf6/n9fm3evFlFRUUKhUJ66aWXNGvWLNXV1fWakoVCIYVCofD8J598Ikl67rncQT8GABgshw7FugIkqpMnT8rn8xncI68MssNY47Vs2TIdOnRI7777br/rjR07VmPHjg3Pl5SUqK2tTevWreu18aqqqtKaNWt6LP/www/1rW99K/rCE1hHR4dyc3PV1tYWcT8deuJcDQznaeA4VwPHuRqYY8eO6cYbb9TIkSNjXQr6YaTxWr58uV5//XXt3btX1157re3tJ02apO3bt/f6WUVFhcrLy8PzwWBQeXl5ysnJ4f+gA5SWlsa5GiDO1cBwngaOczVwnKv+XRxiNH+OzsvMje8kXpdkWZaWL1+uV155RXV1dSooKHD0PU1NTfL7/b1+5vV65fV6eyz3eDyO9gUAwFB08e8ef//im6uN19KlS/Xyyy/rtddeU2pqqgKBgCTJ5/MpJSVF0oXE6tixY9q2bZskqbq6WmPGjNH48ePV2dmp7du3q6amRjU1NW6WCgAAHOFxEna42nht3LhRkjRjxoyI5Vu2bNGiRYskSe3t7WptbQ1/1tnZqVWrVunYsWNKSUnR+PHjtXPnTs2ePXtA+/R6vXrsscd6TcEQiXM1cJyrgeE8DRznauA4VwPDeRoaPBa/OwUAADZ1dHR8/evJxyUlG9jjWUn/S8FgcEjf68dLsgEAAAzhJdkAACAK3ONlB4kXAACAITReAAAAhiRE43Xq1CktWLBAPp9PPp9PCxYs0BdffNHvNosWLZLH44mYJk2aZKZgg5599lkVFBQoOTlZRUVFeuedd/pdv76+XkVFRUpOTtZ1112nTZs2Gao09uycq7q6uh7Xj8fj0UcffWSwYvP27t2rOXPmKCcnRx6PR6+++uoltxmu15TdczVcr6mqqirdeuutSk1NVWZmpubOnasjR45ccrvhdl05OU/mrqlug9PQlxCN13333afm5mbt3r1bu3fvVnNzsxYsWHDJ7e644w61t7eHp127dhmo1pwdO3ZoxYoVevTRR9XU1KSpU6eqrKws4vEd39TS0qLZs2dr6tSpampq0urVq/Xwww8Pi2eo2T1XFx05ciTiGrr++usNVRwbZ86c0c0336wNGzYMaP3hfE3ZPVcXDbdrqr6+XkuXLlVjY6Nqa2vV1dWl0tJSnTlzps9thuN15eQ8XTTcrql4N+QfJ3H48GHdeOONamxs1MSJEyVJjY2NKikp0UcffRTx3sdvWrRokb744osB/Rf7UDVx4kTdcsst4eepSdK4ceM0d+5cVVVV9Vj/5z//uV5//XUdPnw4vGzJkiX6z//8TzU0NBipOVbsnqu6ujrNnDlTp06d0lVXXWWw0vjh8Xj0yiuvaO7cuX2uM5yvqW8ayLnimrrg888/V2Zmpurr63t9P6/EdSUN7Dy5fU39z+MkVsvc4ySe4HESsdbQ0CCfzxduuqQL73b0+Xzat29fv9vW1dUpMzNTN9xwgxYvXqzjx4+7Xa4xnZ2dOnjwoEpLSyOWl5aW9nleGhoaeqx/++2368CBAzp37pxrtcaak3N10YQJE+T3+zVr1iy9/fbbbpY5JA3Xayoaw/2aCgaDkqT09PQ+1+G6Gth5umi4X1PxZsg3XoFAQJmZmT2WZ2Zmhl9R1JuysjL9/ve/11tvvaXf/OY32r9/v2677TaFQiE3yzXmxIkT6u7uVlZWVsTyrKysPs9LIBDodf2uri6dOHHCtVpjzcm58vv92rx5s2pqavTHP/5RY8eO1axZs7R3714TJQ8Zw/WacoJr6sL7fcvLyzVlyhQVFhb2ud5wv64Gep7MXVMXX5Lt9sRLsl1VWVmpNWvW9LvO/v37JfX+QlDLsvp9Uej8+fPD/1xYWKji4mLl5+dr586duvvuux1WHX/+9hxc6rz0tn5vyxORnXM1duzYiGHskpIStbW1ad26dX3G/sPVcL6m7OCakpYtW6ZDhw7p3XffveS6w/m6Guh54pqKT3HbeC1btkz33ntvv+uMGTNGhw4d0meffdbjs88//7zHfxH1x+/3Kz8/X0ePHrVdazzKyMhQUlJSj8Tm+PHjfZ6X7OzsXtcfMWKErrnmGtdqjTUn56o3kyZN0vbt2we7vCFtuF5Tg2U4XVPLly/X66+/rr179+raa6/td93hfF3ZOU+9ceea4gGqdsRt45WRkaGMjIxLrldSUqJgMKj3339fP/jBDyRJ7733noLBoCZPnjzg/Z08eVJtbW3y+/2Oa44no0aNUlFRkWpra3XXXXeFl9fW1urOO+/sdZuSkhK98cYbEcv27Nmj4uJijRw50tV6Y8nJuepNU1NTwlw/g2W4XlODZThcU5Zlafny5XrllVdUV1engoKCS24zHK8rJ+epN8Phmop3cdt4DdS4ceN0xx13aPHixXruueckST/96U/1ox/9KCJi/e53v6uqqirdddddOn36tCorKzVv3jz5/X598sknWr16tTIyMiL+8A515eXlWrBggYqLi1VSUqLNmzertbVVS5YskSRVVFTo2LFj2rZtm6QLvwrasGGDysvLtXjxYjU0NOj555/XH/7wh1gehhF2z1V1dbXGjBmj8ePHq7OzU9u3b1dNTU1C/5xdkk6fPq2PP/44PN/S0qLm5malp6crLy+Pa+ob7J6r4XpNLV26VC+//LJee+01paamhpMsn8+nlJQUSfy7SnJ2nsxdUyRetlgJ4OTJk9ZPfvITKzU11UpNTbV+8pOfWKdOnYpYR5K1ZcsWy7Is67//+7+t0tJSa/To0dbIkSOtvLw8a+HChVZra6v54l32zDPPWPn5+daoUaOsW265xaqvrw9/tnDhQmv69OkR69fV1VkTJkywRo0aZY0ZM8bauHGj4Ypjx865evLJJ61vf/vbVnJysnX11VdbU6ZMsXbu3BmDqs16++23LUk9poULF1qWxTX1TXbP1XC9pno7R9/897VlcV1ZlrPz5PY1FQwGv66j3JIqDEzlliQrGAwO2jHEwpB/jhcAADDvf57jVS7Ja2CPIUnrh/xzvIb8UCMAAIil8zLzqIfEeJzEkH+OFwAAwFBB4gUAAKLAzfV2kHgBAAAYQuIFAACiQOJlB4kXAACAISReAAAgCiRedpB4AQAAGELiBQAAonBeZtIonuMFAAAAG0i8AABAFHhyvR0kXgAAAIbQeAEAABjCUCMAAIgCj5Owg8QLAADAEBIvAAAQBRIvO0i8AAAADCHxAgAAUSDxsoPECwAAwBASLwAAEAUSLztIvAAAAAwh8QIAAFHgJdl2kHgBAAAYQuMFAABgCEONAAAgCudlZhiQoUYAAADYQOIFAACi0C0zOQ6PkwAAAIANJF4AACAKJF52kHgBAAAYQuMFAACi0G1wcs/atWs1efJkXX755brqqqt6Xae1tVVz5szRFVdcoYyMDD388MPq7Oy0tR+GGgEAwLDX2dmpe+65RyUlJXr++ed7fN7d3a2///u/1+jRo/Xuu+/q5MmTWrhwoSzL0r/9278NeD80XgAAIAqJ8RyvNWvWSJK2bt3a6+d79uzRhx9+qLa2NuXk5EiSfvOb32jRokVau3at0tLSBrQfhhoBAAAuoaGhQYWFheGmS5Juv/12hUIhHTx4cMDfQ+IFAACi0GV0Px0dHRFLvV6vvF6v63sPBALKysqKWHb11Vdr1KhRCgQCA/4eGi8AAGDbqFGjlJ2drUDg/xrb55VXXqnc3NyIZY899pgqKyt7Xb+ysjI8hNiX/fv3q7i4eED793g8PZZZltXr8r7QeAEAANuSk5PV0tJi+1d90eityekv7Vq2bJnuvffefr9zzJgxA9p3dna23nvvvYhlp06d0rlz53okYf2h8QIAAI4kJycrOTk51mX0KSMjQxkZGYPyXSUlJVq7dq3a29vl9/slXbjh3uv1qqioaMDfQ+MFAACGvdbWVv31r39Va2ururu71dzcLEn6zne+oyuvvFKlpaW68cYbtWDBAv3617/WX//6V61atUqLFy8e8C8aJcljWZbl0jEAAAAMCYsWLdKLL77YY/nbb7+tGTNmSLrQnD300EN66623lJKSovvuu0/r1q2zdXM/jRcAAIAhPMcLAADAEBovAAAAQ2i8AAAADKHxAgAAMITGCwAAwBAaLwAAAENovAAAAAyh8QIAADCExgsAAMAQGi8AAABDaLwAAAAMofECAAAw5P8D33/vr/wnYLYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 700x700 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[1,0,2],[0,1,2],[2,1,0]])\n",
    "B = [[6,3,-6],[12,6,-12],[18,36,-18]]\n",
    "\n",
    "X = linalg.solve(A, B)\n",
    "\n",
    "print('A = ',A,'\\n')\n",
    "print('b = ',B,'\\n')\n",
    "print('x = ', X)\n",
    "\n",
    "fig = plt.figure(figsize = (7,7))\n",
    "ax = fig.add_axes([.1,.1,.8,.8])\n",
    "\n",
    "im = ax.imshow(X, cmap='seismic')\n",
    "fig.colorbar(im, ax=ax)\n",
    "im.set_clim(vmin=-10, vmax=20)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df405a3",
   "metadata": {},
   "source": [
    "## Matrix Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd85aed4",
   "metadata": {},
   "source": [
    "Scipy provides several functions to compute functions of matrices, with either matrix, or scalar output. Common examples include computing determinant, trace, transpose and inverse of matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `scipy.linalg.inv(A)` computes the inverse of a matrix `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "de7afdb5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [[0.79565993 0.93963841 0.80108429]\n",
      " [0.91848384 0.56806837 0.47002508]\n",
      " [0.20451638 0.21208673 0.79074078]] \n",
      "\n",
      "A^-1 =  [[-1.39220359  2.28288218  0.05344425]\n",
      " [ 2.51010537 -1.85353984 -1.44117489]\n",
      " [-0.31316408 -0.09329933  1.6373556 ]]\n"
     ]
    }
   ],
   "source": [
    "A = rnd.random((3,3))\n",
    "print('A = ', A,'\\n')\n",
    "\n",
    "print('A^-1 = ',linalg.inv(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `scipy.linalg.expm` computes the determinant of a matrix `A`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exp(A) =  [[3.48308471 2.44793675 2.62414311]\n",
      " [2.31375521 2.8641687  1.95757257]\n",
      " [0.76268791 0.72523388 2.60711291]]\n"
     ]
    }
   ],
   "source": [
    "print('exp(A) = ', linalg.expm(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd5853e",
   "metadata": {},
   "source": [
    "The function `scipy.linalg.det(A)` computes the determinant of a matrix `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "7ca79cdb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "det(A) =  -0.2510471512872222\n"
     ]
    }
   ],
   "source": [
    "print('det(A) = ',linalg.det(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "527fc382",
   "metadata": {},
   "source": [
    "Another function, `scipy.linalg.norm(A, ord=2)` computes the modulus of order `ord`. If `ord` $=$ `np.inf`, then the maximum value in $A$ is given. If `ord` $=$ `-np.inf`, then the minimum value in `A` is given."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "07743c27",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm(A,ord=1) = 4 + 3 + 1 + 2 + 2 + 1 + 1 =  14.0\n",
      "norm(A,ord=2) = (4**2 + 3**2 + 1**2 + 2**2 + 2**2 + 1**2 + 1**2)**(1/2) =  6.0\n",
      "norm(A,ord=3) = (4**3 + 3**3 + 1**3 + 2**3 + 2**3 + 1**3 + 1**3)**(1/3) =  4.791419857062784\n",
      "norm(A,ord=-np.inf) =  1.0\n",
      "norm(A,ord=np.inf) =  4.0\n"
     ]
    }
   ],
   "source": [
    "A = np.array([4,3,1,2,2,1,1])\n",
    "\n",
    "print('norm(A,ord=1) = 4 + 3 + 1 + 2 + 2 + 1 + 1 = ',linalg.norm(A, ord=1))\n",
    "print('norm(A,ord=2) = (4**2 + 3**2 + 1**2 + 2**2 + 2**2 + 1**2 + 1**2)**(1/2) = ',linalg.norm(A))\n",
    "print('norm(A,ord=3) = (4**3 + 3**3 + 1**3 + 2**3 + 2**3 + 1**3 + 1**3)**(1/3) = ',linalg.norm(A, ord=3))\n",
    "print('norm(A,ord=-np.inf) = ',linalg.norm(A, ord=-np.inf))\n",
    "print('norm(A,ord=np.inf) = ',linalg.norm(A, ord=np.inf))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c6119f",
   "metadata": {},
   "source": [
    "## Spectral Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05722290",
   "metadata": {},
   "source": [
    "We can represent a matrix in terms of it's eigenvectors and eigenvalues, or in other words, perform a spectral decomposition; `scipy.linalg.eig(A)` can be used to achive this by computing the eigenvalues and eigenvectors of `A`, giving a tuple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "802e92ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues =  [ 2.1571515+0.j -0.2511262+0.j  0.0255913+0.j] \n",
      "\n",
      "Eigenvector Matrix =  [[-0.48012849 -0.63035729 -0.46959699]\n",
      " [-0.67539884  0.68623736 -0.16116989]\n",
      " [-0.55974373 -0.36294348  0.86804547]] \n",
      "\n",
      "First Eigenvector =  [-0.48012849 -0.67539884 -0.55974373] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = rnd.random((3,3))\n",
    "\n",
    "eigval, eigvec = linalg.eig(A)\n",
    "\n",
    "print('Eigenvalues = ',eigval, '\\n')\n",
    "print('Eigenvector Matrix = ',eigvec, '\\n')\n",
    "\n",
    "print('First Eigenvector = ',eigvec[:,0], '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbc6fff9",
   "metadata": {},
   "source": [
    "## Special Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ba90699",
   "metadata": {},
   "source": [
    "Scipy supports the creation of special matrices that are commonly used. A full list can be found at https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html. \n",
    "\n",
    "For example, we can create a Hadamard matrix of order `n` using `scipy.linalg.hadamard(n)`. The integer `n` must be a power of $2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "af3f7195",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1  1]\n",
      " [ 1 -1]]\n"
     ]
    }
   ],
   "source": [
    "H = linalg.hadamard(n=2)\n",
    "\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "4d92ff8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAJSCAYAAAA1TUoVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0oklEQVR4nO2df3zV1ZnnP/dyCSQYQgKXXxKZmALya1lGO1TGWiyIsouFMgO2bqdYmZn+EGcosC4u0xE6rVE7/phZa2fsskhL0cqsaF2dViyCdaysIq5YWgZjWrCTTEgCN0BCMMl3//iGayPfe7j33Pt9vuec+3m/XvelcHK+z3me55x8D/ec53linud5IIQQQgghFyQe9QAIIYQQQmyBGydCCCGEkCzhxokQQgghJEu4cSKEEEIIyRJunAghhBBCsoQbJ0IIIYSQLOHGiRBCCCEkS7hxIoQQQgjJEm6cCCGEEEKyhBsnQgghhJAs4caJEEIIIdby0ksv4YYbbsDYsWMRi8Xw1FNPhSqPGydCCCGEWMvp06cxY8YMPPTQQyLyEiJSCCGEEEJCYMGCBViwYIGYPH7jRAghhBCSJfzGiRBCCCEZOXPmDM6ePSsmz/M8xGKxfn83aNAgDBo0SGwMKrhxIoQQQkggZ86cQWlNEmg6JSbzoosuwqlT/eXdeeed2LBhg9gYVHDjRAghhJBAzp4962+ajn4VGCrwjU97F05VP4CjR49i6NCh6b825dsmgBsnQgghhFyI8kFA+eDw5Xj+f4YOHdpv42QS3DgRQgghxFpOnTqFd955J/3nhoYGvPnmm6iqqsIll1xScHncOBFCCCFEjRfzPxJycuT111/HNddck/7z6tWrAQDLly/Ho48+WqiRpeHGiRBCCCHWMmfOHHieJyaPGydCCCGEqPGQvn8UuhzDYQJMQgghhJAs4TdOhBBCCFFj8B0nafiNEyGEEEJIlvAbJ0IIIYSo4R2nNPzGiRBCCCEkS7hxIoQQQgjJEh7VEUIIIUQNL4en4TdOhBBCCCFZwm+cCCGEEKKGl8PT8BsnQgghhJAs4TdOhBBCCFHDO05p+I0TIYQQQkiW8BsnQgghhKjhHac0/MaJEEIIISRL+I0TIYQQQtTwG6c0/MaJEEIIISRLuHEihBBCCMkSHtURQgghRA3TEaThN06EEEIIIVnCb5wIIYQQoobfOKXhN06EEEIIIVnCb5wIIYQQoobpCNLwGydCCCGEkCzhN06EEEIIUcM7Tmn4jRMhhBBCSJbwGydCCCGEqOEdpzT8xokQQgghJEu4cSKEEEIIyRIe1RFCCCFEDS+Hp+E3ToQQQgghWcJvnAghhBCihpfD0/AbJ0IIIYSQLOE3ToQQQghRwztOafiNEyGEEEJIlvAbJ0IIIYRcGAvuH0nAb5wIIYQQQrKE3zgRQgghRA3vOKXhN06EEEIIIVnCb5wIIYQQooZ5nNLwGydCCCGEkCzhxokQQgghJEt4VEcIIYQQNbwcnobfOBFCCCGEZAm/cSKEEEKIGl4OT8NvnAghhBBCsiTSjdPDDz+MmpoaDB48GJdffjl+9rOfRTkcQgghhARx7o6TxMdwIts4/fCHP8SqVauwfv167N+/Hx//+MexYMECHDlyJKohEUIIIYQoiWzjdP/992PFihX40z/9U0yePBkPPvggqqur8Z3vfCeqIRFCCCEkCE/wYziRXA4/e/Ys9u3bh3Xr1vX7+/nz5+OVV1457+e7urrQ1dWV/nNvby/a2towfPhwxGLmf61HCCGEFALP83Dy5EmMHTsW8TivKUdBJBunlpYW9PT0YNSoUf3+ftSoUWhqajrv5+vq6rBx40ap4RFCCCFGc/ToUYwbN05OIPM4pYl0u/rhb4s8zwv8BumOO+5AKpVKf9L3oL73d0CiDfjvzwBtdwDz3gDGHQF+9Q3gnx5hG9vYxja2sc2tNvjvv/Ly8rBf0SQDMc/zxE8Uz549i7KyMmzfvh2f/vSn03//l3/5l3jzzTexZ88eZf/29nZUVFQAb9wDPHsF0B0HeuLAY9OArgTw8/8JnEkAj09jG9vYxja2sc2dtvoEgAqkUikMHTo05Lf177xvf/VNoHxw6PJw8gxw2Xox/XSIRyG0pKQEl19+OXbu3Nnv73fu3InZs2dn/6Da48AdL/sT6htXAw2VwIBeIO4BH2ljG9vYxja2sc2ttnhvgd/IWcLL4WkSUQlevXo1/uRP/gRXXHEFrrzySjzyyCM4cuQIvvSlL2X/kH1jgHv/E/Cvw/1J9je7gKuOAMM7gdcuBv76GraxjW1sYxvb3GmbXA8sCe/dTC5MZBunG2+8Ea2trfj617+OxsZGTJs2Dc899xzGjx+f/UM6BgK/HAH0xoBLUsAfHvU/zUP8v3++1m8b0Ou3T2gDEr3AafZjP/ZjP/ZjPwv7TWoN78WsgpfD00Ryxylf0meuzeuB9yv9SRX3/B16Wylww2eBXw8DWss+6DS8A6g5ATyzDag8A7SWsh/7sR/7sR/72dWvpA1orJa/43TwLrk7TlP+O+84hcapEuCtUf4FurEngQMjgRcu9SdbDMD17/hnxDHPn6STjwEDe4GT7Md+7Md+7Md+FvabeuzC78Yw8AQ/hmP3xumNMcCnPgs8Nt2/PPfX1wA3L/Z36Ff8G/Cjx4DPHvB38V/fBTz6lL9jf30s+7Ef+7Ef+7Gfff22P1GAlyfJB7uP6oLSEdRX+T/0kTbgM2/7Z8gDeoHPvg3Utvlt71QFhHiyH/uxH/uxH/sZ3i/ZCFTcLX9U93ad3FHdtDuMPqqze+OUWgcMKQU2zPFDNeOe/3VmT9z/b9wDNuwG/uol//zY67tgB/g/w37sx37sx37sZ1O/VBdQVceNU4Qkoh5AXgSlI/i9E8BfLABqjgN/8yIwsdWfbH/1Sf/M+H/8M/Bu5fkhnuzHfuzHfuzHfqb321sOoC7/92eueJC5fyQhI0/s3jgFpSOY1OJPtiktwPx6fyL2xIGjFcChEcCRisyhoQN7gYveBya3+P2OVACvVPtfm3YmgsNGKY/yKI/yKI/ypOTtqYr6zVv02H1UF5SOYGCPn/9iYK9/8e4crWX+xLtlEXB06PkhngN7+/r1+M+Jef5Xqa2lQLLDT3+f7Dg/bJTyKI/yKI/yKE9KXlcHIim5ckDwqG662Ud18agHkBdB6Qj2/J6fF8MD8OOP+BfvAKBhGPD2SH+SxpAhNLQHKD8L/Id/B0q7gcaLgOnNwDUN/p+DwkYpj/Ioj/Ioj/Kk5M2rRyScS4Ap8TEcuzdOGdMRlGqGhnZqho1SHuVRHuVRHuUJyGM6gsix+6guKB2Bsqp0324/KMRzUDfwsT/1d/a5hI1SHuVRHuVRHuVJyYsqHcH/u1vuqG7GOh7VhUbt8dwrTgNAfSVQd5X/5796yb+81xP3vyLM1AboVbimPMqjPMqjPMormDzzj7JcJxH1APIiKB2Bqqq0KsTzXy7xv1K9/N9yCxulPMqjPMqjPMqTkhdZOoKYzP0jC+442b1xCkpHcKFq1JlCPP+lGhh5Ojhs9HAV8H5cr/o15VEe5VEe5VFeoeQdqoz6zVv02H3HKSgdgaqqtCrEs6oTeOYx/18BHw4bfT/uT+L3B5wfbkp5lEd5lEd5lCclr7kbmLJe/o7T/nvk7jjN/G+84xQaQekILliNuidDiOe7/iT1cH7Y6C+T/qLQqn5NeZRHeZRHeZRXIHkHkyDRYvfGKWM6AgtCSimP8iiP8iiP8nKVt3RZAV6eGngQyuMUjXq5YPdRXVA6AltCSimP8iiP8iiP8nKVV59AJJnD37gHKC8NXR5OdgK/z6O68MiUjgCwIKSU8iiP8iiP8igvR3nxvmdGgSfwsYBE1APIi6B0BLaElFIe5VEe5VEe5eUqb3I9sCT81yvJjN0bp6B0BKpq1A3DZENKXxrvjzOo7dAIf1FkemaiN/dq29SP+lE/6kf93NZvUmt471SSFXbfcQpKR6CqRt3Zt08s7ZYJKe3qkzeo+/y29wf0yYuf/8yaE8Az2/zikLlU26Z+1I/6UT/q57Z+JW1AY7X8Had99wIXCdxxOtUJXH477ziFRlA6AlU16sPDgWSHXEjpNQ3+oghq+8SvgeODFc/szb3aNvWjftSP+lE/t/WbegwkWuzeOGVMR6AKDdVt0w0p1W3TDX2lftSP+lE/6uesftufKMDLUwNP8GM4dh/VBaUjUIWG6rbphpTqtumGvlI/6kf9qB/1c1u/ZCNQcbf8Ud3rgkd1V/CoLjwypSOIe3pho6o2QC+kVLcNKLwO1I/6UT/qR/0s1y9WmPdnrvAbpzSJqAeQF0HpCFShobptuiGlum26oa/Uj/pRP+pH/dzWb285gDr59y1JY/fGKSgdwYVCQzOFqYZVNbs7Llttm/pRP+pH/aifu/odqozmfevF/I+EHMOx+45TUDoCVWioKkw1jCrWDcOAG27yoyMkqm1TP+pH/agf9XNbv+ZuYMp6+TtO//dbcnec/uC/8o5TaASlI7hgaGivcNXsuHC1bepH/agf9aN+zup3MIlI8AQ/hmP3xiljOgJVaKgqTFU3pFT3mbohs9SP+lE/6kf9ilK/pcsK8PIk+WD3UV1QOgJVaKiy4rSin3YVa82x6Ia+Uj/qR/2oH/VzW7/6BIAK+aO6vYJHdbPMPqqze+OUWgcMKQU2zPFDQ+Oe//VpT9z/b9wDNuz2wzh7+y62nQsd7Yln3+9oBTB7hf81aqGeqepXnQJe2eR/TVyoZ1I/6kf9qB/1s18/pIDeSvmN06t/K7dx+thaozdOiagHkBdB6QhUoaGqMNUwqljrjkU39JX6UT/qR/2on9v6Ta4HlhT+dUqyx+6NU1A6AlX160Mj/HBPqSrWl7X4//rJNWT26t/48nKt7k39qB/1o37Uz239ZjYU5v2ZK0xHkMbuo7qgdASq6tdHKoBbFgFHh+YWUqpbxbo6BWx+Gqhuzy1kdlC3/3NdifPbqB/1o37Uj/oVr35dp6MpufLz++SO6q5cY/RRXTzqAeRFUDoCVfXrt0f6iyIGmSrW05v9RRE0FlXI7Is1/mLMtbo39aN+1I/6UT+39dt5KSLBE/wYjt0bp4zpCHRDUVUhpbptuiGzum3Uj/pRP+pH/ZzVj+kIIsfuo7qgdAQmVbF2vUo39aN+1I/6UT9Z/aJKR/Av98sd1f3hah7VhUbtccOrWPeZ19kq3dSP+lE/6kf9RPWL9+b/7iR5kYh6AHkRlI7ApCrWrlfppn7Uj/pRP+onq19U6Qg8yNw/kpCRJ3ZvnILSEZhUxdr1Kt3Uj/pRP+pH/WT1m9Qa9Zu36LH7jlNQOgKTqli7XqWb+lE/6kf9qJ+sfiVtQGO1/B2nl++Xu+N0Fe84hUdQOgKTqli7XqWb+lE/6kf9qJ+sflOPIRLOJcCU+BiO3RunjOkIdENKVaGhqmfqhszqPpP6UT/qR/2oX1Hqt/2JArw8ST7YfVQXlI7ApCrWrlfppn7Uj/pRP+onq1+yMZrM4S89IHdUd/VXjT6qs3vjlFoHDCk1t4q161W6qR/1o37Uj/rJ6pfqAqrquHGKkETUA8iLoHQEJlWxdr1KN/WjftSP+lE/Wf32lgOoK/z79EJI3T+y4I6T3RunoHQEqurXDcP0QkNVz3ShSrfqmap+NSf0Qnvrq4DORO52of/oP/qP/it2/+2pCu+dSrLC7qO6oHQEqurXnX37xNLu3EJDVc90oUq36pmqfoP79DuTyO2ZyQ7/bD/ZkZtd6D/6j/6j/4rdf10diKTkyu4H5Y7q5qwy+qguHvUA8iIoHYGq+vXh4f4vi6A2VWio6pkuVOlWPVPVb2Ir0FKmN87S7tztQv/Rf/Qf/Vfs/ptXDxItdm+cMqYjUIV42tKmGzJrS5tuiK4tbfSf3W30n91tDvsvqnQE5+44SXwMx+6juqB0BGFUzZZuk67SLd0mXfWc/qP/6D/6zzQ/6LZFlY7gxb+TO6q75i95VBcatcflqmZLtwGyVbql2wBzbE3/0X/0n5l+oP8C2iL6RsYT/BhOIuoB5EVQOoIwqmZLt0lX6ZZuk656Tv/Rf/Qf/WeaH3TbokpHQNLYvXEKSkcQVtVs6X7SVbql+0lXPaf/6D/6j/4zzQ86/Q5VRv3mLXrsvuMUlI4gjKrZ0v2kq3RL95Ouek7/0X/0H/1nmh90+zV3A1PWy99x+unfy91xmvsXvOMUGkHpCEKrmi3dL1Mo6rshVemW7pchBNk4P9B/9J/JfqD/is5/B5Mg0WL3xiljOgJVKKoqDNekfqowVVt00A15tkUH+s+MsdB/9F8R+W/psgK8PDXxBD4WYPdRXVA6gjCqZkv3k67SLd1Puuo5/Uf/0X/0n006KPslEEnm8J/+PTBE4KjutPlHdXZvnFLrfEdumBNu1WzpftJVuqX7SVc9p//oP/qP/rNJB1U/pIDeSvmN0wv/Q27jNO82ozdOiagHkBdB6QjCqJot3U+6Srd0P+mq5/Qf/Uf/0X826aDqN7keWFK41yjJHbs3TkHpCMKomq1bTVy6SvdlLf6/JnMNQa5t8wtOSsmTrnpO/9F/9B/954r/ZjbIvF8/jAeZO0gSMvLE7qO6oHQEYVTN1q0mLl2luzoFbH4aqG7PLQT5WJl/J6ClTEaedNVz+o/+o//oP1f813U6mpIrOwWP6q41+6guHvUA8iIoHUFYVbN1qolLV+me3uwvwqCxXKgK+ZmErDzpquf0H/1H/9F/Lvhv56WIBC8m9zEcuzdOGdMR6Iba6oawhiFPt016nLrydNvoPzPsQv+ZMU76r+j8F2U6AgLA9qO6oHQEYVTNDiO016Sq59J2ka56Tv/Rf/Qf/WfaOLXlJRBJOoLnH5I7qpu/kkd1oVF7XK5qNmBHlW7pcerK022j/8ywC/1nxjjpv+LzX7z3wu/GMPAEP4aTiHoAeRGUjiCMqtlhhPaaVPVc2i7SVc/pP/qP/qP/TBunrjymI4gcuzdOQekIwqqarVtNXLpKt/Q4paus039m2IX+M2Oc9F/x+W9SazTvW6mL2xZcDrf7jlNQOoIwqmbrVhOXrtItPU7pKuv0nxl2of/MGCf9V5z+K2kDGqvl7zj9+Ntyd5yuv5V3nEIjKB1BaFWzM4SGGlelW3qcuvJ6LRkn/Uf/mTxO+q/o/Df1GCLBE/wYjt0bp4zpCGwJmVXJUz3TpHHqytPVnf4zwy70nxnjpP+Kzn/bnyjAy5Pkg91HdUHpCGwJmTWp6rm0XaSrntN/9B/9R/+ZNk5decnGaDKHP/ew3FHdf/oKj+pCI1M6AsD+UFTVM00ap648Xd3pPzPsQv+ZMU76rwj9FwOJlkTUA8iLoHQEtoTMmlT1XNou0lXP6T/6j/6j/0wbp668veUA6sJ7r2bCg8z9IwkZeWL3xikoHYGq2nbDMNlQVN0q3apnqvrVnJANmdWtsq56Jv1H/9F/9B/9l9kPhyrDe6eSrCj4HacNGzZg48aN/f5u1KhRaGpqAgB4noeNGzfikUcewfHjxzFr1ix8+9vfxtSpU7OWoUxHoKq23dm3TyztlglF1a3SrXqmqt/gbv//zyRkQmZ1q6yrnkn/0X/0H/1H/2X2Q3M3MGW9/B2nZwXvOP3nIrzjNHXqVDQ2NqY/Bw4cSLfde++9uP/++/HQQw/htddew+jRo3Httdfi5MmTuQsKSkegqrZ9eDiQ7JALRdWt0q16pqrfxFagpUwuZFa3yrrqmfQf/Uf/0X/0X2Y/HEwiEs4lwJT4GE4oG6dEIoHRo0enP8mk72jP8/Dggw9i/fr1WLJkCaZNm4YtW7ago6MD27Zty11QxnQEqhBP3TbdUFRb2nRDZm1po//sbqP/7G6j/wrWtnRZ7u9KUlBCOar71re+hYqKCgwaNAizZs3CXXfdhUsvvRTvvvsuamtr8cYbb2DmzJnpPosWLcKwYcOwZcuWrGQo0xGYVBXcljbpKuvSbfSfGX6g/+i/qG3tgv/qEwAq5I/qnvmO3FHdDV8urqO6WbNm4Xvf+x5+8pOf4Lvf/S6ampowe/ZstLa2pu85jRo1ql+f370DFURXVxfa29v7fQBkTkcQ98ypCm5LGyBbZV26jf6zu43+s7uN/itcW7xX+/1MCkOi0A9csGBB+v+nT5+OK6+8ErW1tdiyZQs+9rGPAQBisf5nmJ7nnfd3v0tdXd15F84BBKcjMKkquC1t0lXWpdvoPzP8QP/Rf1Hb2gX/Ta4HluT3ntZC6v6RBXecCr5x+jBDhgzB9OnTcfjwYSxevBgA0NTUhDFjxqR/prm5+bxvoX6XO+64A6tXr07/ub29HdXV1cHpCEyrCq5bpVu6n3SVdel+9J8ZfqD/6D+T/WCD/ya1hva+JtkResmVrq4u1NbW4s///M/xta99DWPHjsVXv/pV3H777QCAs2fPYuTIkbjnnnvwxS9+MatnKtMRmFQVXLdKt3Q/6Srr0v3oPzP8QP/Rfyb7wRb/lbQBjdXyd5ye/ge5O06LvlRcd5zWrl2LPXv2oKGhAXv37sUf//Efo729HcuXL0csFsOqVatw1113YceOHXj77bdx8803o6ysDDfddFPuwoLSERhXFTxuST/dUGJde0r3o//M8AP9R/+Z7AcL/Df1GEi0FHzj9N577+Gzn/0sJk2ahCVLlqCkpASvvvoqxo8fDwC4/fbbsWrVKnzlK1/BFVdcgd/+9rd4/vnnUV5enruwjOkIVKGvqrBR3ZBZ3Wea1E8VFmuLDvSfGWOh/+g/+i+0sWx/Ivd3ZSE4d8dJ4mM4oR/VhYEyHYFJVcHDGItJVcFN0oH+o/9M1oH+o/8KNZZkI1Bxt/xR3VP/KHdUt/iLRh/V2b1xSq3zHblhjh+qGff8rzN74v5/4x6wYbcfNtrbt5M9FxrbE8++39EKYPYK/2vbQj3TpH7VKeCVTf7XxFGPhf6j/+g/+3Sg/+T8l+oCqurkN047BDdOnzZ745SIegB5EZSOwKSq4GGMxaSq4CbpQP/RfybrQP/Rf4Uay95yAHWFe4+SnLF74xSUjkBZVXqEH14qVaX7shb/XyO5hrDWtvkFJ6XkaVfp1rRnfZVfMFRKHv1H/9F/9J8r/ttTFf67lSix+6guKB2Bqqr0kQrglkXA0aG5hczqVumuTgGbnwaq23MLYT1W5t8JaCmTkadbpVvXnskO/05AskNGHv1H/9F/9J8r/uvqQCQlV558RO6obsmfG31UF496AHkRlI5AVVX67ZH+JI1Bpkr39GZ/EQaN5UJVs88kZOXpVOnOx56l3bLy6D/6j/6j/1zw37x6kGixe+OUMR2BbmioKmRWty2MENYw5Om2SY+T/qP/TB4n/Uf/hTzOyNIRCH4Mx+6juqB0BCZV6Q4jhNWoKt3C46T/6D+Tx0n/0X8S44wqHcH/Fjyq+yMe1YVH7XHzq3QDstXEdeXptkmPk/6j/0weJ/1H/4U+zlh278dCw2+c0iSiHkBeBKUjMKlKdxghrCZV6ZYeJ/1H/5k8TvqP/pMYJ9MRRI7dG6egdASmVenWrZqtO07pKuvS46T/6D+Tx0n/0X9hj/NQZXTvXC8WnWyDsPuOU1A6ApOqdOtWzdYdp3SVbulx0n/0n8njpP/oP4lxNncDU9bL33H6p0eAsrILd8iXjg7gj3nHKTyC0hEYV6U7QwhraOPUlddryTjpP/rP5HHSf/RfyOM8mEQkeIIfw7F745QxHUEYoaiqZ9oSMqurg0njpP/oP5PHSf/RfyGPc+myArw8ST7YfVQXlI7ApCrdJoXMmlT1XNou9J8Z46T/6D+Tx2mL/+oTiCRz+HbBo7qlPKoLj0zpCIDCh6KqnqnqF/fMCZnV1cGkcdJ/9J/J46T/6L+wxxnve6Y0XkzuYziJqAeQF0HpCEyq0m1SyKxJVc+l7UL/mTFO+o/+M3mctvhvcj2wJLzXKrkwdm+cgtIRqKpRNwzTC0VVPTOMqtm649SVp6u77jh1Q3R1q6yrnqnqV3NCdpy68nR1z8eemcKvuf64/rj+wh3npNZCv0mzQ+rbIAu+cbL7jlNQOgJVNerOvn1iaXduoaiqZ4ZRNVt3nLrydHXXHaduiK5ulXXVM1X9BvfZ80xCZpy68nR11x2nKvya64/rj+sv3HGWtAGN1fJ3nH74P+XuON34p7zjFBpB6QhU1agPDweSHbmHoqqeGUbVbN1x6srT1V13nLohurpV1lXPVPWb2Aq0lMmNU1eeru552bNXbr5w/XH9cf198MypxxAJnuDHcOzeOGVMR6AKKWWbGW26IbpsM6NNN6SbbWa0cf1Z27b9iQK8PEk+2H1UF5SOIIwq1mwzoyo428xok64izzauP7Z90JZsBCrulj+qe2yT3FHdZ1fwqC40MqUjiHuyFa7ZZn6VdbaZUbWebWa0cf1Z3BaL5n1L0iSiHkBeBKUjCKOKNdvMqArONjPapKvIs43rj20ftO0tB1An/771IHP/SEJGnti9cQpKRxBWFWv2i74qeD5V1tnPjKr1tsxP1/tx/dnb71Bl1G/eosfuO05B6QjCqGLNfmZUBdetss5+ZlStt2V+ut6P68/ufs3dwJT18nectgnecbqJd5zCIygdQWhVrNkv8qrgeVVZZz8zqtbbMj9d78f1Z22/g0lEwrkEmBIfw7F745QxHYEq1FYVRs1+RlcFN06HYu5XatBY2I/rr4j6LV1WgJcnyQe7j+qC0hGEUcWa/QypCm6QDsXcT7qKPPtx/bHf7/RLAKiQP6rb+r/kjuo+d4vRR3V2b5xS64AhpcCGOX6obdzzv47uifv/jXvAht1+GGdv31eA50Jxe+LsF1W/oxXA7BX+19JRj4X9cu9XnQJe2eQfC0U9Fvbj+iu2fkgBvZXcOEVIIuoB5EVQOoIwqliznxlVwU3SoZj7SVeRZz+uP/b7oN/kemBJYV6hOSF1/8iCO052b5yC0hGoqokfGuGHe+Yaoltf5RfUzLV6OeUVtir4ZS3+tx25hsjXtvkFPCmvMPKu/o3vP1fnp+vyuP7sljezIZr3LUlj91FdUDoCVTXxIxXALYuAo0NzC9FNdvh3OpId54f9Ul7u8nSrglengM1PA9XtuYXIHyvz73S0lFFeIeQN6vZ/rivh5vx0XR7Xn93yuk5HU3Lle5vljuo+/wWjj+riUQ8gL4LSEaiqib890v8lE0PuVaxLu3OvXk55ha0KPr3Z/yUTNJYLVa0/k6C8QsrrSrg9P12Xx/Vnr7ydl4JEi90bp4zpCHRDpXVDdClPrmK463axRZ5um+t2sUUe15+18qJKR3DujpPEx3DsPqoLSkcQRlXwMEKzi1mebpvrdrFFnm6b63axRR7Xn+XyEogkHcGWR+WO6pbfzKO60Kg9LlcVHJCtQu66PN021+1iizzdNtftYos8rj975cV7z38XSuAJfgwnEfUA8iIoHUEYVcHDCM0uZnm6ba7bxRZ5um2u28UWeVx/dsuLKh0BSWP3xikoHUFYVcF1q8FTXmGribtuF1vkJXrtGCflcf25Jm9Sa9Rv3qLH7jtOQekIwqgKrlsNnvIKWxXcdbvYIq/mBPDMNj/ay+RxUh7Xn4vyStqAxmr5O06bt8jdcfrCct5xCo2gdAShVQUfIFyF3HV5cUvGSXnB8notGSflcf05Jm/qMZBosXvjlDEdgcOhqM7IUz3TpHFSXmF967pdbJHH9WetvO1PFODlqYEn+DEcu4/qgtIROB+K6oA83argrtvFFnm6vnXdLrbI4/qzW16yMZrM4f9L8KjuFh7VhUemdASAu6GorshTPdOkcVJeYX3rul1skcf1Z7G8GCKhEIkts/0YTiLqAeRFUDoC10NRXZCnWxXcdbvYIk/Xt67bxRZ5XH92y9tbDqAu3HcrUWL3xikoHYGqmnjDMNlQVN0q5Kpwb1f0y/RMVb+aE7IhwVf/xtcvVz+owr1d0S/TM12Zn1x/ds9Pl9ffocrw3qkqPMjcP5KQkSd233EKSkegqibe2bdPLO2WCUXVrUKuCvd2QT/VM1X9Bnf7/38mIRMSPKhPXlciNz+owr1d0E/1TBfmJ9ef3fPT9fXX3A1MWS9/x+l/fk/ujtOffp53nEIjKB2Bqpr44eFAskMuFFW3Crkq3NsZ/TI8U9VvYivQUiYXEvxijf9LLVc/qMK9ndEvwzOdmZ9cf3bPT4fX38EkIoF3nNLYvXHKmI5AFeKp26YbiqrbphvCaot+0m2u21NaP+k21+3J9We3PQX1W7qsAC9Pkg92H9UFpSMwqaq7dBVyW/STbnPdntL6Sbe5bk+uP7vtKa5fAkCF/FHdd78vd1T3Z3/Co7rQyJSOIO6ZU9Vdtw2wo2q9LW2u21NaP+k21+3J9We3PSX1i/cW/l2aDZ7gx3ASUQ8gL4LSEZhU1V26Crkt+km3uW5Paf2k21y3J9ef3faU1m9yPbAkmlcu8bF74xSUjsC0qu7SVcht0U+6n+v2lNZPup/r9uT6s9uekvpNao3ohRsTurgtISM/7L7jFJSOwKSq7tJVyG3RT7qf6/aU1k+6n+v25Pqz257S+pW0AY3V8necHtkKlArccersAP78c7zjFBpB6QiMq+oed7xqva5+0v1ct6e0ftL9XLcn15/d9hTUb+oxRIIn+DEcuzdOGdMRqEJRVWGjuiGsus/UDWF1QT/pfq7bU1o/6X6u25Prz257Cuq3/YkCvDxJPth9VBeUjsCkqu7SVcht0U+6n+v2lNZPup/r9uT6s9ue0volG4GKu+WP6v7hB3JHdV/6L0Yf1dm9cUqtA4aUAhvm+KGacc//OrMn7v837gEbdvtho719F9vOhar2xLPvd7QCmL3C/9q2UM9U9atOAa9s8r+2LdQzTdJPup/r9pTWT7qf6/bk+rPbntL6pbqAqjpunCIkEfUA8iIoHYFJVd2lq5Dbop90P9ftKa2fdD/X7cn1Z7c9pfXbWw6grrDv0mzwIHP/SEJGnti9cQpKR6CsKj3CDy+VqkJ+WYv/rxGxqtmW6Ffb5hfUpD3t1K++yi9oS3vaqR/Xn9367akq/LuU5ITdR3VB6QhUVaWPVAC3LAKODs0thFW3Cnl1Ctj8NFDdnlsIq27VbFv0O1bm30FoKaM9bdQv2eHf6Uh20J426sf1Z7l+HYik5Mp3BI/qvmz2UV086gHkRVA6AlVV6bdH+pM0Bpkq5NOb/UUYNJYwqmbbpN+ZBO1ps36l3bSnzfpx/dmr37x6RIIXk/sYjt0bp4zpCHRDQ1UhrLptuiGsum226BdGKHEx21NaP2m/u25Prj+77SmoH9MRRI7dR3VB6QiKumq2JfqFEUpczPaU1k/a767bk+vPbntK6xdVOoJvb5M7qrv1Jh7VhUbtcVbNtlU/3XHSnmboJ+131+3J9We3PUX1M/8oy3USUQ8gL4LSERRz1Wxb9AsjlLiY7Smtn7TfXbcn15/d9pTWL7J0BDGZ+0cW3HGye+MUlI6gmKtm26Sf7jhpTzP0k/a76/bk+rPbnpL6HaqM+s1b9Nh9xykoHUExV822RT/dcdKeZugn7XfX7cn1Z7c9pfVr7gamrJe/4/Q/HpO743TbZ3nHKTSC0hEUc9Vsq/TTHSftaYZ+0n533Z5cf3bbU1C/g0mQaLF745QxHUEYobaqZ4YREqw7Flv00x0n7WmGftJ+d92eXH9221NQv6XLCvDy1ODcHSeJj+HYfVQXlI6gmKtm26KfdIi16/aU1k/a767bk+vPbnuK65dAJJnD//5xuaO6v/iM0Ud1dm+cUuuAIaWsmm2bfrrjpD3N0E/a767bk+vPbntK64cU0Fspv3H6O8GN01+avXFKRD2AvAhKR1DMVbNt0U86xNp1e0rrJ+131+3J9We3PaX1m1wPLAnnlUqyw+6NU1A6AlU16oZheqGoqmeaVDVbtyq46pmqfjUn9EJtdcepK09Xd1uqrOvOXd31oDtOrj+uP66//NfDzIbCvkdJzth9VBeUjkBVjbqzb59Y2p1bKKrqmUZVze7TL9eq4KpnqvoN7tPvTEJmnLrydHW3pcq67tzVXQ+64+T64/rj+vugn+566DodTcmVB38od1S36kajj+riUQ8gL4LSEaiqUR8eDiQ7cg9FVT3TpKrZulXBVc9U9ZvYCrSUyY1TV56u7rZUWdedu7rrQXecXH9cf1x/+a+HnZeCRIvdG6eM6QhUIaW2tOmGsLLNjDbdkGe2mdHG9Wd3m8PrL7J0BIIfw7H7qC4oHUEYVayl26SrgrPN7irrbOP6Y1vxrL+o0hE8IHhU91XHjupeeukl3HDDDRg7dixisRieeuqpfu2e52HDhg0YO3YsSktLMWfOHPziF7/o9zNdXV247bbbMGLECAwZMgSf+tSn8N577+U++trjclWspdsA2argbLO7yjrbuP7YVhzrL96b+7uyEPAbpzSJXDucPn0aM2bMwBe+8AX80R/90Xnt9957L+6//348+uijmDhxIr7xjW/g2muvxaFDh1BeXg4AWLVqFZ555hk8/vjjGD58ONasWYOFCxdi3759GDBgQPaDCUpHEEYVa+k26argbLO7yjrbuP7YVjzrj+kIIifnjdOCBQuwYMGCwDbP8/Dggw9i/fr1WLLE9+yWLVswatQobNu2DV/84heRSqWwadMmfP/738e8efMAAFu3bkV1dTVeeOEFXHfdddkPJigdQVhVrKX7SVcFZz+7q6yzH9cf+xXH+pvUmutruzB4MYiUQ5GQkSd53XGKxWLYsWMHFi9eDAB49913UVtbizfeeAMzZ85M/9yiRYswbNgwbNmyBbt27cLcuXPR1taGysrK9M/MmDEDixcvxsaNG8+T09XVha6urvSf29vbUV1dHZyOIIwq1tL9pKuCs5/dVdbZj+uP/T7o5/r6K2kDGqvl7zjd94TcHac1y9y646SiqakJADBq1Kh+fz9q1Kh0W1NTE0pKSvptmj78Mx+mrq4OFRUV6U91dbXfEJSOILQq1tL9MoTohlYVnP3srrLOflx/7FcU62/qMUSCJ/gxnIJunM4Ri/X/qs3zvPP+7sOofuaOO+5AKpVKf44ePeo3ZExHoAolVoWpmtRPN9TWJB2KuV+pQWNhP66/Yuvn8Prb/oTyXUrCx4qjug+jTEcQRhVr6X7SVcHZz+4q6+zH9cd+xbP+ko3RZA7/W8GjurVmH9UVdOPkeR7Gjh2Lr371q7j99tsBAGfPnsXIkSNxzz33pC+HJ5NJbN26FcuWLQMANDY2Yty4cXjuueeyuhyedmRqHTCkNPwq1tL9pKuCs19h+0lXWWc/rj/2K571l+oCqurkN07f2i63cfqvS43eOCVy7XDq1Cm888476T83NDTgzTffRFVVFS655BKsWrUKd911FyZMmIAJEybgrrvuQllZGW666SYAQEVFBVasWIE1a9Zg+PDhqKqqwtq1azF9+vR0lF3WBKUjCKOKtXQ/6arg7Gd3lXX24/pjv+JZf3vLAdTl9q4kBSXnjdPrr7+Oa665Jv3n1atXAwCWL1+ORx99FLfffjs6Ozvxla98BcePH8esWbPw/PPPp3M4AcADDzyARCKBZcuWobOzE3PnzsWjjz6aWw4nIDgdQRhVz+ur/IKMpldZ160KXtvmF/CkvMLI062ybsv8dF0e15/d8lxff3uqcn1tFwamI0hjd8mVoHQEYVQ9T3b4Z+bJjvPDRk2qsq5bFfxYmX+no6WM8gohT7fKui3z03V5XH92y3N9/XV1IJKSK/f+k9xR3e1/bPRRXTzqAeRFUDqCsKqel3abX2U9n6rgZxKUF3WVdZvmp+vyuP7slufy+ptXj0jwBD+GY/fGKWM6At1QVN0Q5DDk6bZJj5PyCtvmul1skcf1Z7c8h9cf0xFEjt1HdUHpCMKoeh5G6KtJVcFtsYst8rSrnjtuF1vkcf3ZLc/19RdVOoJ7/gkYPCR0eThzGvhvPKoLj9rjclXPATuqrEuPk/IK2+a6XWyRx/Vntzyn118s+H1IxEhEPYC8CEpHEEbV8zBCX02qCm6LXWyRp9vmul1skcf1Z7c819dfVOkIPMjcP5KQkSd2b5yC0hGEVfVct9q2dJV16XFSXmGrpbtuF1vkcf3ZLc/l9XeoMuo3b9Fj9x2noHQEYVQ91622LV1lXXqclFfYqueu28UWeVx/dstzff01dwNT1svfcbpb8I7TOt5xCo+gdAShVT0fYEmVdelxUl5hq6W7bhdb5HH92S3P4fV3MIlIOJcAU+JjOHZvnDKmI7AgpNS4quC22MUWebq+dd0utsjj+rNbnsPrb+myArw8ST7YfVQXlI7AlpBSk6qC22IXW+Tp+tZ1u9gij+vPbnmur7/6BCLJHH7X/5Y7qvvvf8SjutDIlI4AsCCk9ALyVM80aZyUV1jfum4XW+Rx/dktz+X1F+97JomMRNQDyIugdAS2hJSaVBXcFrvYIk/Xt67bxRZ5XH92y3N9/U2uB5aE+2oNROr+kQV3nOzeOAWlI1BVo24YJhtSqltlXfVMVb+aE7KhtrpVyFXh3q7ol+mZrszPTOHerujH9We3fi6vv0mt4b1TSVbYfccpKB2Bqhp1Z98+sbRbJqRUt8q66pmqfoO7/f8/k5AJtdWtQq4K93ZBP9UzXZifqnBvF/Tj+rNbP9fXX0kb0Fgtf8fpG0/K3XH6qyW84xQaQekIVNWoDw8Hkh1yIaW6VdZVz1T1m9gKtJTJhdrqViFXhXs7o1+GZzozP3sd14/rz279HF5/U4+BRIvdG6eM6QhUoaG6bbohpba06Yba2mJPaf2k21y3J9ef3fbk+itY2/YnCvDy1ODcHSeJj+HYfVQXlI7ApKrZtrTZUoXcFv2k21y3J9ef3fbk+iusPZONQMXd8kd1f7ND7qjua5/mUV1oZEpHEPfMqZptSxtgRxVyW/STbnPdnlx/dtuT66+AbRF9I+MJfgwnEfUA8iIoHYFJVbNtabOlCrkt+km3uW5Prj+77cn1V1h77i0HUBfNO5cAsH3jFJSOwLSq2bpV1qX72VKF3Bb9pPu5bk+uP7vtyfVXOHseqoz6zVv02H3HKSgdgUlVs3WrrEv3s6UKuS36Sfdz3Z5cf3bbk+uvsPZs7gamrJe/47TxKbk7Tncu5h2n0AhKR2Bc1ey4Jf10Q21tsae0ftL9XLcn15/d9uT6K5g9DyZBosXujVPGdASq0FBVWKxuSKnuM03qpxtqa4s9pfWT7ue6Pbn+7LYn11/BxrJ0WQFenhp4gh/DsfuoLigdgUlVs8MYi0lV3W2xp7R+0v1ctyfXn9325Por8FgSACrkj+o2PCV3VLfB7KM6uzdOqXXAkFJgwxw/NDTu+V+f9sT9/8Y9YMNuP4yzN+afH58LHe2JZ9/vaAUwe4X/NWqhnmlSv+oU8Mom/2tpibFI21NaP+l+rtuT689ue3L9FXYsSAG9lfIbpzuflts4bVxk9MYpEfUA8iIoHYFJVbPDGItJVd1tsae0ftL9XLcn15/d9uT6K+xYJtcDSwr7KiW5YffGKSgdgar69aERfrinVBXry1r8f23lGqJb2+YX1JSSp1uF3BZ7SutXX+UXDKU97dSP689u/VxffzMbwnmfXggPMvePJGTkid1HdUHpCFTVr49UALcsAo4OzS2kVLeKdXUK2Pw0UN2eW4jusTL/zLylTEaebhVyW+wprV+yw7/TkeygPW3Uj+vPbv1cX39dp6MpufLXT8sd1X3d7KO6eNQDyIugdASq6tdvj/QXRQwyVaynN/uLImgsF6oKfiYhK0+nCrlN9pTWr7Sb9rRZP64/u/Vzef3tvBSR4MXkPoZj98YpYzoC3dBXVUipblsYIbphyNNts8We0vpJ+911e3L92W1Prr+CjTOqdAQkjd1HdUHpCEyqYh1GiK5JVchtsae0ftJ+d92eXH9225Prr8DjTCCSdARf+5HcUd3ffIpHdaFRe9zwKtZ95pWstq0rT7fNFntK6yftd9ftyfVntz25/go3znhv7u/KQuAJfgwnEfUA8iIoHYFJVazDCNE1qQq5LfaU1k/a767bk+vPbnty/RV2nExHEDl2b5yC0hGYVMU6n6rguuNklXUz9JP2u+v25Pqz255cf4Ub56TWaN63Uhe3Lbgcbvcdp6B0BCZVsdatCq47TlZZN0M/ab+7bk+uP7vtyfVX2HGWtAGN1fJ3nNY/I3fH6Zs38I5TaASlIzCpinVeVcF1x6krr9f8quBW6Sftd9ftyfVntz25/go2zqnHEAme4Mdw7N44ZUxHEEYIq+qZtoQ86+rggj2l9ZP2u+v25Pqz255cfwUb5/YnCvDyJPlg91FdUDoCk6pY2xJqa0tVcFv0k/a76/bk+rPbnlx/hR1nsjGazOF3/B+5o7q6hTyqC41M6QiAwoeGqp6p6hf37Ai11dXdFntK6yftd9ftyfVntz25/go4zhhItCSiHkBeBKUjMKmKtS2htrZUBbdFP2m/u25Prj+77cn1V9hx7i0HUBfOO1WFB5n7RxIy8sTujVNQOgJV9euGYXqhoapnhlEVXHecuvJ0dbfFntL6Sftd2p7SVeRVz1T1qzmhF3quO05debq622JPaf2k/S5tzz1VhX6Tkhyx+45TUDoCVfXrzr59Yml3bqGhqmeGURVcd5y68nR1t8We0vpJ+13antJV5FXPVPUb3KffmYTMOHXl6epuiz2l9ZP2u/j87EAkJVfWCd5xupt3nMIjKB2Bqvr14eFAsiP30FDVM8OoCq47Tl15urrbYk9p/aT9Lm1P6Sryqmeq+k1sBVrK5MapK09Xd1vsKa2ftN+l7TmvHtEQ+yAJZpgfxCLSL3vs3jhlTEegCvFkG9vYln+bbkg329jGtrzamI4gcuw+qgtKRxBGlW62sY1tH7RJV5FnG9vYFn06gv/2LDBI4Kiu6zRwz3/mUV1oZEpHEPdkK3izjW3F1gbIVpFnG9vYxnQEhpCIegB5EZSOIIwq3WxjG9s+aJOuIs82trHNgHQEMbDIr4/dG6egdARhVelmP/Zjv+iqyLMf+7FfXzqQyqjfvEWP3XecgtIRhFGlm/3Yj/0+6CddRZ792I/9PujX3A1MWS9/x+m/Pid3x+lb/4l3nEIjKB1BaFW62Y/92C+6KvLsx37sh4ZhwMEkSLTYvXHKmI5AFSrdIVzdm/3Yz8V+nQaNhf3Yr4j6LV1WgJenBufuOEl8DMfuo7qgdARhVOlmP/Zjv+xCpW3Rgf3Yz9p+CUSSOXztP8sd1f3tAqOP6uzeOKXWAUNKgQ1z/FDpuOcfJ/TE/f/GPWDDbj+Ms7dvJ3sulLonzn7sx346/Y5WALNX+McKUY+F/div2PohBfRWym+c1ghunO4ze+OUiHoAeRGUjiCMKt3sx37sl12otC06sB/72dpvcj2wpBAvUKKL3RunoHQEYVR1r6/yC7BKVZGnPMozWZ5uFfnLWoDqVObQ7PfjwTrUtvkFUSmP8igPmNkQ7XuXWH5UF5SOIIyq7skOv/xEsuP8sG3Ko7xik6dbRb46BWx+GqhuP38s78f7+g04X4djZf6dqpYyyqM8yus6HU3JldU/ljuqu/96o4/q4lEPIC+C0hGEVdW9tFu2ijzlUZ7J8nSqyE9v9l8SQWP5ZdJ/SQTp8GKN/2fKozzKA3ZeChItdm+cMqYjKBUOsaY8yis2ebptrtuF8igvZHmRpSMQ/BiO3Ud1QekIwqjqrgoNpTzKKzZ5um2u24XyKE9EXgKRpCP4quBR3QM8qguP2uNyVd0B2SrylEd5JsvTbXPdLpRHeWHLi/ciEviNU5pE1APIi6B0BGFUdVeFhlIe5RWbPN021+1CeZQnIY/pCCLH7o1TUDqCsKq6ZwoNpTzKK0Z5utXgXbcL5VFe2PImtUbzvvViECmHIiEjT+y+4xSUjiCMqu6q0FDKo7xik6db1d11u1Ae5UnIK2kDGqvl7zj95U/k7jj93XW84xQaQekIQqvqPkC4ijzlUZ7J8uKWjJPyKM8xeVOPIRLOfeMk8TEcuzdOGdMROByKSnmUZ4Q81TNNGiflUZ5j8rY/UYCXJ8kHu4/qgtIROB+KSnmUF7E83arurtuF8ihPQl6yMZrM4X/xvNxR3d/P51FdaGRKRwC4G4pKeZRngjzVM00aJ+VRnnPyIjrK8gQ/hpOIegB5EZSOwPVQVMqjvKjl6VZ1d90ulEd5EvL2lgOoC/vtShTYvXEKSkegqgbfMEw2FFW3ivyENv+r2Vyr3VM/6ielX6ZnqvrVnJAN6b76N75+ufpBlW6B+lG/qPU7VBneO1UF0xGksfuOU1A6AlU1+M6+fWJpt0woqm4V+ZoTwDPb/OKsuVS7p37UT0I/1TNV/QZ3+/9/JiET0j2oT15XIjc/qNItUD/qF7V+zd3AlPXyd5xW7pS74/TQtbzjFBpB6QhU1eAPDweSHXKhqLpV5Ccf8xdMrtXuqR/1E9MvwzNV/Sa2Ai1lciHdL9b4L6Vc/aBKt0D9qF/U+h1MIhI8wY/h2L1xypiOQBXiqdumG4qq26Ybwkr9qJ+EftJtrtuT+lG/LNuWLivAy5Pkg91HdUHpCFQhnrpttlSRp37UT0I/6TbX7Un9qF9O+iUAVMgf1X3lBbmjuofn8aguNDKlI4h7euGfqjbAjiry1I/6Segn3ea6Pakf9cu2Ld5buHco0SIR9QDyIigdgSrEU7fNliry1I/6Segn3ea6Pakf9culbXI9sCSC960HmftHEjLyxO6NU1A6gguFZmcKE3elijz1o34S+kn3c92e1I/6ZdtvUmvUb96iJ+c7Ti+99BK+9a1vYd++fWhsbMSOHTuwePHidPvNN9+MLVu29Osza9YsvPrqq+k/d3V1Ye3atXjsscfQ2dmJuXPn4uGHH8a4ceOyGoMyHYEqNFsVJu5CFXnqR/0k9JPu57o9qR/1y6VfSRvQWC1/x+nLgnecvuPYHafTp09jxowZeOihhzL+zPXXX4/Gxsb057nnnuvXvmrVKuzYsQOPP/44Xn75ZZw6dQoLFy5ET09PboMJSkdwwdDsXseryFM/6iehn3Q/1+1J/ahflv2mHkMknEuAKfExnJw3TgsWLMA3vvENLFmS+ZB10KBBGD16dPpTVVWVbkulUti0aRPuu+8+zJs3DzNnzsTWrVtx4MABvPDCC7kNJmM6AlUoqipsVDeEVfeZuiGs1I/6Ra2fdD/X7Un9qF+W/bY/kdt7khScvNIRxGKxwKO6p556CiUlJRg2bBg+8YlP4Jvf/CZGjhwJANi1axfmzp2LtrY2VFZWpvvNmDEDixcvxsaNGy8oV5mOQLdyuwtV5Kkf9ZPQT7qf6/akftQvl37JRqDibvmjui/+VO6o7h/nGn1UV/CN0w9/+ENcdNFFGD9+PBoaGvC1r30N3d3d2LdvHwYNGoRt27bhC1/4Arq6uvo9a/78+aipqcE//uM/nienq6ur38+3t7ejuroaSK0DhpQCG+b4oZpxz/86syfu/zfuARt2+2GjvX1fAZ4LVe2JZ9/vaAUwe4X/tW2hnqnqV50CXtnkf21bqGdSP+onMRbak/pRv3D1S3UBVXXcOEVIotAPvPHGG9P/P23aNFxxxRUYP348nn32WeXxnud5iMVigW11dXXB30QFpSPQrdzuQhV56kf9JPST7ue6Pakf9cul395yAHUZ36WhYsH9IwkKvnH6MGPGjMH48eNx+PBhAMDo0aNx9uxZHD9+vN9RXXNzM2bPnh34jDvuuAOrV69O/zn9jVNQOgJlVekRfnipVBX5y1r8f42IVc2mftRPQL/aNr8gKu1J/aifvH57qrJ7+ZLQKPhR3YdpbW3FxRdfjEceeQSf//znkUqlkEwmsXXrVixbtgwA0NjYiHHjxuG5557Dddddd0G5ynQEqqrSRyqAWxYBR4fmFsKqW0W+OgVsfhqobs8thFW3ajb1o34S+h0r8++QtJTRntSP+onr14FISq78+U+BkotCl4ezp4BHzD6qi+fa4dSpU3jzzTfx5ptvAgAaGhrw5ptv4siRIzh16hTWrl2Ln//85/j1r3+N3bt344YbbsCIESPw6U9/GgBQUVGBFStWYM2aNfjpT3+K/fv343Of+xymT5+OefPm5TiYgDBOVVXpt0f6kzQGmSry05v9RRg0ljCqZlM/6iel35kE7Un9qF8U+s2rB4mWnDdOr7/+OmbOnImZM2cCAFavXo2ZM2fir//6rzFgwAAcOHAAixYtwsSJE7F8+XJMnDgRP//5z1FeXp5+xgMPPIDFixdj2bJl+MM//EOUlZXhmWeewYABA3IbTMZ0BLqhoaoQVt023RBW3TbqR/0k9AsjFLyY7Un9qF+W44wqHYEXYx6nPvI6qosKZTqCoq6aTf2on4B+YYSCF7M9qR/1y2WcUaUj+LNdckd13/2kW0d1RlF7nFWzqR/1i0I/3XHSntSP+uWpX0TfyHiCH8NJRD2AvAhKR1DMVbOpH/WT0C+MUPBitif1o365jDPKdAQEgO0bp6B0BMVcNZv6UT8p/XTHSXtSP+qXn36HKqN+8xY9dt9xCkpHUMxVs6kf9ZPQT3ectCf1o37569fcDUxZL3/HacWLcnecNl3DO06hEZSOoJirZlM/6iemn+44aU/qR/3y0u9gEiRa7N44ZUxHEEaotOqZYYR0646F+lE/Cf10x0l7Uj/ql5d+S5cV4OWpgSf4MRy7j+qC0hEUc9Vs6kf9bE05UMz2pH7UL6dnJhBJ5vBbBI/q/pfZR3V2b5xS64AhpayaTf2on6R+uuOkPakf9ctfP6SA3kr5jdMXdsttnDbPMXrjlIh6AHkRlI6gmKtmUz/qJ6GfdIoD1+1J/ahfLs+cXA8sCeOFSrLF7o1TUDoCVTXqhmF6oaiqZ5pUNVu3Krjqmap+NSf0Qm3rq4DORO52of/M8J/uOHXl6equO87LWvxvH6Tmi+7c1V0PtW3A4G76z1b/zWwI5316ITzI3D+SkJEndh/VBaUjUFWj7uzbJ5Z25xaKqnqmUVWz+/TLtSq46pmqfoP79DuTyO2ZyQ7/LkGyIze70H9m+E93nLrydHXXHWd1Ctj8tF/wVWK+6M5d3fVwrMy/W9NSRv/Z6L+u09GUXLl5t9xR3aNmH9XFox5AXgSlI1BVoz483H9ZB7WpQlFVzzSparZuVXDVM1X9Jrb6v3x1xlnanbtd6D9z/KczTl15urrrjnN6s//SlZovunNXdz28WOP/mf6z0387L0UkeDG5j+HYvXHKmI5AFVJqS5tuCKstbbohwba0ue4/19s4P+1uc9h/UaUjIGnsPqoLSkcQRhVr6TbpquDSbdJVyOk/tnF+Fs/8dN5/CUSSjmD5Hrmjui2f4FFdaNQel6tiLd0GyFYFl24DzLE1/cc2zk8z/UD/nd8W75V7x/4unuDHcBJRDyAvgtIRhFHFWrpNuiq4dJt0FXL6j22cn8UzP133H9MRRI7dG6egdARhVbGW7iddFVy6n3QVcvqP/Tg/i2d+uuy/Sa3RvG+lLm5bcDnc7jtOQekIwqhiLd1Puiq4dD/pKuT0H/txfhbP/HTdfyVtQGO1/B2nP3lJ7o7T96/mHafQCEpHEFoVa+l+GUJ0Q6sKLt0vQ6itcX6g/4qzH+enGX6g/87rN/UYIsET/BiO3RunjOkIVKG2qjBVk/rphtqapINuSLAtOhSz/1zvx/lpvg5F6r/tTxTg5Unywe6juqB0BGFUsZbuJ10VXLqfdBVy+o/9OD+LZ3667r9kYzSZw//Lz+SO6n7wcaOP6uzeOKXWAUNKw69iLd1Puiq4dD/pKuT0H/txfhbP/HTdf6kuoKqOG6cISUQ9gLwISkcQRhVr6X7SVcGl+0lXIaf/2I/zs3jmp+v+21sOoK5Qb9Hs8SBz/0hCRp7YvXEKSkcQRtX6+iq/IGOu1a915UlXBdetli5dhZz+M8N/rsvj/KT/TPbfniq5dywJxO6juqB0BGFUrU92+GfmyY7zw0bDkCddFVy3Wrp0FXL6zwz/uS6P85P+M9p/HYik5MpNgkd128w+qotHPYC8CEpHEFbV+tLu3Ktf5yNPuiq4TrV06Srk9J85/nNdHucn/Weq/+bVIxLOJcCU+BiO3RunjOkIdENRdUN0w5Cn2yY9Tl15um30n912sUUe56cZ46T/mI7APOw+qgtKRxBGVfAwQl9NqgoubRfpqu70H+VxfhbP/HTdf1GlI/jMy3JHdY9fxaO60Kg9LlcVHLCjCrn0OHXl6bbRf3bbxRZ5nJ9mjJP+C2gz/yjLdRJRDyAvgtIRhFEVPIzQV5OqgkvbRbqqO/1HeZyfxTM/XfdfZOkIYjL3jyy442T3xikoHUFYVcF1q21LVyGXHqd0FXL6z2672CKP89OMcdJ/AWkMKqN+8xY9dt9xCkpHEEZVcN1q29JVyKXHKV2FnP6z2y62yOP8NGOc9F+wvOZuYMp6+TtOy/5F7o7TE3/IO06hEZSOILSq4AMsqUIuPU5deb2WjJP+K055nJ9mjJP+O0/ewSRItNi9ccqYjsCCkFLjqoJL20VXd/rPbrvYIo/z04xx0n/nPXPpsgK8PHWICeVwikWkX/bYfVQXlI7AlpBSk6qCS9tFuqo7/Ud5nJ/FMz+d918CkWQOX/YKMPCi0OXh/VPAE7N5VBcamdIRABaElF5AnuqZJo1TV56u7vSf3XaxRR7npxnjpP/Of2a875nSeIIfw0lEPYC8CEpHYEtIqUlVwaXtIl3Vnf6jPM7P4pmfrvtvcj2wJMwXK7kQdm+cgtIRqKpRNwyTDSnVrUKueqaqX80J2VBb3SrkqmfSf3b4L1M4uyv6cX7Sf6b6b1JreO9UkhV233EKSkegqkbd2bdPLO2WCSnVrUKueqaq3+Bu///PJGRCbXWrkKueSf+Z7z9VOLsL+nF+0n8m+6+kDWislr/j9Mc/l7vj9E9X8o5TaASlI1BVoz48HEh2yIWU6lYhVz1T1W9iK9BSJhdqq1uFXPVM+s8S/8Ud14/zk/4z1H9Tj4FEi90bp4zpCFShobptuiGltrTphtra0kb/2W1Pzk+72+i/grVtf6IAL08NPMGP4dh9VBeUjsCkqtm2tElXIZduo//stifnpxnzjP4zw3/JRqDibvmjuj8SPKr73zyqC49M6QjinjlVs21pA2SrkEu30X9225Pz0+42+q+AbbFw36uZ4DdOaRJRDyAvgtIRmFQ125Y26Srk0m30n9325Pw0Y57Rf2b4b285gLqo3roEtm+cgtIRmFY1W7cKuXQ/6Srk0v3oP7vtyflpxjyj/6L336HKaN63Xsz/SMgxHLvvOAWlIzCparZuFXLpftJVyKX70X9225Pz04x5Rv+Z4b/mbmDKevk7Tp9+Ve6O046P8Y5TaASlIzCuanbckn66oba69pTuR//ZbU/OTzPmGf0Xuf8OJhEJ575xkvgYjt0bp4zpCFShoaqwWN2QUt1nmtRPN9TWJB3oP3ftyflpvg70n8hYli4rwMuT5IPdR3VB6QhMqpodxlhMqnpukg70n/lV5G3Rj/OT/jPafwkAFfJHdYv3yh3VPTXL6KM6uzdOqXXAkFJgwxw/NDTu+V+f9sT9/8Y9YMNuP4yzt+8rwHOhoz3x7PsdrQBmr/C/Ri3UM03qV50CXtnkfy0d9VjoP/P9J21Pzk/zdaD/5PyHFNBbKb9xWiS4cXra7I1TIuoB5EVQOgKTqmaHMRaTqp6bpAP9Z34VeVv04/yk/0z23+R6YEkhX6QkV+zeOAWlI1BVvz40wg/3lKpifVmL/6+tXEN0a9v8gppS8nSrkOvas77KL6gpJY/+s9uenJ+cn/TfBzrMbAj3vZoJqYvbFlwOt/uoLigdgar69ZEK4JZFwNGhuYWU6laxrk4Bm58GqttzC9E9VuafmbeUycjTrUKua89kh3/nIdkhI4/+s9uenJ+cn/TfBzp0nY6m5Mqn/q/cUd2P/sDoo7p41APIi6B0BKrq12+P9BdFDDJVrKc3+4siaCwXqgp+JiErT6cKeT72LO2WlUf/2W1Pzk/OT/rP12HnpYgET/BjOHZvnDKmI9ANfVWFlOq2hRGiG4Y83TbpcdJ/ZvhP2p6cn2aMk/6L3H9MRxA5dh/VBaUjMKmKdRghuiZVIZceJ/1nhv+k7cn5acY46T9D/JdAJOkIFr4md1T3fz7Ko7rQqD1ueBXrPvNKVtvWlafbJj1O+s8M/0nbk/PTjHHSf9H7L96b/TuShEIi6gHkRVA6ApOqWIcRomtSFXLpcdJ/ZvhP2p6cn2aMk/4zw39RpSPwIHP/SEJGnti9cQpKR2BSFet8qoLrjlO6Crn0OOk/M/wnbU/OTzPGSf9F779JrVG/eYseu+84BaUjMKmKtW5VcN1xSlchlx4n/WeG/6TtyflpxjjpPzP8V9IGNFbL33H6z4J3nJ7lHafwCEpHYFIV67yqguuOU1deryXjpP/M8J+0PTk/zRgn/Re5/6YeQyScS4Ap8TEcuzdOGdMRhBHCqnqmLSHBujqYNE76zwz/SduT89OMcdJ/kftv+xMFeHmSfLD7qC4oHYFJVaxtCbWVrgoubRf6z257cn6aMU76zwz/JRujyRy+4HW5o7p/voJHdaGRKR0BUPjQUNUzVf3inh2htrq6mxQSTP+5a0/OTzPGSf8Z4D/zj7JcJxH1APIiKB2BSVWsbQm1la4KLm0X+s9ue3J+mjFO+s8M/+0tB1AX1ls1M1L3jyy442T3xikoHYGq+nXDML3QUNUzw6gKrjtOXXm6uuuOUzckWLcKueqZqn41J2THqStPV3db7Cmtn7Tfpe2pGyJvy+8X6d+70u+HQ5WFfpOSHLH7jlNQOgJV9evOvn1iaXduoaGqZ4ZRFVx3nLrydHXXHaduSLBuFXLVM1X9BvfZ80xCZpy68nR1t8We0vpJ+13anroh8rb8fpH+vSv9fmjuBqasl7/jdN0+uTtOP7mcd5xCIygdgar69eHhQLIj99BQ1TPDqAquO05debq6645TNyRYtwq56pmqfhNbgZYyuXHqytPV3RZ7Susn7Xdpe+qGyNvy+0X69670++FgEiRa7N44ZUxHoArxZJsZbbohwWxjG9vyb9MNkWdb5G1LlxXg5anB7+ZZCvtjOHYf1QWlIwijSjfbzKh6zja2sS3/Nt0QebaZ0VafAFAhf1Q3/w25o7rnf59HdaGRKR1B3JOt4M0286uQs41tbMs/RJ5t0bfFe6N533qCH8NJRD2AvAhKRxBGlW62mVH1nG1sY1v+bboh8mwzo21yPbAk6pdvcWP3xikoHUFYVbrZz4yq57pVyNmP/dgv/1Qhtvx+cbnfpNao37xFj913nILSEYRRpZv9zKh6rluFnP3Yj/0+6KebKsSW3y+u9ytpAxqr5e84zXsDGFgeujy8fxJ4gXecwiMoHUFoVbrZz4yq53H2Yz/2y7tfhhB5435PsN95/aYeA4kWuzdOGdMRqELdVWG47Gd+1XOTdGA/9rO1X6lBY2G/nPptf6IAL09NPIGPBdh9VBeUjiCMKt3sZ0bVc5N0YD/2s7WfbqoQk3Qo5n7JRqDi7miO6hICR3Xd5h/V5bRxqqurw5NPPolf/epXKC0txezZs3HPPfdg0qRJ6Z/xPA8bN27EI488guPHj2PWrFn49re/jalTp6Z/pqurC2vXrsVjjz2Gzs5OzJ07Fw8//DDGjRuX1TjSjkytA4aUAhvm+KGacc//OrMn7v837gEbdvthuL0x//z/XCh8T5z9oup3tAKYvcI/Voh6LOzHfsXWrzoFvLLJPxaKeizsl3u/VBdQVSe/cZq7X27j9NOZRm+cErn88J49e3Drrbfiox/9KLq7u7F+/XrMnz8fBw8exJAhQwAA9957L+6//348+uijmDhxIr7xjW/g2muvxaFDh1Be7ht91apVeOaZZ/D4449j+PDhWLNmDRYuXIh9+/ZhwIAB2Q8oKB1BGFW62c+Mqucm6cB+7GdrP91UISbpUMz99pYDqMv+PUkKTk4bpx//+Mf9/rx582aMHDkS+/btw9VXXw3P8/Dggw9i/fr1WLJkCQBgy5YtGDVqFLZt24YvfvGLSKVS2LRpE77//e9j3rx5AICtW7eiuroaL7zwAq677rrsBxSUjkBZVXqEH66ba4h8fZVfADLnKtaUV9Aq8pe1+P9azjXEurbNL4hKeZRHecDVv/HXn6u/X1yXt6cql9d24fAgcwdJQkae5HXH6Z133sGECRNw4MABTJs2De+++y5qa2vxxhtvYObMmemfW7RoEYYNG4YtW7Zg165dmDt3Ltra2lBZWZn+mRkzZmDx4sXYuHHjeXK6urrQ1dWV/nN7ezuqq6uD0xGoqkofqQBuWQQcHZpbiHyyw78TkOw4P2yU8nKXp1tFvjoFbH4aqG7PLcT6WJl/p6qljPIoj/IGdfs/15Vw8/eL6/K6OhBJyZVP7pc7qttl9lFdXLej53lYvXo1rrrqKkybNg0A0NTUBAAYNWpUv58dNWpUuq2pqQklJSX9Nk0f/pkPU1dXh4qKivSnurrabwhKR6CqKv32SH+SxpB71fPS7tyrWFNeYavIT2/2XxJBY7lQFfIzCcqjPMo7J68r4fbvF5flzatHJHgxuY/haG+cVq5cibfeeguPPfbYeW2xWH/FPc877+8+jOpn7rjjDqRSqfTn6NGjfkPGdAS6oba6IfKUZ35Vd1vsQnmUJyFPt811u1ggL8p0BASA5lHdbbfdhqeeegovvfQSampq0n8f1lHdh1GmI9ANtQ2jmjjlmVHV3Ra7UB7lScjTbXPdLrbIiyodwTVvyh3Vvfgf3Tmq8zwPK1euxJNPPoldu3b12zQBQE1NDUaPHo2dO3em/+7s2bPYs2cPZs+eDQC4/PLLMXDgwH4/09jYiLfffjv9M1lTezz3itOAbDVxyjOjqrstdqE8ypOQp9vmul2skBfRUZYn+DGcRC4/fOutt2Lbtm14+umnUV5enr6TVFFRgdLSUsRiMaxatQp33XUXJkyYgAkTJuCuu+5CWVkZbrrppvTPrlixAmvWrMHw4cNRVVWFtWvXYvr06ekou6wJSkegG2obRjVxyjOjqrstdqE8ypOQp9vmul1skcd0BJGT08bpO9/5DgBgzpw5/f5+8+bNuPnmmwEAt99+Ozo7O/GVr3wlnQDz+eefT+dwAoAHHngAiUQCy5YtSyfAfPTRR3PL4QQEpyO4UDXqTCGeYVUTpzwzqrrbYhfKozwJeYleO8ZJeQFpDCpze08WCqmL2xZcDre75EpQOgJVVWlViGcY1cQpz4yq7rbYhfIoT0JezQngmW1+tJfJ46S8YHnN3cCU9fJ3nD7x/+TuOO2Z4c4dJ+MISkdwwWrUPcLVxCnPjKruttiF8ihPQl6vJeOkvPPkHUwiEjzBj+HYvXHKmI7AgpDSopeneqZJ46Q8ynNRnu7adN0uFshbuqwAL0+SD3Yf1QWlI7AlpLSY5UlXdbfFLpRHeRLydNem63axRV59ApFkDv/4W3JHdT/7DzyqC41M6QgAC0JKi1ye6pkmjZPyKM9Febpr03W72CAv3vdMEhmJqAeQF0HpCGwJKS1medJV3W2xC+VRnoQ83bXpul1skTe5HlgS9ss1AA8y948kZOSJ3RunoHQEqmrUDcNkQ0pfGu+PM9eq56pwYVf0y/RMVb+aE7KhxLpV5FXpFqgf9TNBv0zPdOX3i8u/Pye1hvdOJVlh9x2noHQEqmrUnX37xNJumZDSrj55g7rPb1NVPVeFC7ugn+qZqn6Du/3/P5OQCSXWrSKvSrdA/ahf1PqpnunC7xfXf3+WtAGN1fJ3nK56S+6O08u84xQeQekIVNWoDw8Hkh1yIaXXNPiLIteq56pwYWf0y/BMVb+JrUBLmVwosW4VeVW6BepH/YzQL8Mznfn94vDvz6nHEAnnEmBKfAzH7o1TxnQEqtBQ3TbdkFLdNt3QV1v0k25z3Z7Uj/qZ3Oa6PQX12/5EAV6eJB/sPqoLSkdgUlVw3TZbqp7b0ua6Pakf9TO5zXV7SuuXbAQq7pY/qpt9QO6o7pXpPKoLjUzpCOKeOVXBddsAO6qe29Lmuj2pH/Uzuc11e4rqZ/5Rluskoh5AXgSlIzCpKrhumy1Vz21pc92e1I/6mdzmuj2l9dtbDqBO/n0rdf/IgjtOdm+cgtIRmFYVPFNYsytVz3X1k+7nuj2pH/UzuZ/r9pTU71Bl1G/eosfuO05B6QhMqgquCmt2oeq5rn7S/Vy3J/Wjfib3c92e0vo1dwNT1svfcfrY23J3nF6dxjtOoRGUjsC4quBxx6ue6+on3c91e1I/6mdyP9ftKajfwSRItNi9ccqYjkAVGqoKU9UNKdV9pm7IrAv6Sfdz3Z7Uj/qZ3M91ewrqt3RZAV6eGpy74yTxMRy7j+qC0hGYVBVcdyy2VD0PYyy0J/WjfvatsWK2p/jvzwSACvmjulm/kDuq2zvV6KM6uzdOqXXAkFJgwxw/NDTu+V+f9sT9/8Y9YMNuP4yzt28ney50tCeefb+jFcDsFf7XqIV6pqpfdQp4ZZP/NXGhnmmSftL9XLcn9aN+Jvdz3Z7S+iEF9FbKb5z+QHDj9H/N3jgloh5AXgSlIzCpKrjuWGypeh7GWGhP6kf97FtjxWxPaf0m1wNLCv0yJblg98YpKB2Bqvr1oRF+uKdUFevLWvx//UhVPbdFv9o2vyAq7Un9qJ+8fvVVfkFb2tNO/WY2FPY9SnLG7qO6oHQEqurXRyqAWxYBR4fmFlKqW8W6OgVsfhqobs8tZFa36rkt+h0r8+8EtJTRntSP+knrl+zw7wAlO2hPG/XrOh1NyZWPHpQ7qnttitFHdfGoB5AXQekIVNWv3x7pL4oYZKpYT2/2F0XQWMKoem6TfmcStCf1o35R6VfaTXvaqt/OS0Gixe6NU8Z0BLqhqKqQUt023ZBZ3TZb9AsjtLeY7Un9qJ/J69Z1ewrqF1k6AsGP4dh9VBeUjsCkKtauV+kOIySY9qR+1M/MdUR7mqFfVOkIrhA8qnudR3XhUXvc8CrWfeZ1tkp3HvrpjpP2pH7UL7p1RHtGr1+8V/+dmQ/8xilNIuoB5EVQOgKTqli7XqU7jJBg2pP6UT8z1xHtaYZ+TEcQOXZvnILSEZhUxdr1Kt356Kc7TtqT+lG/6NYR7Rm9fpNaI3rhxiBTDkVCRn7YfccpKB2BSVWsXa/Srauf7jhpT+pH/fLXT3rdum5Paf1K2oDGavk7Tpf/EhggcMep5ySwbzLvOIVGUDoCk6pYu16lOy/9dMdJe1I/6hfdOqI9I9dv6jFEgheT+xiO3RunjOkIwgjRVT0zjNBe3bHYol+BQ3SL3p7Uj/qZvG5dt6egftufKMDLk+SD3Ud1QekITKpi7XqVbltCpV23J/WjfiavW9ftKa1fsjGazOG//yu5o7o3LjP6qM7ujVNqHTCk1Nwq1q5X6dbVT3ectCf1o37RrSPa0wz9Ul1AVZ38xmmm4MZpv9kbp0TUA8iLoHQEJlWxdr1Kty2h0q7bk/pRP5PXrev2lNZvbzmAujDeqCRL7N44BaUjUFW/bhimFxqqeqYLVbpVz1T1qzmhF9qrO05debq621LVXXfu6q6H2jZgcDf9Z6v/dMfJ359m/P7cU1X4d2k2SF3ctuByuN1HdUHpCFTVrzv79oml3bmFhqqe6UKVbtUzVf0G9+l3JiEzTl15urrbUtVdd+7qrodjZf5dkJYy+s9G/+mOk78/Dfn92YFISq78x0NyR3VvTjL6qC4e9QDyIigdgar69eHhQLIj99BQ1TNdqNKteqaq38RW/+UpNU5debq621LVXXfu6q6HF2v8P9N/9vpPZ5z8/WnG78959YgET/BjOHZvnDKmI1CFeNrSphsyyzYz2nRDkG1pc31+uu4/19scnp9MRxA5dh/VBaUjCKNqtnSbdJVutplR9dyWNtfnp+v+c73N9fkZVTqC//Cvckd1b03kUV1o1B6Xq5ot3QbIVulmW/RVz21qc31+uu4/19ucnp8x2fcsOY9E1APIi6B0BGFUzZZuk67SzTYzqp7b0ub6/HTdf663uT4/o0pH4EHm/pGEjDyxe+MUlI4grKrZ0v2kq3SzX/RVzzk/zennuv9c7+fy/DxUGfWbt+ix+45TUDqCMKpmS/eTrtLNfmZUPef8NKOf6/5zvZ/r87O5G5iyXv6O03TBO04HeMcpPILSEYRWNVu6X6ZQ1HdDqtLNfpFXPef8NKif6/5zvZ/D8/NgEpFwLgGmxMdw7N44ZUxHoApFVYUZm9RPFaZqiw7F3E83VNokHYp5frruP9f7OTw/ly4rwMuT5IPdR3VB6QjCqJot3U+6Sjf7mVH13CQdinl+uu4/1/u5Pj/rE4gkc/jUw3JHdb+YYPRRnd0bp9Q6YEhp+FWzpftJV+lmv8L2k67qzvlJ/7Ff8cxPpIDeSm6cIiQR9QDyIigdQRhVs6X7SVfpZj8zqp6bpEMxz0/X/ed6P9fn5+R6YEmhXqI5IHX/yII7TnZvnILSEYRRNbu+yi9waXqVbt2q7rrV7inPjKrunJ/0H+UVz/yc2SD7niXnYfdRXVA6gjCqZic7/DsPyY7zw2JNqtKtW9Vdt9o95ZlR1Z3zk/6jvOKZn12noym5MvkduaO6X37E6KO6eNQDyIugdARhVc0u7Ta/Snc+Vd11qt1TnjlV3Tk/6T/KK475ufNSkGixe+OUMR2BbiixbghrGPJ026THSXlmVHXn/KT/KK8o5mdU6QjO3XGS+BiO3Ud1QekIwqiaHUbosklV3W2xiy3ypKu6c37Sf5RXPPMzqnQEl9XLHdX9qpZHdaFRe1yuajZgR5Vu6XFSnhlV3Tk/6T/KK475Ge/N/E4ME0/wYziJqAeQF0HpCMKomh1G6LJJVd1tsYst8qSrunN+0n+UVzzzM6p0BCSN3RunoHQEYVXN1q2WLl2lW3qclGdGVXfOT/qP8opjfk5qjfrNW/TYfccpKB1BGFWzdaulS1fplh4n5ZlR1Z3zk/6jvOKZnyVtQGO1/B2nie/K3XH610t5xyk0gtIRhFY1e4AlVbqlx0l5ZlR15/yk/yivKObn1GMg0WL3xiljOgJbQmZV8lTPNGmclGdGVXfOT/qP8opifm5/ogAvTw08wY/h2H1UF5SOwJaQWZOquttiF1vkSVd15/yk/yiveOZnsjGazOET3pU7qjvMo7rwyJSOALA/FFX1TJPGSXmF9a0tdnF9frruP9flOT0/Y4iEQiS2zPZjOImoB5AXQekIbAmZNamquy12sUWedFV3zk/6j/KKZ37uLQdQF+ablVwAuzdOQekIVNW2G4bJhqLqVulWPVPVr+aEbMisbhV5VbiwK/pleibnJ/1ngv8ypVtwRT+X5+ehyvDeqSo8yNw/kpCRJ3bfcQpKR6Cqtt3Zt08s7ZYJRdWt0q16pqrf4G7//88kZEJmdavIq8KFXdBP9UzOT/ovav+p0i24oJ/r87O5G5iyXv6OU22D3B2n+hrecQqNoHQEqmrbh4cDyQ65UFTdKt2qZ6r6TWwFWsrkQmZ1q8irwoWd0S/DMzk/6T8j/NfruH4Oz8+DSUSCB6E7TtGolwt2b5wypiNQhXjqtumGotrSphsya4s9pfWTbuP8tLvN9fXA+VmwtqXLCvDyJPlg91FdUDoCk6qe29JmSxV5W/STbuP8NMMPtvjP9fXu+vysTwCokD+qu7QBGBC+PPS0A+/yqC48MqUjiHvmVD23pQ2wo4q8LfpJt3F+2t3m+nrg/CxcW7w3vHfqhfAEPhaQiHoAeRGUjsCkque2tNlSRd4W/aTbOD/N8IMt/nN9vbs+PyfXA0uieukSwPaNU1A6AtOqnutW6ZbuZ0sVeVv0k+7H+WmGH2zxn+vr3eX5Oak1mvftucvbEnIMx+47TkHpCEyqeq5bpVu6ny1V5G3RT7of56cZfrDFf66vd9fnZ0kb0Fgtf8ep5tdAPHx56G0HGn6Pd5xCIygdgXFVz+OW9NMNJbbFntL6Sffj/DTDD7b4z/X17vD8nHoMkeAJfgzH7o1TxnQEqtBXVdiobsis7jNN6qcKi3XBntL6Sffj/DRfB5P85/p6d3h+bn+iAC9Pkg92H9UFpSMwqep5GGMxqSq4LfaU1k+6H+en+TqY5D/X17vr8zPZCFTcLX9UN/43ckd1vxlv9FGd3Run1DpgSCmwYY4fqhn3/K8ze+L+f+MesGG3Hzba23ex7VxobE88+35HK4DZK/yvbQv1TJP6VaeAVzb5XxNLjEXantL6Sffj/DRfB5P85/p6d31+prqAqjpunCIkEfUA8iIoHYFJVc/DGItJVcFtsae0ftL9OD/N18Ek/7m+3l2fn3vLAdQV8k2aHR5k7h9JyMgTuzdOQekIlFWlR/jhpVJVui9r8f81kmsIa22bX3BSSp52lW5L7CmtX32VXxDVVXu6Pj9d95/r6931+bmnKrx3KsmKnI7q6urq8OSTT+JXv/oVSktLMXv2bNxzzz2YNGlS+mduvvlmbNmypV+/WbNm4dVXX03/uaurC2vXrsVjjz2Gzs5OzJ07Fw8//DDGjRuX1TiU6QhUVaWPVAC3LAKODs0tZFa3Snd1Ctj8NFDdnlsI67Ey/05AS5mMPN0q3bbYU1q/ZId/5yHZ4aY9XZ+frvvP9fXu+vzs6kAkJVcu+Y3cUd0Rs4/q4rn88J49e3Drrbfi1Vdfxc6dO9Hd3Y358+fj9OnT/X7u+uuvR2NjY/rz3HPP9WtftWoVduzYgccffxwvv/wyTp06hYULF6Knpye30QelI1BVlX57pD9JY5Cp0j292V+EQWO5UNXsMwlZeTpVum2yp7R+pd1u29P1+em6/1xf7y7Pz3n1iIRzCTAlPoaT08bpxz/+MW6++WZMnToVM2bMwObNm3HkyBHs27ev388NGjQIo0ePTn+qqj74ajGVSmHTpk247777MG/ePMycORNbt27FgQMH8MILL+Q2+ozpCHRDQ1Uhs7ptYYSwhiFPt80We0rrJ+13zk8zxmmL/1xf7w7PT6YjiJy8oureeecdTJgwAQcOHMC0adMA+Ed1Tz31FEpKSjBs2DB84hOfwDe/+U2MHDkSALBr1y7MnTsXbW1tqKysTD9rxowZWLx4MTZu3HhBucp0BCZV6Q4jhNWoKt2W2FNaP2m/c36aMU5b/Of6end9fkaVjmDcEbmjuvcuceeo7nfxPA+rV6/GVVddld40AcCCBQvwgx/8ALt27cJ9992H1157DZ/85CfR1dUFAGhqakJJSUm/TRMAjBo1Ck1NTYGyurq60N7e3u8DAKg9bn6VbkC2mriuPN02W+wprZ+03zk/zRinLf5zfb07PT/NP8pynYRux5UrV+Ktt97Cyy+/3O/vb7zxxvT/T5s2DVdccQXGjx+PZ599FkuWLMn4PM/zEIsFT4i6urrgb6KC0hGYVKU7jBBWk6p022JPaf2k/c75acY4bfGf6+vd9fkZWTqCmMz9IwvuOGltnG677Tb86Ec/wksvvXTBSLgxY8Zg/PjxOHz4MABg9OjROHv2LI4fP97vW6fm5mbMnj078Bl33HEHVq9enf5ze3s7qqurg9MRmFalW7dqtu44Xa96bot+0n7n/DRjnLb4z/X17vL8PFQZ+J4kcuR0x8nzPNx2223YsWMHdu/ejQkTJlywT2trKy6++GI88sgj+PznP49UKoVkMomtW7di2bJlAIDGxkaMGzcOzz33HK677roLPlOZjsCkKt26VbN1x+l61XNb9JP2O+enGeO0xX+ur3fX52dzNzBlvfwdp7FH5e44/Vu1O3ecbr31VmzduhXbtm1DeXk5mpqa0NTUhM7OTgDAqVOnsHbtWvz85z/Hr3/9a+zevRs33HADRowYgU9/+tMAgIqKCqxYsQJr1qzBT3/6U+zfvx+f+9znMH36dMybNy+30QelIzCuSvcA4WriuvJ6Lal6bot+0n7n/DRjnLb4z/X17vD8PJgEiRgvB/BB0vV+n82bN3ue53kdHR3e/PnzvWQy6Q0cONC75JJLvOXLl3tHjhzp95zOzk5v5cqVXlVVlVdaWuotXLjwvJ9RkUqlfNn/9IiHgd0evr7LQ/dGDwv+1UO81wM8D9cf9nD26x6+ttvDgB4Pz33fw2/u9zC2PbitZ4MH704P//z93J6p6nfJCQ+//dvcn6k7Tl15urrbYk9p/aT9zvlpxjht8Z/r6931+Rk74QHwUqlULq9vbdLv2zHvebi4PfzPmPdE9dPB7iK/QekITKrSbVLIrAtVz23RT9rvnJ9mjNMW/7m+3l2fn/UJRJI5fMx7ckd1jeOMPqqze+OUWgcMKTW3Srdu1Wzdcbpe9dwW/aT9zvlpxjht8Z/r6931+YkU0Fspv3EaLbhxajJ745SIegB5EZSOwKQq3SaFzLpQ9dwW/aT9zvlpxjht8Z/r6931+Tm5HlgS1kuVZIPdG6e2UmDnpf6OfEAvcHkjMLXZr4yd7ACufdf/+564X57lXJXuoLZfjPQrorfk+ExVv3NVs3N9pu44deXp6m6LPaX1k/Y756cZ47TFf66vd9fn55UNub0nScGx8qgulUph2LBhOHr0qLFf5RFCCCGF5lwewxMnTvhHaALyKioqgFG/lTuq+/eLeVRXaE6ePAkAfhJMQgghpMg4efKkyMaJnI+VG6exY8fi4MGDmDJlCr91Colz/6qhfQsPbRsutG940Lbhko19Pc/DyZMnMXbsWNnBeX0fCTmGY+XGKR6P4+KLLwYADB06lAs4RGjf8KBtw4X2DQ/aNlwuZF9+0xQtVm6cCCGEECIIi/ymiUc9AEIIIYQQW7D2G6dBgwbhzjvvxKBBg6IeipPQvuFB24YL7RsetG24GG1f3nFKY2U6AkIIIYSETzodwYh/k0tH0DKW6QgIIYQQYjG845SGd5wIIYQQQrKEGydCCCGEkCyxduP08MMPo6amBoMHD8bll1+On/3sZ1EPyTo2bNiAWCzW7zN69Oh0u+d52LBhA8aOHYvS0lLMmTMHv/jFLyIcsbm89NJLuOGGGzB27FjEYjE89dRT/dqzsWVXVxduu+02jBgxAkOGDMGnPvUpvPfee4JamMuF7HvzzTefN5c/9rGP9fsZ2jeYuro6fPSjH0V5eTlGjhyJxYsX49ChQ/1+hvNXj2xsa83c9QQ/hmPlxumHP/whVq1ahfXr12P//v34+Mc/jgULFuDIkSNRD806pk6disbGxvTnwIED6bZ7770X999/Px566CG89tprGD16NK699tp0yRvyAadPn8aMGTPw0EMPBbZnY8tVq1Zhx44dePzxx/Hyyy/j1KlTWLhwIXp6eqTUMJYL2RcArr/++n5z+bnnnuvXTvsGs2fPHtx666149dVXsXPnTnR3d2P+/Pk4ffp0+mc4f/XIxrYA5651eBbyB3/wB96XvvSlfn932WWXeevWrYtoRHZy5513ejNmzAhs6+3t9UaPHu3dfffd6b87c+aMV1FR4f3DP/yD0AjtBIC3Y8eO9J+zseWJEye8gQMHeo8//nj6Z37729968Xjc+/GPfyw2dhv4sH09z/OWL1/uLVq0KGMf2jd7mpubPQDenj17PM/j/C0kH7at55k/d1OplP89UGWTh+Ed4X8qmzwAXiqVCl03Xaz7xuns2bPYt28f5s+f3+/v58+fj1deeSWiUdnL4cOHMXbsWNTU1OAzn/kM3n33XQBAQ0MDmpqa+tl50KBB+MQnPkE750g2tty3bx/ef//9fj8zduxYTJs2jfbOkt27d2PkyJGYOHEi/uzP/gzNzc3pNto3e1KpFACgqqoKAOdvIfmwbc/BuWsX1m2cWlpa0NPTg1GjRvX7+1GjRqGpqSmiUdnJrFmz8L3vfQ8/+clP8N3vfhdNTU2YPXs2Wltb07aknfMnG1s2NTWhpKQElZWVGX+GZGbBggX4wQ9+gF27duG+++7Da6+9hk9+8pPo6uoCQPtmi+d5WL16Na666ipMmzYNAOdvoQiyLWDZ3OX9JgAW53GKxfrnevA877y/I2oWLFiQ/v/p06fjyiuvRG1tLbZs2ZK+nEg7Fw4dW9Le2XHjjTem/3/atGm44oorMH78eDz77LNYsmRJxn60b39WrlyJt956Cy+//PJ5bZy/+ZHJtpy79mHdN04jRozAgAEDzttpNzc3n/cvIpIbQ4YMwfTp03H48OF0dB3tnD/Z2HL06NE4e/Ysjh8/nvFnSPaMGTMG48ePx+HDhwHQvtlw22234Uc/+hFefPFFjBs3Lv33nL/5k8m2QRg7d88lwJT4GI51G6eSkhJcfvnl2LlzZ7+/37lzJ2bPnh3RqNygq6sLv/zlLzFmzBjU1NRg9OjR/ex89uxZ7Nmzh3bOkWxsefnll2PgwIH9fqaxsRFvv/027a1Ba2srjh49ijFjxgCgfVV4noeVK1fiySefxK5du1BTU9OvnfNXnwvZNgjOXQuI5k56fjz++OPewIEDvU2bNnkHDx70Vq1a5Q0ZMsT79a9/HfXQrGLNmjXe7t27vXfffdd79dVXvYULF3rl5eVpO959991eRUWF9+STT3oHDhzwPvvZz3pjxozx2tvbIx65eZw8edLbv3+/t3//fg+Ad//993v79+/3fvOb33iel50tv/SlL3njxo3zXnjhBe+NN97wPvnJT3ozZszwuru7o1LLGFT2PXnypLdmzRrvlVde8RoaGrwXX3zRu/LKK72LL76Y9s2CL3/5y15FRYW3e/dur7GxMf3p6OhI/wznrx4Xsq0NczcdVVfx7x6GdYb/qfh346PqrNw4eZ7nffvb3/bGjx/vlZSUeL//+7/fL7yTZMeNN97ojRkzxhs4cKA3duxYb8mSJd4vfvGLdHtvb6935513eqNHj/YGDRrkXX311d6BAwciHLG5vPjii4FXHZcvX+55Xna27Ozs9FauXOlVVVV5paWl3sKFC70jR45EoI15qOzb0dHhzZ8/30smk97AgQO9Sy65xFu+fPl5tqN9gwmyKwBv8+bN6Z/h/NXjQra1Ye5y43Q+Mc/zLLnHTgghhBBJ2tvbUVFRAVT8OxAbGr5Arx1IjUIqlcLQoQLyNLA2qo4QQgghQngxAAIXt3k5nBBCCCHEHfiNEyGEEELUSF3qseDyEL9xIoQQQgjJEn7jRAghhBA1vOOUht84EUIIIYRkCb9xIoQQQoga3nFKw2+cCCGEEEKyhN84EUIIIUQN7zil4TdOhBBCCCFZwm+cCCGEEKKGd5zS8BsnQgghhJAs4caJEEIIISRLeFRHCCGEEDW8HJ6G3zgRQgghhGQJv3EihBBCiBpeDk/Db5wIIYQQQrKE3zgRQgghRA3vOKXhN06EEEIIIVnCb5wIIYQQooZ3nNLwGydCCCGEkCzhN06EEEIIUcM7Tmn4jRMhhBBCSJZw40QIIYQQkiU8qiOEEEKIGl4OT8NvnAghhBBCsoTfOBFCCCFEDS+Hp+E3ToQQQgghWcJvnAghhBCihnec0vAbJ0IIIYSQLOHGiRBCCCFqvJjcJ0S++c1vYvbs2SgrK8OwYcO0nsGNEyGEEEKKgrNnz2Lp0qX48pe/rP0M3nEihBBCyIWx4P7Rhdi4cSMA4NFHH9V+Br9xIoQQQgjJEn7jRAghhJAL0C4qp729v7xBgwZh0KBBQmNQw2+cCCGEEBJISUkJRo8eDaAaQIXApxoXXXQRqqurUVFRkf7U1dVlHOOGDRsQi8WUn9dff71gNuE3ToQQQggJZPDgwWhoaMDZs2fFZHqeh1isf3Sd6tumlStX4jOf+Yzymb/3e79XiKEB4MaJEEIIIQoGDx6MwYMHRz2MjIwYMQIjRowQk8eNEyGEEEKKgiNHjqCtrQ1HjhxBT08P3nzzTQDARz7yEVx00UVZPSPmeZ4DAYaEEEIIIWpuvvlmbNmy5by/f/HFFzFnzpysnsGNEyGEEEJIljCqjhBCCCEkS7hxIoQQQgjJEm6cCCGEEEKyhBsnQgghhJAs4caJEEIIISRLuHEihBBCCMkSbpwIIYQQQrKEGydCCCGEkCzhxokQQgghJEu4cSKEEEIIyRJunAghhBBCsoQbJ0IIIYSQLPn/4QF+UKOwvZYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 700x700 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "H = linalg.hadamard(n=2**8)\n",
    "\n",
    "fig = plt.figure(figsize = (7,7))\n",
    "ax = fig.add_axes([.1,.1,.8,.8])\n",
    "\n",
    "im = ax.imshow(H, cmap='winter')\n",
    "fig.colorbar(im, ticks=[-1,1])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d462820d",
   "metadata": {},
   "source": [
    "## Sparse Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f56421",
   "metadata": {},
   "source": [
    "A sparse matrix (loosely defined) is a matrix with a many zero elements. When dealing with large sparse matrices it is extrememly useful to take advantage of their sparseness and the memory requirements, speeding-up our code. To convert from a memory dense matrix to a memory sparse matrix we use `scipy.sparse.csr_matrix(A)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "78ba43e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0 1]\n",
      " [2 0 0]\n",
      " [0 3 0]] \n",
      "\n",
      "<class 'scipy.sparse._csr.csr_matrix'> \n",
      "\n",
      "  (0, 2)\t1\n",
      "  (1, 0)\t2\n",
      "  (2, 1)\t3\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[0,0,1],[2,0,0],[0,3,0]])\n",
    "print(A,'\\n')\n",
    "\n",
    "A = scipy.sparse.csr_matrix(A)\n",
    "print(type(A),'\\n')\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e26518b",
   "metadata": {},
   "source": [
    "Notice that each tuple corrseonds to the element location. We can conceivably represent any matrix in this form, but clearly this would be less efficient for non-sparse matrices. We can perform the standard matrix operations with sparse matrices as shown below. To convert back to a dense matrix, we use the `A.todense()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "15937567",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 0 1]\n",
      " [2 0 0]\n",
      " [0 3 0]] \n",
      "\n",
      "[[0 0 2]\n",
      " [4 0 0]\n",
      " [0 6 0]] \n",
      "\n",
      "[[0 3 0]\n",
      " [0 0 2]\n",
      " [6 0 0]] \n",
      "\n",
      "[[0 3 0]\n",
      " [0 0 2]\n",
      " [6 0 0]] \n",
      "\n",
      "  (0, 2)\t1\n",
      "  (1, 0)\t2\n",
      "  (2, 1)\t3 \n",
      "\n",
      "[[3]\n",
      " [2]\n",
      " [6]]\n"
     ]
    }
   ],
   "source": [
    "print(A.todense(),'\\n')\n",
    "\n",
    "print((A+A).todense(),'\\n')\n",
    "\n",
    "print((A*A).todense(),'\\n')\n",
    "\n",
    "print((A.dot(A)).todense(),'\\n')\n",
    "\n",
    "print(A, '\\n')\n",
    "\n",
    "vec = np.array([[1],[2],[3]])\n",
    "\n",
    "print((A*vec))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63ab4df2",
   "metadata": {},
   "source": [
    "It is also useful to initialise an empty sparse matrix and then populate it with non-zero elements with `scipy.sparse.lil_matrix((n,m))`. This saves us from creating a very large matrix with many zeros in the first place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "f9af87d1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 0)\t1.0\n",
      "  (0, 1)\t2.0\n",
      "  (2, 0)\t7.0\n",
      "  (2, 2)\t1.0 \n",
      "\n",
      "[[1. 2. 0.]\n",
      " [0. 0. 0.]\n",
      " [7. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "A = scipy.sparse.lil_matrix((3,3)) #Empty 3x3 sparse matrix\n",
    "\n",
    "A[0,0] = 1\n",
    "A[0,1] = 2\n",
    "A[2,0] = 7\n",
    "A[2,2] = 1\n",
    "\n",
    "\n",
    "print(A,'\\n')\n",
    "print(A.todense())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3023f511",
   "metadata": {},
   "source": [
    "## Matrix Exponential"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b40a37ae",
   "metadata": {},
   "source": [
    "The matrix exponential of a matrix $A$ is written as\n",
    "\n",
    "$$ \\exp (A) = e^{A} = \\sum_{k=0}^{\\infty} \\frac{1}{k!}A^k . $$\n",
    "\n",
    "This seems like a monumental task but it can easily be achieved with `scipy.linalg.expm(A)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "6e937f1d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5  0.5  0.5  0.5]\n",
      " [ 0.5 -0.5  0.5 -0.5]\n",
      " [ 0.5  0.5 -0.5 -0.5]\n",
      " [ 0.5 -0.5 -0.5  0.5]] \n",
      "\n",
      "[[ 2.13068123  0.5876006   0.5876006   0.5876006 ]\n",
      " [ 0.5876006   0.95548004  0.5876006  -0.5876006 ]\n",
      " [ 0.5876006   0.5876006   0.95548004 -0.5876006 ]\n",
      " [ 0.5876006  -0.5876006  -0.5876006   2.13068123]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "n=4\n",
    "A = (1/np.sqrt(n))*linalg.hadamard(n)\n",
    "print(A,'\\n')\n",
    "\n",
    "expA = linalg.expm(A)\n",
    "\n",
    "print(expA,'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98215844",
   "metadata": {},
   "source": [
    "Also we can verify that \n",
    "$$ \\exp \\left(A^T\\right) = \\exp \\left(A\\right)^T $$\n",
    "and\n",
    "$$ \\exp \\left(A^\\dagger \\right) = \\exp \\left(A \\right)^\\dagger $$\n",
    "up to a computational error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "c65b2301",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ True  True  True  True]\n",
      " [ True  True  True  True]\n",
      " [ True  True  True  True]\n",
      " [ True  True  True  True]] \n",
      "\n",
      "[[ True  True  True  True]\n",
      " [ True  True  True  True]\n",
      " [ True  True  True  True]\n",
      " [ True  True  True  True]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = rnd.random((n,n)) + 1j*rnd.random((n,n))\n",
    "\n",
    "Q = np.round(linalg.expm(A).T,3) == np.round(linalg.expm(A.T),3)\n",
    "print(Q,'\\n')\n",
    "\n",
    "Q = np.round((np.matrix(linalg.expm(A))).H,3) == np.round(linalg.expm(np.matrix(A).H),3)\n",
    "print(Q,'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00268ad6",
   "metadata": {},
   "source": [
    "One key point with matrices is that if $A$ and $B$ do not commute,\n",
    "$$AB \\neq BA \\iff [A,B]= AB -BA \\neq0.$$\n",
    "then\n",
    "$$ e^A e^B \\neq e^{A+B}.$$\n",
    "\n",
    "But if \n",
    "$$[A,B]= AB -BA =0$$\n",
    "then\n",
    "$$e^Ae^B = e^{A+B}. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "9f048e18",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ True  True]\n",
      " [ True  True]] \n",
      "\n",
      "[[False False]\n",
      " [False False]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[0,1],[0,0]])\n",
    "B = np.eye(2)\n",
    "# These matrices commute\n",
    "\n",
    "EA = linalg.expm(A)@linalg.expm(B) == linalg.expm(A+B)\n",
    "print(EA,'\\n')\n",
    "\n",
    "A = np.array([[0,1],[0,0]])\n",
    "B = np.array([[0,1],[1,0]])\n",
    "# These matrices do not commute\n",
    "\n",
    "EA = linalg.expm(A)@linalg.expm(B) == linalg.expm(A+B)\n",
    "print(EA,'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23a839ed",
   "metadata": {},
   "source": [
    "## Unitary Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "405bba2d",
   "metadata": {},
   "source": [
    "A Unitary Matrix is defined as a matrix that satisfies the following equation,\n",
    "$$AA^{\\dagger} =  A^{\\dagger}A = \\mathbb{1}$$\n",
    "where $\\mathbb{1}$ is the identity matrix. In other words the Hermitian conjugate of the matrix is it's multiplicative inverse,\n",
    "$$A^{-1} = A^\\dagger$$\n",
    "The Hadamard gate is unitary as we see below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "b3090022",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5  0.5  0.5  0.5]\n",
      " [ 0.5 -0.5  0.5 -0.5]\n",
      " [ 0.5  0.5 -0.5 -0.5]\n",
      " [ 0.5 -0.5 -0.5  0.5]] \n",
      "\n",
      "[[ True  True  True  True]\n",
      " [ True  True  True  True]\n",
      " [ True  True  True  True]\n",
      " [ True  True  True  True]] \n",
      "\n",
      "[[1. 0. 0. 0.]\n",
      " [0. 1. 0. 0.]\n",
      " [0. 0. 1. 0.]\n",
      " [0. 0. 0. 1.]] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "n=4\n",
    "A = (1/np.sqrt(n))*linalg.hadamard(n)\n",
    "print(A,'\\n')\n",
    "A = np.matrix(A)\n",
    "\n",
    "print(A*A.H==A.H*A,'\\n')\n",
    "\n",
    "print(A*A.H,'\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c86a55b",
   "metadata": {},
   "source": [
    "We can also generate a random unitary matrix, of dimension `n`, using `scipy.stats.unitary_group.rvs(n)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "012519a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.-0.j  0.+0.j  0.+0.j -0.+0.j]\n",
      " [ 0.-0.j  1.-0.j -0.-0.j -0.-0.j]\n",
      " [ 0.-0.j -0.+0.j  1.+0.j  0.-0.j]\n",
      " [-0.-0.j -0.-0.j  0.+0.j  1.+0.j]] \n",
      "\n",
      "[[ 0.084-0.002j -0.39 +0.416j  0.013+0.118j -0.779+0.215j]\n",
      " [-0.546+0.053j  0.525+0.222j  0.39 -0.018j -0.3  -0.364j]\n",
      " [-0.099+0.489j -0.31 -0.098j  0.576+0.501j  0.212+0.13j ]\n",
      " [-0.341-0.572j -0.311-0.384j  0.409-0.287j -0.054+0.244j]] \n",
      "\n",
      "[[ 1.-0.j  0.+0.j -0.-0.j  0.-0.j]\n",
      " [ 0.-0.j  1.+0.j -0.-0.j -0.-0.j]\n",
      " [-0.-0.j -0.+0.j  1.-0.j  0.+0.j]\n",
      " [ 0.+0.j -0.-0.j  0.-0.j  1.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "U = stats.unitary_group.rvs(n)\n",
    "\n",
    "print(np.round(np.dot(U,U.conj().T)),'\\n')\n",
    "\n",
    "U = np.matrix(U)\n",
    "\n",
    "print(np.round(U,3),'\\n')\n",
    "\n",
    "print(np.round(U.H*U))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "183e1082",
   "metadata": {},
   "source": [
    "## Tensor Product"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a770bb2",
   "metadata": {},
   "source": [
    "In the function `scipy.sparse.kron(A,B)` the matrices/vectors `A` and `B` can be sparse or dense. The tensor product is computed as a sparse matrix. The tensor product of two $2\\times2$ matrices $A$ and $B$ is as follows\n",
    "\n",
    "$$ A = \n",
    "\\begin{pmatrix}\n",
    "a_{11} & a_{12} \\\\ \n",
    "a_{21} & a_{22}\n",
    "\\end{pmatrix}, \n",
    "\\ \\ \\ \\ \n",
    "B = \n",
    "\\begin{pmatrix}\n",
    "b_{11} & b_{12} \\\\ \n",
    "b_{21} & b_{22}\n",
    "\\end{pmatrix},$$\n",
    "\n",
    "$$ A \\otimes B = \n",
    "\\begin{pmatrix}\n",
    "a_{11} B & a_{12} B \\\\ \n",
    "a_{21} B & a_{22} B\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "a_{11} \\begin{pmatrix}\n",
    "b_{11} & b_{12} \\\\ \n",
    "b_{21} & b_{22}\n",
    "\\end{pmatrix} & a_{12} \\begin{pmatrix}\n",
    "b_{11} & b_{12} \\\\ \n",
    "b_{21} & b_{22}\n",
    "\\end{pmatrix} \\\\ \n",
    "a_{21} \\begin{pmatrix}\n",
    "b_{11} & b_{12} \\\\ \n",
    "b_{21} & b_{22}\n",
    "\\end{pmatrix} & a_{22} \\begin{pmatrix}\n",
    "b_{11} & b_{12} \\\\ \n",
    "b_{21} & b_{22}\n",
    "\\end{pmatrix}\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "a_{11} b_{11} & a_{11} b_{12} & a_{12} b_{11} & a_{12} b_{12} \\\\ \n",
    "a_{11} b_{21} & a_{11} b_{22} & a_{12} b_{21} & a_{12} b_{22} \\\\ \n",
    "a_{21} b_{11} & a_{21} b_{12} & a_{22} b_{11} & a_{22} b_{12} \\\\ \n",
    "a_{21} b_{21} & a_{21} b_{22} & a_{22} b_{21} & a_{22} b_{22}\n",
    "\\end{pmatrix}.$$\n",
    "The sucture is similar for differently shaped matrices and vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "112427a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result dense = \n",
      " [[ 1  2 -2 -4]\n",
      " [ 3  4 -6 -8]\n",
      " [-1 -2  2  4]\n",
      " [-3 -4  6  8]] \n",
      "\n",
      "Result sparse = \n",
      " [[ 1  2 -2 -4]\n",
      " [ 3  4 -6 -8]\n",
      " [-1 -2  2  4]\n",
      " [-3 -4  6  8]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, -2], [-1, 2]])\n",
    "B = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "result_dense = scipy.sparse.kron(A,B)\n",
    "\n",
    "print('Result dense = \\n',result_dense.todense(),'\\n')\n",
    "\n",
    "A = scipy.sparse.csr_matrix(A)\n",
    "B = scipy.sparse.csr_matrix(B)\n",
    "\n",
    "result_sparse = scipy.sparse.kron(A,B)\n",
    "\n",
    "print('Result sparse = \\n',result_sparse.todense())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1b5b0d0",
   "metadata": {},
   "source": [
    "## Least Square Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca455468",
   "metadata": {},
   "source": [
    "The package `scipy.linalg` also supports least-square fitting, a method that fits data points to a curve. The relavent function is `scipy.linalg.lstsq(A, c)`. To demonstrate how this works, we are going to use this function to approximate the values of $c_1$ and $c_2$ in the following function\n",
    "$$ y = c_1 x^2 + c_2 x.$$\n",
    "We will generate data points by adding noise to this function for certain values of $c_1$ and $c_2$. We will then then try to fit the noisy data with `scipy.linalg.lstsq(A, c)`, comparing our approximate values of $c_i$ to the actual values.\n",
    "\n",
    "Lets set \n",
    "$$ c_1 = 0.5$$ and $$ c_2 = 2.0$$\n",
    "$A$ is a matrix of the set of functions, of $x$, describing $y$. So in this case\n",
    "$$ A = [x^2,x]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "0108a613",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHACAYAAABDKXcJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABEP0lEQVR4nO3de1xUZeLH8c+AMoABpuUgiYqFaV7S1EzUsBTKzM2uW3ax3C1NTck2L7mV9UtIKyojLasty3XtpmYXS8TCyrXISxdrBcuUSiQTGVQcdOb8/jgxiaB5gTnDzPf9es1rX8+ZC99odb495znnsRmGYSAiIiISxEKsDiAiIiJiNRUiERERCXoqRCIiIhL0VIhEREQk6KkQiYiISNBTIRIREZGgp0IkIiIiQU+FSERERIKeCpGIiIgEvQZWB/AFj8fDL7/8QlRUFDabzeo4IiIi4gOGYVBWVkZcXBwhIX8yB2RYbP/+/caUKVOM1q1bG+Hh4UZCQoLxwAMPGG632/saj8dj3H///Ubz5s2N8PBwIzk52fjmm2+O+mcUFhYagB566KGHHnroEYSPwsLCP+0Kls8QTZ8+nWeeeYa5c+fSoUMHvvjiC2655RZiYmIYN24cADNmzCAzM5OXXnqJtm3b8tBDD5GSksLGjRuJior6059R+ZrCwkKio6Pr9J9HRERE/IPT6SQ+Pv6ouoLNMKzd3PXSSy/F4XDwwgsveI9deeWVREZG8sorr2AYBnFxcaSlpTFx4kQAXC4XDoeD6dOnM2LEiD/9GU6nk5iYGEpLS1WIREREgsSxfP9bvqi6T58+5OTkkJ+fD8CXX37JJ598wiWXXALA5s2bKSoqIjU11fseu91OcnIyq1atqvEzXS4XTqezykNERETkcCw/ZTZx4kRKS0tp164doaGhuN1upk2bxnXXXQdAUVERAA6Ho8r7HA4HW7ZsqfEzMzIyeOCBB+o2uIiIiAQMy2eIXn31VebNm8f8+fNZu3Ytc+fO5dFHH2Xu3LlVXnfo1WGGYRz2irHJkydTWlrqfRQWFtZZfhEREan/LJ8huvvuu5k0aRLXXnstAJ06dWLLli1kZGQwbNgwYmNjAXOmqHnz5t73FRcXV5s1qmS327Hb7XUfXkRERAKC5TNEe/furXZvgNDQUDweDwAJCQnExsaSnZ3tfb6iooLc3FySkpJ8mlVEREQCk+UzRIMHD2batGm0bNmSDh06sG7dOjIzMxk+fDhgnipLS0sjPT2dxMREEhMTSU9PJzIykqFDh1qcXkRERAKB5YXoqaee4t5772XUqFEUFxcTFxfHiBEjuO+++7yvmTBhAuXl5YwaNYqSkhJ69uzJsmXLjuq+AiIiIiJ/xvL7EPmC7kMkIiISfOrVfYhERERErKZCJCIiIkFPhUhERER86vHsfGbmFNT43MycAh7PzvdxIhUiERER8bHQEBuZNZSimTkFZGbnExpS842X65LlV5mJiIhIcBnbPxGAzN9ngsb2T/SWofEpbb3P+5IKkYiIiPjcwaUoa8UmKtwey8oQ6JSZiIiIWGRs/0TCQkOocHsICw2xrAyBCpGIiIhYZGZOgbcMVbg9h11o7Qs6ZSYiIiI+d+iaocoxoDVEIiIiEvhqWkBd00JrX1IhEhEREZ9ye4waF1BXjt0e3+8qpr3MREREJCBpLzMRERGRY6BCJCIiIkFPhUhERESCngqRiIiIBD0VIhEREQl6KkQiIiIS9FSIREREJOipEImIiEjQUyESERGRoKdCJCIiItYp3wUet9UpVIhERETEIkXfwLPnQ86DVidRIRIRERELfP0GvJACu7bAt4vBVWZpHO12LyIiIr7jPgA5U2HVU+b49AvhyhfAHmVpLBUiERER8Y09v8Ebt8DmXHPcZzxc+E8ICbU2FypEIiIi4gvbvoQFN0DpVmjYCIbMgg5DrE7lpUIkIiIidevLV+HtsXBgHzRpA3/9NzjOsjpVFSpEIiIiUjfc+2HZvfDZbHOcmApXPAcRjS2NVRMVIhEREal9u3+F12+GLZ+Y4/MnQL/JEOKfF7irEImIiEjt+nkNvHojOH+GsCi4/Blof6nVqY5IhUhERERqz7p/wzt3gtsFTRPh2n/DqWdanepPqRCJiIjIiTtQAR/cA3nPmeMzLzFnhsJjrM11lFSIRERE5MSUbYfXh8HW/5rjfvfA+Xf77XqhmlietHXr1thstmqP0aNHA2AYBlOnTiUuLo6IiAj69evHhg0bLE4tIiIiABTmwZxkswzZo+G6V6HfxHpVhsAPClFeXh7btm3zPrKzswG4+uqrAZgxYwaZmZlkZWWRl5dHbGwsKSkplJVZu+eJiIhI0FvzErx0CZRtg1POhFs/hDMvtjrVcbG8EJ166qnExsZ6H++88w6nn346ycnJGIbBE088wZQpU7jiiivo2LEjc+fOZe/evcyfP9/q6CIiIsHpgAveHmc+3BXQfjDcmgOnnGF1suNmeSE6WEVFBfPmzWP48OHYbDY2b95MUVERqamp3tfY7XaSk5NZtWrVYT/H5XLhdDqrPERERKQWOLfBS4PM2SFs0P8+uOYVyzdnPVF+VYgWL17Mrl27uPnmmwEoKioCwOFwVHmdw+HwPleTjIwMYmJivI/4+Pg6yywiIhI0tq421wv9lGdePXb9G9D3LrDZrE52wvyqEL3wwgsMHDiQuLi4Ksdth/yiDcOoduxgkydPprS01PsoLCysk7wiIiJBwTAg73lzZmj3dmh2Ftz2ESQOsDpZrfGby+63bNnC8uXLWbhwofdYbGwsYM4UNW/e3Hu8uLi42qzRwex2O3a7ve7CioiIBIv9++C9u2DdPHPc4XL4SxbYT7I2Vy3zmxmiF198kWbNmjFo0CDvsYSEBGJjY71XnoG5zig3N5ekpCQrYoqIiASP0p/gxYFmGbKFQMqDcNWLAVeGwE9miDweDy+++CLDhg2jQYM/ItlsNtLS0khPTycxMZHExETS09OJjIxk6NChFiYWEREJcD9+at5scc+vEHEyXPUvOP1Cq1PVGb8oRMuXL2fr1q0MHz682nMTJkygvLycUaNGUVJSQs+ePVm2bBlRUfV7NbuIiIhfMgz47FlYNgU8ByC2E/x1Hpzc2upkdcpmGIZhdYi65nQ6iYmJobS0lOjoaKvjiIiI+Kf95fB2Gny1wBx3uhoGz4SwSEtjHa9j+f73ixkiERERsdiurfDqDbDtS7CFQur/wXmjAuKS+qOhQiQiIhLsfsiF12+G8p0Q2RSufgkSzrc6lU+pEImIiAQrw4D/Pg3Z94LhgeZdzPVCjYPvhsYqRCIiIsGoYi8suQO+ecMcn30dXPo4NIywNpdFVIhERESCTcmPsOAG2P41hDSAizLg3FuDZr1QTVSIREREgsmmHHjzb1BeAo1OhavnQuveVqeynAqRiIhIMDAM+PQJyHnQXC90Wjdzl/qY06xO5hdUiERERAKdaze8NRq+XWyOu94AlzwGDcMtjeVPVIhEREQC2W/fm/cXKv4WQhrCwOnQfXhQrxeqiQqRiIhIoCrINtcL7SuFkxzmKbKWPa1O5ZdUiERERAKNxwOfPAYrpgEGtDgXrnkZoptbncxvqRCJiIgEElcZLBoJ/3vHHHcfDhdPhwZh1ubycypEIiIigWJHASy4HnZshNAwuORR6DbM6lT1ggqRiIhIINi4FBbeBi4nRMXBX1+BFt2tTlVvqBCJiIjUZx4P5E6H3IfNccte5s0WoxzW5qpnVIhERETqq707YeGtsGm5OT73NkidpvVCx0GFSEREpD76ZR28ehOUboUG4XDpE9DlOqtT1VsqRCIiIvXNmrnw3t3gdsHJCeZ6odhOVqeq11SIRERE6ov95fDeP2DdPHPcdiBc/gxENLY0ViBQIRIREakPSn6E126CbV+CLQQumAJ9xkNIiNXJAoIKkYiIiL8ryIY3/w77dkFkU7jyBTj9AqtTBRQVIhEREX/lvaR+OmBA3DnmFhyN461OFnBUiERERPzRoZfUdx8OFz8MDezW5gpQKkQiIiL+5pf18NqNsEuX1PuKCpGIiIg/WfsyvPuP3y+pbw1/nadL6n1AhUhERMQf7N/3+yX1r5jjthf/fkn9ydbmChIqRCIiIlYr2WKeItv2JWCDC6dAn7uOeEn949n5hIbYGNs/sdpzM3MKcHsM7kxpW4ehA4tuXiAiImKlguUwJ9ksQxFN4MaFcP7df3p/odAQG5nZ+czMKahyfGZOAZm/lyU5epohEhERsYLHAytnwEcPczyX1FfODGVm53vHlWVofErbGmeO5PBUiERERHxt705YeBtsyjbH3W6BgdOP+ZL6g0tR1opNVLg9KkPHSafMREREfOmX9eYpsk3Z5iX1Q2bD4CeO+/5CY/snEhYaQoXbQ1hoiMrQcVIhEhER8ZW1r8ALqeb9hU5uDX/Lhi5DT+gjZ+YUeMtQhdtTbU2RHB2dMhMREalr+/fB0rvNewxBrV1Sf+iaocoxoJmiY6RCJCIiUpdKtvy+S/16wGbuUt/3yJfUH42aFlDXtNBajo4KkYiISF0pWA4L/w7lJeYl9Vc+D2f0r5WPdnuMGhdQV47dHqNWfk6w8Is1RD///DM33HADTZs2JTIyki5durBmzRrv84ZhMHXqVOLi4oiIiKBfv35s2LDBwsQiIiJH4PHAR9Ph31eZZSiuK4zIrbUyBHDnEa4mG9s/UTdlPEaWF6KSkhJ69+5Nw4YNWbp0Kd9++y2PPfYYjRs39r5mxowZZGZmkpWVRV5eHrGxsaSkpFBWVmZdcBERkZrs3Qnzr4GP0gHDvKR++AfQuKXVyeQIbIZhWDqnNmnSJD799FM+/vjjGp83DIO4uDjS0tKYOHEiAC6XC4fDwfTp0xkxYkS197hcLlwul3fsdDqJj4+ntLSU6OjouvkHEREROXSX+kGZ0PV6q1MFLafTSUxMzFF9/1s+Q7RkyRK6d+/O1VdfTbNmzejatSvPPfec9/nNmzdTVFREamqq95jdbic5OZlVq1bV+JkZGRnExMR4H/HxR3fXTxERkeN28CX1jVuZl9SrDNUblheiH374gdmzZ5OYmMgHH3zAyJEjGTt2LC+/bF6aWFRUBIDD4ajyPofD4X3uUJMnT6a0tNT7KCwsrNt/CBERCV7798GSO2DJGHC7IPEic71Q885WJ5NjYPlVZh6Ph+7du5Oeng5A165d2bBhA7Nnz+amm27yvs5mq7pJnWEY1Y5Vstvt2O3Hd8dPERGRo1ZHl9SL71n+b6x58+acddZZVY61b9+erVu3AhAbGwtQbTaouLi42qyRiIiIz3h3qV9v3mDxhjcg+c93qRf/ZPm/td69e7Nx48Yqx/Lz82nVqhUACQkJxMbGkp2d7X2+oqKC3NxckpKSfJpVRESk5kvqV8IZA6xOJifA8lNmd955J0lJSaSnp3PNNdfw+eefM2fOHObMmQOYp8rS0tJIT08nMTGRxMRE0tPTiYyMZOjQE9v/RURE5Jjs3QmLRkDBMnPc7Wa4eDo0DLc0lpw4ywtRjx49WLRoEZMnT+bBBx8kISGBJ554guuv/2Nl/oQJEygvL2fUqFGUlJTQs2dPli1bRlRUlIXJRUQkqGz7El69QZfUByjL70PkC8dyHwIREZFq1s2Dd8abV5E1bgV/naeryOqBY/n+t3yGSERExG/t3wdLJ8DaueY48SK44tkT3qVe/I8KkYiISE12bYVXb9Ql9UFChUhERORQm5bDm5W71J/8+y71uooskKkQiYiIVPJ44ONH4cPfN2aN6wrXvKyNWYOACpGIiAj8fkn9SCj4wBzrkvqgokIkIiLy0xp4/WYorbyk/jHoeoPVqcSHVIhERCR4GQZ8Pgc+mAKe/XBya/MUWfOzrU4mPqZCJCIiwWmf09yl/tvF5rjdpTBkFoTHWBpLrKFCJCIiwafoa3htGOz8HkIaQMr/wXm3g81mdTKxiAqRiIgED8OAda/Ae3fDgX0QfRoLWj1I8d6zGVtDGZqZU4DbY3BnSlsLwoov6e5SIiISHCr2wOJR5mmyA/vgjBQY8THFjc8mMzufmTkFVV4+M6eAzOx8QkM0axQMNEMkIiKB79d8eO0m+PU7sIWYd53uMx5CQhjbvykAmdn5AIztn+gtQ+NT2jK2f6KVycVHVIhERCSwff0GLBkL+/dAo2Zw1QuQcH6Vl1SWnszsfLJWbKLC7VEZCjI6ZSYiIoFp/z5zh/o3/2aWodZ9YeQn1cpQpbH9EwkLDaHC7SEsNERlKMioEImISODZuRn+lQpfvGCOz78bbnoLohyHfcvMnAJvGapwe6qtKZLAplNmIiISWL57x1w87SqFiCZwxXOQeOSNWQ9dM1Q5BjRTFCRUiEREJDC498PyqfDfLHPc4ly4+kWIaXHEt9W0gPrgNUUHjyVwqRCJiEj9V/oTvH4L/PS5Oe41BgZMhdCGf/pWt8eocQF15djtMWo7rfghm2EYAf9v2ul0EhMTQ2lpKdHR0VbHERGR2rRpObx5K5TvBHsMDHka2g+2OpX4gWP5/tcMkYiI1E8eN3yUASsfBQyI7QzXzIUmbaxOJvWQCpGIiNQ/ZdvNy+l//Ngcdx8OF2VAw3Brc0m9pUIkIiL1y4+fwBvDYfd2aNgIBj8Jna+2OpXUcypEIiJSP3g88OnjsOIhMDxwanvzFNmpZ1qdTAKACpGIiPi/vTth0QgoWGaOz74OBj0GYY2szSUBQ4VIRET8W2EevH4zOH+CBuFwySPQ9UawaRd6qT0qRCIi4p8MA1bPhux7wXPAvHrsmpchtpPVySQAqRCJiIj/2VcKb42G7942x2ddBn/JgnDdS07qhgqRiIj4l21fwmvDoGQzhDSEi6bBubfpFJnUKRUiERHxD4YBa16CpRPB7YKYeLj6JWjR3epkEgRUiERExHqu3fDuePjqVXOceBFc/gxENrE2lwQNFSIREbFW8f/gtZtgx0awhUL/eyFpHISEWJ1MgogKkYiIWOfLV+GdNNi/F06Khav+Ba17W51KgpAKkYiI+N7+ffD+RHPNEEBCMlz5PJzUzNJYErxUiERExLd++x5eHwZFXwM2SJ4AyRMhJNTqZBLELD9BO3XqVGw2W5VHbGys93nDMJg6dSpxcXFERETQr18/NmzYYGFiERE5bt++BXP6mWUosinc8CZccI/KkFjOL2aIOnTowPLly73j0NA//mDMmDGDzMxMXnrpJdq2bctDDz1ESkoKGzduJCoqyoq4IiJyiMez8wkNsTG2f2K152bmFGAcqGCc5xX4bLZ5MP48c71QzGk+TipSM78oRA0aNKgyK1TJMAyeeOIJpkyZwhVXXAHA3LlzcTgczJ8/nxEjRvg6qoiI1CA0xEZmdj5AlVI0M6eABdmrWHjqHCj7xjyYNBb63wehDa2IKlIjy0+ZARQUFBAXF0dCQgLXXnstP/zwAwCbN2+mqKiI1NRU72vtdjvJycmsWrXqsJ/ncrlwOp1VHiIiUnfG9k9kfEpbMrPzmZlTAJhlaG3Oq+ScdC+xZd9AeAxc+x9I/T+VIfE7ls8Q9ezZk5dffpm2bduyfft2HnroIZKSktiwYQNFRUUAOByOKu9xOBxs2bLlsJ+ZkZHBAw88UKe5RUSkqsqZoczsfGav2MgY26u8FLYEDgBxXc27Tp/c2sqIIodl+QzRwIEDufLKK+nUqRMDBgzg3XffBcxTY5Vsh+xfYxhGtWMHmzx5MqWlpd5HYWFh3YQXEZEqxvZP5LTQUl4MfYjRDZaYB3vcCsM/UBkSv2Z5ITpUo0aN6NSpEwUFBd51RZUzRZWKi4urzRodzG63Ex0dXeUhIiJ1b9Ebr7C4wUTOC/mO3UY4S9ulw6BHoYHd6mgiR+R3hcjlcvHdd9/RvHlzEhISiI2NJTs72/t8RUUFubm5JCUlWZhSRESqcB8g74U0Lvv6Dk61OaFZBxZ2n8ft61t71xSJ+DPL1xD94x//YPDgwbRs2ZLi4mIeeughnE4nw4YNw2azkZaWRnp6OomJiSQmJpKenk5kZCRDhw61OrqIiACU/sQvLwylh/NLsAHdboGLM7ipYQS7IgtqvPpMxN9YXoh++uknrrvuOnbs2MGpp57Keeedx+rVq2nVqhUAEyZMoLy8nFGjRlFSUkLPnj1ZtmyZ7kEkIuIPNi6FxbcTV15CRWgjwi5/Cjpe6X26sgS5PYZVCUWOis0wjID/f6nT6SQmJobS0lKtJxIRqQ0HKmD5VFj9tDlu3gWufhGatLEylUgVx/L9b/kMkYiI1DM7N8Mbt8Av68zxeaNgwFQtnJZ6TYVIRESO3jcL4e1x4HJCeGMYMhvaXWJ1KpETpkIkIiJ/bn85fHAPfPEvcxzfE658ARrHW5tLpJaoEImIyJH9mm+eItv++15kfcabO9Rr+w0JICpEIiJyeOv/A++Oh/17odGpcPmzcEZ/q1OJ1DoVIhERqc61G967G76cb44TzocrnoOoWGtzidQRFSIREamq6BvzFNmOfLCFQL/J0PcuCAm1OplInVEhEhERk2HAmhdh6SRwuyCqublwunVvq5OJ1DkVIhERgX2l5uX0GxaZ4zNS4PJnoNEp1uYS8REVIhGRYPfzGnhjOJT8CCENoP/90GsMhPjd/t8idUaFSEQkWBkGrJ4F2feDZz/EtDS332jR3epkIj6nQiQiEoz27oTFoyB/qTluPxj+8hREnGxtLhGLqBCJiASbravNU2TOnyE0DC5Khx5/B5vN6mQillEhEhEJFh4PfPo4rJgGhhuanG6eImt+ttXJRCynQiQiEgx2F8PC2+CHD81xp2vg0kywR1mbS8RPqBCJiAS6Hz4yy9Du7dAgAgY9Cl2u1ykykYOoEImIBCr3AcidDisfAQw4tb15iqxZe6uTifgdFSIRkUDk/AXe/Dts+dQcn3MTXDwdwiKtzSXip1SIREQCTf4HsGgklO+EsJNg8JPQ6SqrU4n4NRUiEZFAcaACVjwIq54yx7Gd4eqXoOnplsYSqQ9UiEREAkHJj+a9hX5eY47PHQGp/wcN7JbGEqkvVIhEROq7b9+Ct+4AVymEx8BlT5t3nhaRo6ZCJCJSX+3fB8umQN7z5rhFD7jqX9C4pbW5ROohFSIRkfpoxyZ4/WbY/rU57j0OLrwXQhtaGkukvlIhEhGpb756Dd5Og/17ILIpXD4HEgdYnUqkXlMhEhGpLyr2wHsTYP08c9yqD1z5PEQ3tzaXSABQIRIRqQ+2f2ueItuxEbBB8kRIngAhoVYnEwkIKkQiIn7k8ex8QkNsjO2faB4wDFg7F5ZOhAP72N3wFE4a+iIknG9tUJEAo0IkIuJHQkNsZGbnAzC2twPeSYNv3gQg192Z/J6PcmtCTwsTigQmFSIRET9SOTO0bPn73JT3DI33/YSHUKbvv4ZGF9zJ2AFnWpxQJDCFWB1AREQO4vEwNvw93gqfSuN9P/GzcQpXue6l0YV3qQyJ1CHNEImI+IuyInNT1h8+JBT4wHMuEyr+TnloNAsr1xSJSJ3QDJGIiD/I/wBm94YfPoQGEeSccQ8jKsZRHhpNhdvDzJwCqxOKBDTNEImIWGn/Plh+P3z2jDl2dGTeaffxz1UHGJ/SlrH9E5mZU/DHQmvNFInUCb+aIcrIyMBms5GWluY9ZhgGU6dOJS4ujoiICPr168eGDRusCykiUlt+3QjPD/ijDPUcydNnPFulDIFZgsantCUzO18zRSJ1xG8KUV5eHnPmzKFz585Vjs+YMYPMzEyysrLIy8sjNjaWlJQUysrKLEoqInKCDAPWvATPJpt7kUU2haGvwcDpVBBWpQxVqixFbo9hTWaRAGczDMPyP127d+/mnHPOYdasWTz00EN06dKFJ554AsMwiIuLIy0tjYkTJwLgcrlwOBxMnz6dESNGHNXnO51OYmJiKC0tJTo6ui7/UUREjqy8BJaMhe+WmOM2/eDyZyEq1tJYIoHoWL7//WKGaPTo0QwaNIgBA6puTrh582aKiopITU31HrPb7SQnJ7Nq1arDfp7L5cLpdFZ5iIhYbssqmN3HLEMhDSDlQbhhkcqQiB+wfFH1ggULWLt2LXl5edWeKyoqAsDhcFQ57nA42LJly2E/MyMjgwceeKB2g4qIHC/3AcidDh8/CoYHmrSBK1+A086xOpmI/M7SGaLCwkLGjRvHvHnzCA8PP+zrbDZblbFhGNWOHWzy5MmUlpZ6H4WFhbWWWUTkmJRsgZcugZUzzDJ09lAYsVJlSMTPWDpDtGbNGoqLi+nWrZv3mNvtZuXKlWRlZbFx40bAnClq3ry59zXFxcXVZo0OZrfbsdvtdRdcRORofLMQ3k4DVynYo+HSx6HTVVanEpEaWFqI+vfvz9dff13l2C233EK7du2YOHEibdq0ITY2luzsbLp27QpARUUFubm5TJ8+3YrIIiJ/rmIPLJ0A6+aZ4xY94Mrn4eTWlsYSkcM7pkJUWFhIfHx8rf3wqKgoOnbsWOVYo0aNaNq0qfd4Wloa6enpJCYmkpiYSHp6OpGRkQwdOrTWcoiI1Jpf1sObf4PfNgE26HsX9JsEoQ2tTiYiR3BMhahdu3aMHz+eSZMm0ahRo7rKVMWECRMoLy9n1KhRlJSU0LNnT5YtW0ZUVJRPfr6IyFHxeGD1LFg+FTz7ISoOrpgDCX2tTiYiR+GY7kO0atUq7rzzTgoLC5k2bRq33HJLXWarNboPkYjUqd3F5qas3+eY43aXwl+egsgm1uYSCXJ1dh+ipKQkPvvsMx5++GHuu+8+unbtykcffXQiWUVE6reC5TA7ySxDDcJhUCb8dZ7KkEg9c1yX3d90003k5+czePBgBg0axOWXX86mTZtqO5uIiP864IL374F/Xwl7foVmHeC2j6DH3+AItwUREf903PchMgyD1NRUbrvtNpYsWULHjh256667tMeYiAS+X/Ph+f6w+mlzfO5tcOsKaNbe2lwictyOaVH1M888Q15eHnl5eXz33XeEhobSuXNnRo8eTZcuXfj3v//NWWedxaJFi+jevXtdZRYRsYZhwLpXYOlE2L8XIprAkFlw5kCrk4nICTqmRdXx8fGcd9553kf37t2r3QAxPT2d+fPn880339R62OOlRdUicsLKd8Hb4+DbxeY4IdnclDW6+ZHeJSIWOpbv/1rf7X779u3ExcXhdrtr82NPiAqRiJyQravhzb9DaaG5KeuF/4SkcRDiF/tji8hhHMv3f63fqbpZs2asWLGitj9WRMT33AfMDVlzp5v7kJ2cYG7K2qLbn79XROqVWi9ENpuN5OTk2v5YERHf2lUIC2+Frf81x52vhUGPgl03hRUJRJbuZSYi4pc2LIa3x8K+UgiLgkszofM1VqcSkTqkQiQiUqliD7w/Cda+bI5P62ZuytqkjbW5RKTOqRCJiABs+8rclHVHPmCDPnfCBfdoU1aRIKFCJCLBzTBg9WxYfj+4KyCquXk5fRuthRQJJipEIhK8dv8Ki2+HTdnm+MxL4C9Z0KiptblExOdUiEQkOG3KMXeo31MMoXa4aBr0+Lv2IRMJUipEIhJcDlRAzgPw3yxzfGp7uOpf4DjL2lwiYikVIhEJHjs2wZvDYduX5rjH3yH1IWgYYW0uEbGcCpGIBD7DgPX/hvcmwP49EHEyXPY0tBtkdTIR8RMqRCIS2Mp3wTt3woaF5rh1X7hiDkTHWRpLRPyLCpGIBK4tq2DhCCjdCrZQuHAK9E6DkFCrk4mIn1EhEpHAc6ACch+GTx43N2Vt3MpcON2iu9XJRMRPqRCJSGD5Nd/clHXbenPc5Qa4OAPCoy2NJSL+TYVIRAKDYcAXL8AH/4QD5RDeGAY/CR2GWJ1MROoBFSIRqf92F8NbY6DgA3Pcph8Mma2F0yJy1FSIRKR+27jULEN7d5h3nB4wFXqOhJAQq5OJSD2iQiQi9VPFHvhgCqx50Rw36wBXPgeODjyenU9oiI2x/ROrvW1mTgFuj8GdKW19HFhE/Jn+E0pE6p+f18Kz5/9RhnqNgVtXgKMDAKEhNjKz85mZU1DlbTNzCsj8vSyJiBxMM0QiUn943PBJJnz0MHgOQFQcXD7bXDN0kMqZoczsfO+4sgyNT2lb48yRiAQ3FSIRqR9KfjRvsli42hyfNQQufRwim9T48oNLUdaKTVS4PSpDInJYOmUmIv7NMGD9f2B2H7MMhUXBkGfg6pcOW4Yqje2fSFhoCBVuD2GhISpDInJYKkQi4r/27oTXb4bFI6GiDOLPg9s/gS7Xge3P1wHNzCnwlqEKt6famiIRkUo6ZSYi/umHj2DR7VD2C4Q0gH6Toc+dR70P2aFrhirHgGaKRKQaFSIR8S/798GK/4P/ZpnjpmeYu9Of1u2oP6KmBdQ1LbQWEamkQiQi/mP7BnjzVijeYI67D4fUhyCs0TF9jNtj1LiAunLs9hi1EldEAofNMIyA/5vB6XQSExNDaWkp0dHa4FHE73g88NkzsHwquF0QeQpclgVnDrQ6mYjUY8fy/a8ZIhGxlvMXWHy7uWYIIPEiswyd1MzSWCISXCy/ymz27Nl07tyZ6OhooqOj6dWrF0uXLvU+bxgGU6dOJS4ujoiICPr168eGDRssTCwitebbt2B2klmGGkTAoEwY+qrKkIj4nOWFqEWLFjz88MN88cUXfPHFF1x44YVcdtll3tIzY8YMMjMzycrKIi8vj9jYWFJSUigrK7M4uYgct31OWDwKXrsJykug+dkwYiX0+NtRXU4vIlLb/HINUZMmTXjkkUcYPnw4cXFxpKWlMXHiRABcLhcOh4Pp06czYsSIGt/vcrlwuVzesdPpJD4+XmuIRPzB1s9g4a2wawtgg77jIXkSNAizOpmIBJhjWUNk+QzRwdxuNwsWLGDPnj306tWLzZs3U1RURGpqqvc1drud5ORkVq1addjPycjIICYmxvuIj4/3RXwRORL3flgxDV682CxDMS3hlveg/30qQyJiOb8oRF9//TUnnXQSdrudkSNHsmjRIs466yyKiooAcDgcVV7vcDi8z9Vk8uTJlJaWeh+FhYV1ml9E/sRv38O/LoKVM8DwQOdrzTtOt0qyOpmICOAnV5mdeeaZrF+/nl27dvHmm28ybNgwcnNzvc/bDllTYBhGtWMHs9vt2O32OssrIkfJMGDtXHh/MuzfC+Ex5oasHa+0OpmISBV+UYjCwsI444wzAOjevTt5eXk8+eST3nVDRUVFNG/e3Pv64uLiarNGIuJn9uyAJXfAxvfMceu+cPkzENPC2lwiIjXwi1NmhzIMA5fLRUJCArGxsWRnZ3ufq6ioIDc3l6QkTbWL+K38ZTCrl1mGQsPMu03ftERlSET8luUzRPfccw8DBw4kPj6esrIyFixYwEcffcT777+PzWYjLS2N9PR0EhMTSUxMJD09ncjISIYOHWp1dBE5VMVeyL4P8p4zx6e2hyufg9hO1uYSEfkTlhei7du3c+ONN7Jt2zZiYmLo3Lkz77//PikpKQBMmDCB8vJyRo0aRUlJCT179mTZsmVERUVZnFxEqvhlvXk5/Q5z81R63g4D7oeGEZbGEhE5Gn55H6Lapr3MROqQxw2rZpqX1Hv2w0mxMGQWnNHf6mQiEuS0l5mI+MaurbBoJGz51By3uxQGz4RGTa3NJSJyjFSIROT4fPUavHsXuJwQdhIMnA5drtfWGyJSL6kQicixKd9lFqFv3jDHLXrAFXOgSRtLY4mInAgVIhE5eps/Nk+ROX8CWygkT4S+d0Go/ioRkfrNL+9DJCLWezw7n5k5BebggAuW3QtzB4PzJ3aFt+A/nZ6DfhNVhkQkIOhvMhGpUWiIjczsfJruKeD6n6fB9m8A+CZ2CNf8+BdGxnS2OKGISO1RIRKRGo29oA3dCl+ix5rZYHNDZFPeaTWJMetOY3xKW8b2T7Q6oohIrVEhEpHqfvseFt9O78LPwAbL3N24v/RWtq2LVhkSkYCkQiQifzAMyHve3H5j/14Ii4KB0xnzxslUuA3CQkNUhkQkIGlRtYiYSn+GeVfAe/8wy1DrvjBqFTN3nustQxVuzx8LrUVEAohmiESCnWHA16+bRWhfKTQIhwFT4dwRzPzwezKz872nyWbmFJCZbe5VppkiEQkkKkQiwWzPb/DunfDtW+Y47hy4/Fk4ta23/By8Zqjyf1WKRCTQqBCJBKuNS2HJWNhTDCENzJss9hnvva+Q22PUuIC6cuz2BPy+0CISRLTbvUiw2eeED+6Bda+Y41PbmbNCcV0sjSUiUtu0272I1Gzzx7B4FJRuBWzQazRceC80DLc6mYiIpVSIRILB/nLI+T9Y/bQ5btwShjwDrXtbm0tExE+oEIkEup/Xmhuy7thojs8ZBhdNA3uUtblERPyICpFIoHLvh5WPwspHwHDDSQ74y1PQ9iKrk4mI+B0VIpFA9OtGWHgbbFtvjjtcDoMyIbKJpbFERPyVCpFIIPF44LPZsPwBcLsgvDEMegw6XWV1MhERv6ZCJBIoSrbAW6Phx4/N8RkD4C9ZEN3c2lwiIvWACpFIfWcYsG4evD8ZKsqgYSO46CHodgvYbFanExGpF1SIROqzsu3w9jjIX2qO48+Dy2dDkzbW5hIRqWdUiETq2OPZ+YSG2Grc92tmTgFuj8GdKW2P/YM3LIZ37oTynRAaBhdMgaQ7ICT0xEOLiASZEKsDiAS60BAbmdn5zMwpqHK8cvPU0JBjPK1VXgJv3gqvDzPLUGwnuO0j6JOmMiQicpw0QyRSx2raIb6mneSPyqYceGsMlP0CthBzM9bkidAgrC6ii4gEDRUiER84uBRlrdhEhdtzbGWoYg9k3wd5z5vjJqebG7LG96ijxCIiwUWnzER8ZGz/RMJCQ6hwewgLDTn6MrT1M3imzx9l6NzbYOQnKkMiIrVIhUjER2bmFHjLUIXbU21NUTUHXOYNFl+8GHb+ANGnwY2L4ZJHICzSJ5lFRIKFTpmJ+MCha4Yqx0DNM0VF38CiEbD9G3N89nVw8cMQ0dh3oUVEgogKkUgdq2kBdU0LrQHwuOHTJ+HDdPDsh8imMPhJaD/YkuwiIsFChUikjrk9Ro0LqCvHbo9hHvjte1h8OxR+Zo7PHASDn4CTmvkwrYhIcLIZhmFYHaKuOZ1OYmJiKC0tJTo62uo4IlUZhrlgOvs+2L8X7NEwcLp5mkxbb4iIHLdj+f7XDJGIlUp/hiVj4PsV5jjhfLhsFjSOtzaXiEiQsfwqs4yMDHr06EFUVBTNmjVjyJAhbNy4scprDMNg6tSpxMXFERERQb9+/diwYYNFiUVqgWHAV6/B7F5mGWoQDgNnwI1vqQyJiFjA8kKUm5vL6NGjWb16NdnZ2Rw4cIDU1FT27Nnjfc2MGTPIzMwkKyuLvLw8YmNjSUlJoayszMLkIsdpz2/mthsLb4V9pRB3Doz4GHqOgBDL/0iKiAQlv1tD9Ouvv9KsWTNyc3M5//zzMQyDuLg40tLSmDhxIgAulwuHw8H06dMZMWLEn36m1hCJ39i4FJaMhT3FENIAkidBnzshVGevRURq27F8//vdf46WlpYC0KRJEwA2b95MUVERqamp3tfY7XaSk5NZtWpVjZ/hcrlwOp1VHiKWKi+BRSPhP9eaZejU9vD3HEi+W2VIRMQP+FUhMgyD8ePH06dPHzp27AhAUVERAA6Ho8prHQ6H97lDZWRkEBMT433Ex2tNhlgo/wOY1Qu+/A9gg6Q7zN3p47pYHExERCr51X+ajhkzhq+++opPPvmk2nO2Qy4/Ngyj2rFKkydPZvz48d6x0+lUKRLfKy+B9++BL+eb46ZnwJDZEH+utblERKQavylEd9xxB0uWLGHlypW0aNHCezw2NhYwZ4qaN2/uPV5cXFxt1qiS3W7HbrfXbWCRI8n/AN4eB2XbMGeFxsAFU6BhhNXJRESkBpafMjMMgzFjxrBw4UJWrFhBQkJClecTEhKIjY0lOzvbe6yiooLc3FySkpJ8HVfkyMp3weJRMP8asww1PQP+tgxSH1IZEhHxY5bPEI0ePZr58+fz1ltvERUV5V0XFBMTQ0REBDabjbS0NNLT00lMTCQxMZH09HQiIyMZOnSoxelFDpK/DN4e+8esUK/RcOE/VYREROoBywvR7NmzAejXr1+V4y+++CI333wzABMmTKC8vJxRo0ZRUlJCz549WbZsGVFRUT5OK1KD8l3wwT2w/t/muOkZ5t2mW/a0NJaIiBw9v7sPUV3QfYikzuQv+32t0C9oVkhExL9oLzORuqZZIRGRgKJCJHKsCrLNu01rVkhEJGCoEIkcrfJd8MEUWD/PHDc5HYbMgpbnWRpLREROnAqRyNGoaVbogikQFml1MhERqQUqRCJHUr4Llk2BdZoVEhEJZCpEIodz6KzQeaPMtUKaFRIRCTgqRCKH2ldqXkF28KzQZU9Dq17W5hIRkTqjQiRysILl5t2mnT+jWSERkeChQiQCNcwKtTHvK6RZIRGRoKBCJLJpublWyDsrdDtceK9mhUREgogKkQSvfaXmfYXWvWKOm7T5fa1QkrW5RETE51SIJDhpVkhERA6iQiTBRbNCIiJSAxUiqRcez84nNMTG2P6J1Z6bmVOA22NwZ0rbI3/IobNCPUdC//s0KyQiIoRYHUDkaISG2MjMzmdmTkGV4zNzCsj8vSwd1j4nLLkD5l1plqGTE+CW92DgwypDIiICaIZI6onKmaHM7HzvuLIMjU9pW+PMEQCbcn6fFfoJzQqJiMjhqBBJvXFwKcpasYkKt+fwZWif09yDbO3L5vjkBHOtUOvePkwsIiL1hU6ZSb0ytn8iYaEhVLg9hIWG1FyGNuXArF5/lKGeI+H2T1WGRETksFSIpF6ZmVPgLUMVbk/VNUX7nObpsXlXmKfITm4NN78HA6dDWCPLMouIiP/TKTOpNw5dM1Q5BhjbautBa4U4aK2QipCIiPw5FSKpF2paQD22fyJhB3YT/eHd0OBD84Untzb3INPpMREROQYqRFIvuD1G9QXU369g5LdjoUGhOT53BAy4X7NCIiJyzFSIjkOt3CRQjkmV3+c+J2TfC2teMscnt/79CrI+VkQTEZEAoEXVx+GEbhIoJyZ/Gcw6748ydO4IuH2VypCIiJwQzRAdh+O+SaAcvz2/wfuT4OvXzLFmhUREpBapEB2nY7pJoBw/w4ANC+G9CbB3B9hC4LxRcMEU3W1aRERqjU6ZnYCjukmgHD/nL7BgKLwx3CxDzc6Cvy2Hi6apDImISK1SIToBR7xJoBw/wzDXCD3dEza+ByENod89cFsutOhmdToREQlAOmV2nI54k0DNFB2/nT+YN1j88WNzfFp3uCwLmrW3NpeIiAQ0FaLjcLibBAJ+X4r89pYBHjesngUrpsGBcmgYCRfeCz1HQEio7/OIiEhQUSE6DjXeJJA/SpDbY1gR66hU3jIAqpa2g0uez23fAG+NgV/WmuOEZBj8JDRJ8H0WEREJSipEx+FIMyj+OjNUya9uGXDABR8/Zj48B8AeYy6Y7noD2HQvJxER8R0VoiDkF7cM+OkLc1bo1+/McbtL4ZJHIbq57zKIiIj8TleZBSnLbhlQsQfevweeH2CWoUanwtUvwV/nqQyJiIhlLC9EK1euZPDgwcTFxWGz2Vi8eHGV5w3DYOrUqcTFxREREUG/fv3YsGGDNWEDiCW3DPjhI5jVC1Y/DRhw9nUw+nPocLlOkYmIiKUsL0R79uzh7LPPJisrq8bnZ8yYQWZmJllZWeTl5REbG0tKSgplZWU+Tho4Dl4zlD9tIONT2ta4N1utKd9lnh57+TLYtQVi4uH6N+HyZyCySd38TBERkWNg+RqigQMHMnDgwBqfMwyDJ554gilTpnDFFVcAMHfuXBwOB/Pnz2fEiBG+jBoQfH7LgO/egXfvgt1F5vjc26D/fWCPqr2fISIicoIsL0RHsnnzZoqKikhNTfUes9vtJCcns2rVqsMWIpfLhcvl8o6dTmedZ60vfHbLgN3F8N7d8O1ic9w0Ef7yFLTqVTufLyIiUov8uhAVFZmzCg6Ho8pxh8PBli1bDvu+jIwMHnjggTrNVl/V+S0DDAO+XGDuTL9vF9hCofc4SJ4IDcNP/PNFRETqgOVriI6G7ZAFt4ZhVDt2sMmTJ1NaWup9FBYW1nVEAdi1Ff59FSweaZah2M5w24cw4H6VIRER8Wt+PUMUGxsLmDNFzZv/cUl2cXFxtVmjg9ntdux2e53nk995PPDFC7B8KlTshlA79JsESXdAaEOr04mIiPwpv54hSkhIIDY2luzsbO+xiooKcnNzSUpKsjCZeO0ogJcugff+YZahlr3g9k+h73iVIRERqTcsnyHavXs3mzZt8o43b97M+vXradKkCS1btiQtLY309HQSExNJTEwkPT2dyMhIhg4damFqwb0fVs2Ej6aD2wVhJ8GAqdD9bxDi1z1bRESkGssL0RdffMEFF1zgHY8fPx6AYcOG8dJLLzFhwgTKy8sZNWoUJSUl9OzZk2XLlhEVpcu2LbPtS3hrNBR9bY7PGACXPgGN4y2NJSIicrxshmH479bstcTpdBITE0NpaSnR0dFWx6m/9pdD7nT4dCYYbog4GS6eDp2v0Z2mRUTE7xzL97/lM0RST2z5LywZA7/9fnqzw+Uw8BE46VRrc4mIiNQCFSI5MlcZLH8A8p4zxyfFwqDHoP2l1uYSERGpRSpEcngF2fB2Gjh/MsfnDIOUByGisZWpREREap0KkVS3dye8Pxm+WmCOT24Ng2dCm2RLY4mIiNQVFSL5g2HAhkXmHmR7d4AtBM4bBRfcA2GNrE4nIiJSZ1SIxOTcZu5Kv/Fdc3xqe7gsC1p0tzaXiIiID6gQBTvDgLUvw7J7wVUKIQ3h/H9An/HQIMzqdCIiIj6hQhTMdm6Gt8fC5pXm+LRu8JcscJxlbS4REREfUyEKRu79sHoWfJgBB8qhQQRc+E8473YICbU6nYiIiM+pEAWbn9fAknGw/fdtN1r3hb/MhCZtrM0lIiJiIRWiYOEqgxUPwWfPAoa57UbqNOgyVNtuiIhI0FMhCgb/ew/e+wc4fzbHna+Fi6ZBo1OszSUiIuInVIgCmfMXWDoBvnvbHJ/cGi59HE6/0NJYIiIi/kaFKBB53PDFv8w9yCrKIKQBJN0B50+AsEir04mIiPgdFaJAs30DvD0Ofsozx6d1h8FPQmxHa3OJiIj4MRWiQLG/HHJnwKqZ4DkAYVEw4H7oPlyX0ouIiPwJFaJA8P2H8M6dULLZHLe7FC55BKLjrM0lIiJST6gQ1Wd7dsAHU/7YlT4qDgY9Cu0GWZtLRESknlEhqo8MA778j1mGyncCNug5wrzbtD3K6nQiIiL1jgpRffPb9/BO2h/7jzk6mYumW3SzNJaIiEh9pkJUXxyogE+fhJWPgNtl7j/WbxL0Gg2hDa1OJyIiUq+pENUHW1ebl9L/+j9zfPqFMCgTmiRYm0tERCRAqBD5s/JdkPOAeZNFgMhT4OKHodNV2n9MRESkFqkQ+SPDgG8Xw9KJsHu7eazrDZDyfxDZxNJoIiIigUiFyN/sKjQ3Ys1/3xw3PQMufQIS+loaS0REJJCpEPkL9wH4/FlYMQ3274GQhtB3PPQZDw3DrU4nIiIS0FSI/MEv6+HtsbDtS3PcMgkGPwGnnmllKhERkaChQmQl1274KANWzwLDA+ExkPIgdL0JQkKsTiciIhI0VIiskr8M3r0LSrea4w5XmFeQRTmszSUiIhKEVIh8rWw7vD8RNiwyxzEt4dJMSEyxNpeIiEgQUyHyFY8H1s6F5ffDvlKwhcB5o+CCeyCskdXpREREgpoKkS/8utG80/TW/5rj5l3M/cfiuliZSkRERH6nQlSX9u+Djx+DTx4Hz35o2Mjckf7c2yBUv3oRERF/oW/lurL5Y3NX+t82meO2F8Mlj0LjeEtjiYiISHX15truWbNmkZCQQHh4ON26dePjjz+2OlLN9u6ExaNh7qVmGTopFq6eC9ctUBkSERHxU/WiEL366qukpaUxZcoU1q1bR9++fRk4cCBbt261OtofDAO+eg2yesD6eeax7n+DMZ9DhyHajFVERMSP2QzDMKwO8Wd69uzJOeecw+zZs73H2rdvz5AhQ8jIyPjT9zudTmJiYigtLSU6Orr2A+7cDO+Oh+9XmONT25uLplv2rP2fJSIiIkflWL7//X4NUUVFBWvWrGHSpElVjqemprJq1aoa3+NyuXC5XN6x0+msm3DuA/Dfp+Cj6XCgHELtkHw3JI2DBmF18zNFRESk1vn9KbMdO3bgdrtxOKrewdnhcFBUVFTjezIyMoiJifE+4uPraO2OzQb/e9csQ637wqj/wvl3qwyJiIjUM35fiCrZDlmDYxhGtWOVJk+eTGlpqfdRWFhYN6FCQs1TY0Nmw7C3oenpdfNzREREpE75/SmzU045hdDQ0GqzQcXFxdVmjSrZ7Xbsdrsv4oGjg/kQERGResvvZ4jCwsLo1q0b2dnZVY5nZ2eTlJRkUSoREREJJH4/QwQwfvx4brzxRrp3706vXr2YM2cOW7duZeTIkVZHExERkQBQLwrRX//6V3777TcefPBBtm3bRseOHXnvvfdo1aqV1dFEREQkANSL+xCdqDq/D5GIiIj4nWP5/vf7NUQiIiIidU2FSERERIKeCpGIiIgEPRUiERERCXoqRCIiIhL0VIhEREQk6KkQiYiISNBTIRIREZGgp0IkIiIiQU+FSERERIJevdjL7ERV7k7idDotTiIiIiK+Uvm9fzS7lAVFISorKwMgPj7e4iQiIiLia2VlZcTExBzxNUGxuavH4+GXX34hKioKm81Wq5/tdDqJj4+nsLBQG8f6iH7nvqffuTX0e/c9/c59ry5/54ZhUFZWRlxcHCEhR14lFBQzRCEhIbRo0aJOf0Z0dLT+8PiYfue+p9+5NfR79z39zn2vrn7nfzYzVEmLqkVERCToqRCJiIhI0FMhOkF2u537778fu91udZSgod+57+l3bg393n1Pv3Pf85ffeVAsqhYRERE5Es0QiYiISNBTIRIREZGgp0IkIiIiQU+FSERERIKeCtEJmjVrFgkJCYSHh9OtWzc+/vhjqyMFrIyMDHr06EFUVBTNmjVjyJAhbNy40epYQSUjIwObzUZaWprVUQLazz//zA033EDTpk2JjIykS5curFmzxupYAevAgQP885//JCEhgYiICNq0acODDz6Ix+OxOlrAWLlyJYMHDyYuLg6bzcbixYurPG8YBlOnTiUuLo6IiAj69evHhg0bfJpRhegEvPrqq6SlpTFlyhTWrVtH3759GThwIFu3brU6WkDKzc1l9OjRrF69muzsbA4cOEBqaip79uyxOlpQyMvLY86cOXTu3NnqKAGtpKSE3r1707BhQ5YuXcq3337LY489RuPGja2OFrCmT5/OM888Q1ZWFt999x0zZszgkUce4amnnrI6WsDYs2cPZ599NllZWTU+P2PGDDIzM8nKyiIvL4/Y2FhSUlK8e5H6hCHH7dxzzzVGjhxZ5Vi7du2MSZMmWZQouBQXFxuAkZuba3WUgFdWVmYkJiYa2dnZRnJysjFu3DirIwWsiRMnGn369LE6RlAZNGiQMXz48CrHrrjiCuOGG26wKFFgA4xFixZ5xx6Px4iNjTUefvhh77F9+/YZMTExxjPPPOOzXJohOk4VFRWsWbOG1NTUKsdTU1NZtWqVRamCS2lpKQBNmjSxOEngGz16NIMGDWLAgAFWRwl4S5YsoXv37lx99dU0a9aMrl278txzz1kdK6D16dOHnJwc8vPzAfjyyy/55JNPuOSSSyxOFhw2b95MUVFRle9Tu91OcnKyT79Pg2Jz17qwY8cO3G43DoejynGHw0FRUZFFqYKHYRiMHz+ePn360LFjR6vjBLQFCxawdu1a8vLyrI4SFH744Qdmz57N+PHjueeee/j8888ZO3Ysdrudm266yep4AWnixImUlpbSrl07QkNDcbvdTJs2jeuuu87qaEGh8juzpu/TLVu2+CyHCtEJstlsVcaGYVQ7JrVvzJgxfPXVV3zyySdWRwlohYWFjBs3jmXLlhEeHm51nKDg8Xjo3r076enpAHTt2pUNGzYwe/ZsFaI68uqrrzJv3jzmz59Phw4dWL9+PWlpacTFxTFs2DCr4wUNq79PVYiO0ymnnEJoaGi12aDi4uJqLVdq1x133MGSJUtYuXIlLVq0sDpOQFuzZg3FxcV069bNe8ztdrNy5UqysrJwuVyEhoZamDDwNG/enLPOOqvKsfbt2/Pmm29alCjw3X333UyaNIlrr70WgE6dOrFlyxYyMjJUiHwgNjYWMGeKmjdv7j3u6+9TrSE6TmFhYXTr1o3s7Owqx7Ozs0lKSrIoVWAzDIMxY8awcOFCVqxYQUJCgtWRAl7//v35+uuvWb9+vffRvXt3rr/+etavX68yVAd69+5d7XYS+fn5tGrVyqJEgW/v3r2EhFT9OgwNDdVl9z6SkJBAbGxsle/TiooKcnNzffp9qhmiEzB+/HhuvPFGunfvTq9evZgzZw5bt25l5MiRVkcLSKNHj2b+/Pm89dZbREVFeWfnYmJiiIiIsDhdYIqKiqq2RqtRo0Y0bdpUa7fqyJ133klSUhLp6elcc801fP7558yZM4c5c+ZYHS1gDR48mGnTptGyZUs6dOjAunXryMzMZPjw4VZHCxi7d+9m06ZN3vHmzZtZv349TZo0oWXLlqSlpZGenk5iYiKJiYmkp6cTGRnJ0KFDfRfSZ9ezBainn37aaNWqlREWFmacc845ugS8DgE1Pl588UWrowUVXXZf995++22jY8eOht1uN9q1a2fMmTPH6kgBzel0GuPGjTNatmxphIeHG23atDGmTJliuFwuq6MFjA8//LDGv7+HDRtmGIZ56f39999vxMbGGna73Tj//PONr7/+2qcZbYZhGL6rXyIiIiL+R2uIREREJOipEImIiEjQUyESERGRoKdCJCIiIkFPhUhERESCngqRiIiIBD0VIhEREQl6KkQiIiIS9FSIREREJOipEImIiEjQUyESERGRoKdCJCIB6T//+Q/h4eH8/PPP3mN///vf6dy5M6WlpRYmExF/pM1dRSQgGYZBly5d6Nu3L1lZWTzwwAM8//zzrF69mtNOO83qeCLiZxpYHUBEpC7YbDamTZvGVVddRVxcHE8++SQff/yxypCI1EgzRCIS0M455xw2bNjAsmXLSE5OtjqOiPgprSESkYD1wQcf8L///Q+3243D4bA6joj4Mc0QiUhAWrt2Lf369ePpp59mwYIFREZG8vrrr1sdS0T8lNYQiUjA+fHHHxk0aBCTJk3ixhtv5KyzzqJHjx6sWbOGbt26WR1PRPyQZohEJKDs3LmT3r17c/755/Pss896j1922WW4XC7ef/99C9OJiL9SIRIREZGgp0XVIiIiEvRUiERERCToqRCJiIhI0FMhEhERkaCnQiQiIiJBT4VIREREgp4KkYiIiAQ9FSIREREJeipEIiIiEvRUiERERCToqRCJiIhI0Pt/Au8ilNGcCOUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c1, c2 = 0.5, 2.0\n",
    "\n",
    "x = np.arange(11)\n",
    "\n",
    "y = c1*x**2 + c2*x\n",
    "\n",
    "noise = y + 0.2 * np.max(y) * rnd.random(len(y))\n",
    "\n",
    "A = np.c_[x**2,x]\n",
    "\n",
    "c, resid, rank, sigma = linalg.lstsq(A, noise)\n",
    "\n",
    "y2 = c[0]*x**2 + c[1]*x\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_axes([.1,.1,.8,.8])\n",
    "\n",
    "ax.plot(x,noise,'x',x,y2)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac97f3a4",
   "metadata": {},
   "source": [
    "## Singular Value Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d60a92b2",
   "metadata": {},
   "source": [
    "Singular value decomposition is an extension of the eigenvalue, eigenvector problem to non-square matrices. The following matrices, and vector, are generated using `scipy.linalg.svd()`\n",
    "$$ A = USv, $$\n",
    "where $A$ is a non-square matrix, $U$ is a square matrix, $S$ is a diagonal matrix, and $v$ is a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "edbc7020",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [[0.87786746 0.51051299]\n",
      " [0.03760617 0.49497429]\n",
      " [0.1804092  0.32190255]] \n",
      "\n",
      "U =  [[-0.89903833  0.40494623 -0.1665792 ]\n",
      " [-0.30997582 -0.85728672 -0.41106504]\n",
      " [-0.30926537 -0.3179277   0.89625717]] \n",
      "\n",
      "Diagonal Elements =  [1.11390874 0.41601034] \n",
      "\n",
      "v =  [[-0.76908259 -0.63914942]\n",
      " [ 0.63914942 -0.76908259]] \n",
      "\n",
      "S =  [[1.11390874 0.        ]\n",
      " [0.         0.41601034]\n",
      " [0.         0.        ]] \n",
      "\n",
      "A =  [[0.87786746 0.51051299]\n",
      " [0.03760617 0.49497429]\n",
      " [0.1804092  0.32190255]]\n"
     ]
    }
   ],
   "source": [
    "A = rnd.random((3,2))\n",
    "print('A = ',A,'\\n')\n",
    "\n",
    "U, D, v = linalg.svd(A)\n",
    "\n",
    "print('U = ',U,'\\n')\n",
    "print('Diagonal Elements = ',D,'\\n')\n",
    "print('v = ',v,'\\n')\n",
    "\n",
    "S = scipy.sparse.lil_matrix((3,2))\n",
    "S[0,0] = D[0]\n",
    "S[1,1] = D[1]\n",
    "\n",
    "S = S.todense()\n",
    "print('S = ', S, '\\n')\n",
    "\n",
    "print('A = ',U*S*v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bad9d3ca",
   "metadata": {},
   "source": [
    "## Schur Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b01f49",
   "metadata": {},
   "source": [
    "The Schur Decomposition of $n\\times n$ matrices gives said matices in the form\n",
    "$$A = Z T Z^\\dagger $$\n",
    "\n",
    "Where $Z$ is a unitary matrix, and $T$ is an upper-triangular(or quasi-upper) matrix. We can request these matrices in a real or complex form using `scipy.linalg.schur(A)` and `scipy.linalg.schur(A,'complex')` or by using `scipy.linalg.rsf2csf(A)` for the complex form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "0b4138e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =  [[0.63015759 0.07866373 0.94700434]\n",
      " [0.8373407  0.96551313 0.74703415]\n",
      " [0.6212803  0.14236477 0.48691623]] \n",
      "\n",
      "A =  [[0.63015759 0.07866373 0.94700434]\n",
      " [0.8373407  0.96551313 0.74703415]\n",
      " [0.6212803  0.14236477 0.48691623]] \n",
      "\n",
      "A =  [[0.63015759+0.j 0.07866373+0.j 0.94700434+0.j]\n",
      " [0.8373407 +0.j 0.96551313+0.j 0.74703415+0.j]\n",
      " [0.6212803 +0.j 0.14236477+0.j 0.48691623+0.j]]\n"
     ]
    }
   ],
   "source": [
    "A = rnd.random((3,3))\n",
    "print('A = ', A, '\\n')\n",
    "\n",
    "T, Z = linalg.schur(A)\n",
    "Tc, Zc = linalg.schur(A, 'complex')\n",
    "\n",
    "T, Z, Tc, Zc = map(np.matrix,(T,Z, Tc, Zc))\n",
    "\n",
    "print('A = ',Z*T*Z.H, '\\n')\n",
    "\n",
    "print('A = ',Zc*Tc*Zc.H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "d1e9d7c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T =  [[ 1.60854984  0.05269232  0.97821817]\n",
      " [ 0.         -0.20330126  0.29601471]\n",
      " [ 0.          0.          0.67733836]]\n"
     ]
    }
   ],
   "source": [
    "print('T = ',T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "ac0e94b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z =  [[ 0.3801452   0.76145013  0.52505554]\n",
      " [ 0.86755609 -0.09672713 -0.48784248]\n",
      " [ 0.3206806  -0.64096612  0.69737105]]\n"
     ]
    }
   ],
   "source": [
    "print('Z = ',Z)"
   ]
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}