{
"cells": [
{
"cell_type": "markdown",
"id": "8d8e8475",
"metadata": {},
"source": [
"# A Computational Exploration of the Numerical Range\n",
"\n",
"(C) 2021 Andreas Kloeckner\n",
"\n",
"\n",
"MIT License\n",
"Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the \"Software\"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.\n",
" \n",
"THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\n",
"\n",
"\n",
"The *numerical range* (also called the *field of values*) of a matrix is the set of values that the Rayleigh \"quotient\" can assume for vectors of 2-norm 1:\n",
"\n",
"$$\\let\\b=\\boldsymbol\\{\\b x^H A \\b x: \\|\\b x\\|_2=1\\}$$\n",
"\n",
"$\\b x^H$ denotes the complex conjugate transpose of $\\b x$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "9a605f72",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import numpy.linalg as la\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "b129b828",
"metadata": {},
"source": [
"Generate a random $3\\times 3$ matrix and a big bunch of random, complex vectors with 2-norm 1:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "a739a86e",
"metadata": {},
"outputs": [],
"source": [
"n = 3\n",
"np.random.seed(18)\n",
"A = np.random.randn(n, n)\n",
"X = np.random.randn(n, 200000) + 1j* np.random.randn(n, 200000)\n",
"X = X/la.norm(X, 2, axis=0)"
]
},
{
"cell_type": "markdown",
"id": "12ee7b68",
"metadata": {},
"source": [
"In the following plot, the eigenvalues are orange dots, the norm is represented by the black circle. As you can see, it's quite easy to manufacture vectors that produce Rayleigh \"quotients\" bigger than the spectral radius. (but not bigger than the norm!)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "c8f72f8a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-2.521459576466905, 2.521459576466905)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD4CAYAAAAn+OBPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoLUlEQVR4nO3deVzU1f4/8NeZgXFXEFEB90ICw1IQEU0TXBBlMFPQrKup19S8erMyTb+V10raVH5ecU2xxXJNByFMUEkMFCa97qLpTRFyAVRkH+b8/tDhorHZfD5zPst5Ph48wgHO523ymnM+53M+50MopeA4Tj00rAvgOM62eOg5TmV46DlOZXjoOU5leOg5TmXsWBy0VatWtFOnTiwOzXGqYDQab1FKnav7GpPQd+rUCRkZGSwOzXGqQAj5vaav8eE9x6kMDz3HqQwPPcepDA89x6kMDz3HqQwPPcepDA89x6kMk+v0nPRQSlFQUICcnBwUFRXBZDLBZDKhvLwcJpMJWq0WdnZ2sLOzg729PXQ6HZydneHs7AyNhvcdcsJDrwKUUmRlZSEzMxM5OTnIzs5GTk7OQ59nZ2eDEAIXFxc0adIE9vb2lSG3s7OD2WyufAMwmUwoLS3FjRs3cOfOHbRu3RouLi5wdXV96L8uLi7o2LEjPD09YWfHf9Wkgv9LKIwl4EajERkZGTAajTAajSCEwNPTE25ubnBxcUG7du3g5+f3UEibNWv22McrKyvDH3/88ac3k7S0NOTk5ODSpUu4evUqvL294ePjAx8fH/j6+sLLy4u/ETBCWOyc4+vrS/kyXGHcu3cP+/fvx9GjRx8KuCVclqC5ubmBEMKkxoKCAhw7dqyyvoyMjD+9EQQFBaFjx45M6lMiQoiRUupb7dd46OUnKysLe/bsgcFgQEpKCvz8/BAQECCJgNdX1TeC9PR07Nu3D25ubtDr9QgNDYWPjw+fK7ACD73MUUpx/PhxGAwGxMbG4vLlyxg2bBj0ej2GDh2KFi1asC7RahUVFUhNTa38O965cwehoaHQ6/UIDAxEo0aNWJcoKzz0MpWeno6YmBjExsZCp9MhLCwMer0effv2Vfz5cGZmJmJjY2EwGHD8+HEEBgYiIiICo0aNgk6nY12e5PHQy0hxcTG+//57REdH49atW5gyZQpGjRqFp556SvJDdrHk5uYiPj4eMTExOH36NKZMmYKpU6eiQ4cOrEuTrNpCD0qpzT98fHwo97ALFy7QOXPmUCcnJxoSEkL37NlDTSYT67Ik58yZM3TWrFm0ZcuWNCwsjO7du5dWVFSwLktyAGTQGvLHZ0oYqqiogMFgQHBwMAICAmBnZ4ejR48iLi4Ow4cPh1arZV2i5Hh6eiIqKgpXrlzB8OHDMXfuXHh4eGDp0qXIy8tjXZ481PRuIOaH2nv68vJyunbtWtqhQwfq7+9Pv/rqK1pcXMy6LFkym830l19+oS+//DJ1cHCgM2fOpH/88QfrspiDmD09IaQ9IeQAIeQMIeQ0IWS2AO9FikQpxfbt29GtWzds3rwZW7ZsQWpqKl555RU0bNiQdXmyRAhBnz598PXXX+PcuXPQarXw8vLCe++9h7t377IuT5pqejeo7wcAFwA9H3zeDEAmAK/afkaNPX1SUhLt1asX7dGjB01ISKBms5l1SYp1+fJl+re//Y22bt2aLlu2jJaUlLAuyeYgZk9PKc2hlP764PMCAGcBuFnbrlL8+uuvGDJkCKZOnYo5c+YgIyMDQ4cOVe1MvC106tQJmzZtQmJiIpKSkuDh4YFNmzahoqKCdWnSUNO7wV/5ANAJwBUAzav52lQAGQAyOnToIP5bHWO//fYbjYiIoC4uLjQ6OpqWlpayLkm1Dh06RPv27Uu7detGY2NjWZdjE6ilpxcy8E0BGAGMqut7lTy8r6iooMuXL6dOTk508eLF9N69e6xL4uj9CT+DwUC7du1Kw8PD6Y0bN1iXJKraQi/IJTtCiD2AHQC+pZTuFKJNObp48SIGDBiA7du3IzU1FQsXLkSTJk1Yl8Xh/oRfaGgojh8/jg4dOqB79+7YsWMH67KYEGL2ngD4EsBZSulS60uSH7PZjKioKPj7+2P06NFITk6Gu7s767K4ajRq1AifffYZdu7ciQULFiAiIgI3b95kXZZNCdHT9wXwCoBAQsjxBx8hArQrC4/27rNnz+Z3h8lAnz59cOzYMXX2+jWN+8X8UMI5fdVz9+XLl/OloDL2yy+/UA8PD0Wd64MvwxVWbm4uhgwZwnt3haja6z/zzDM4dOgQ65JExX9TH9OpU6fQq1cv+Pj44ODBg/zcXSEs5/oxMTEYPXo01q1bx7ok0fDQP4bdu3dj4MCBWLx4MT755BN+Q4wCDRkyBCkpKVi6dClmzpyJ8vJy1iUJjoe+Hiil+PDDDzFz5kzEx8dj/PjxrEviROTu7o60tDRcvnwZQ4cORW5uLuuSBMVDX4fCwkKMHTsWe/bswdGjR9GrVy/WJXE20KJFCxgMBvj5+cHPzw+nTp1iXZJgeOhrceXKFTz33HNo2LAhDh48CBcXF9YlcTak1WoRGRmJf/3rXxg4cCB2797NuiRB8NDXwGg0wt/fHy+//DJiYmL4ra8qNn78eMTHx2PmzJn4/PPPWZdjNWXvrvgXpaamIiwsDGvXrsXIkSNZl8NJQK9evZCWloagoCAUFhbivffek+2dkjz0j/j5558xevRobNq0CcOGDWNdDichbm5uSE5OxqBBg1BSUoKPP/5YlsHnw/sqEhMT8eKLL+K7777jgeeq1aZNGxw4cAB79+7FnDlzLHeYygoP/QMHDhzAuHHjsHPnTgQFBbEuh5OwVq1aISkpCYcPH8abb74pu+Dz0AM4fPgwwsPDsW3bNjz33HOsy+FkwNHREXv37sWBAwewcOFC1uU8FtWf0x89ehQvvPACvv32Wzz//POsy+FkxNHRET/99BMGDhyIRo0aySb8qg792bNnERoaig0bNmDIkCGsy+FkyNnZGYmJiRgwYABatGiBf/zjH6xLqpNqQ5+fnw+9Xo/IyEiMGDGCdTmcjLVt2xY//fQT/P394eHhIfkORJXn9CaTCREREQgNDcWrr77KuhxOATp27IgtW7bglVdewYULF1iXUytVhn7u3LkghODTTz9lXQqnIP3798fixYsRFhaGO3fusC6nRqoL/caNG7Fnzx58//33in/cM2d7U6dOxcCBAzF+/HjJ7rOvqtCnpqbinXfewe7du+Ho6Mi6HE6hli9fjsLCQixYsIB1KdVSTeizsrIwevRoxMTEwNPTk3U5nILZ29tj27Zt2LZtGzZv3sy6nD9Rxfi2uLgYI0eOxOzZsxESopqNejmGWrVqhV27diEwMBDu7u6S2odBFT39O++8gyeffBJvv/0261I4FfH29saaNWswduxYFBYWsi6nkuJDn5ycjB07dmDVqlWyvCOKk7dRo0ahX79+mDdvHutSKik69IWFhZg0aRJWr17NJ+44ZpYvX44ffvgBBw8eZF0KAIWHfv78+ejbty9CQ0NZl8KpmKOjI1avXo3JkydLYpiv2NBbhvVRUVGsS+E4jBgxQjLDfEWGng/rOSmSyjBfkZfs+LC+esv2ZQIA0i7lwr+LU+Xr241X0c6xMQDAv4sT0i7d3+d9y2t9sGxf5kN/BoCINamVn3P1V3WYf+LECWaPMScsdv3w9fWlGRkZorSdnJyMl156CadOnVJ9Lx+xJrXy82NX8lFW8Xj/1joteehnNARootOioLQCbg4NcXheEPpGJqGdY2P4d3HCG4O7Cla7kk2YMAHNmzfHihUrRDsGIcRIKfWt9mtKCr3JZEL37t3x8ccfq3YX24g1qUj/bx4AwGzjf9pH3wj4aKB6+fn58PLyQnx8PHr06CHKMVQT+g0bNmDTpk04ePCg6q7JW4bhRy7nMa2j6uigd+eWyMovwt3icni5tgAAPiJ4YOXKlYiNjUVCQoIo7asi9CUlJejatSu2bNmCPn3U08N4v5+A5o3sce12CetS6q1ZAy28XFuo+g2grKwMnp6eWL9+PQYOHCh4+7WFXjETeStXrkTPnj1VE3jv9+/3EAWlFSgoleYtnDUpKK3Akct5OHI5r3JSMe1SrqpOB3Q6HRYvXox58+YhLS3NpiNTRVyyu3PnDj755BN89NFHrEuxia4L4mUZ9uocuZyHqKQLzE9LWBg7dixKS0vxww8/2PS4goSeELKBEHKDEMLk0Z6fffYZhg8fjm7durE4vM0s25eJLvPjHnsWXi76RibB+/0EeL+fUHl50fJfJdJoNFiyZAkWLFgAk8lku+MK1E4MgGCB2nosOTk5WLVqFRYtWsTi8DYVlXTB5jPytnTtdknlCCYq6QL6RiYhKun+fnN9I5MYVyeO4OBgtGnTBps2bbLZMQUJPaX0ZwBMxmeLFy/GxIkT0aFDBxaHt5lO8+JYl2BzlsnJLvPjZDVR+TgIIYiMjMQHH3yA4uJimxzTZuf0hJCphJAMQkjGzZs3BWnz0qVL2Lp1K+bPny9Ie1K0bF+mKgNflWV04/1+ArouiFfckN/f3x++vr6Ijo62yfEEu2RHCOkEYA+l9Om6vleoS3ZvvfUWNBqNone19X4/QRETdkLr3bkljl3JR+ZHytgJKSMjA6NHj8Zvv/0GrVZrdXu1XbKT7ex9cXExNm3ahGnTprEuRTQ88DU7cjkPZRUUy/ZlKqLn9/X1RZs2bfDjjz+KfizZhn7Lli3w8/NDly5dWJciiog1qbUGXq9JQYpuFi41eAkpulnQa1JsWJ10RCVdQFTSBSzbl1m5dkGuZsyYYZMhvlCX7L4DkArAgxCSRQiZLES7tYmOjsaMGTPEPgwzlvXz1dFrUhBpvx7tNLegIUA7zS1E2q9XbfCB++EvKK1A38gk2fb84eHhSE9Px2+//SbqcWS5DDc9PR3h4eG4ePGiIOc/UtN1QXyt1+JTdLPQTnPrT69nmVuhX9n/E7M02ZDrUt+3335bkKcvKe6cftWqVZg2bZoiAx+xJrXOxTeu5M+Bv/96rhglyZJlqa9l6C8X06ZNw8aNG0W9fCe70Ofl5eGHH37ApEmTWJfCTDZtVcPrTtW+rnZyCv4TTzwBX19fbNu2TbRjyC70MTExGDFiBJydnVmXIor6rEH/1BSOIqp76LUiqsOnpnCxypI9S/DlEH6xJ/RkdZcdpRSrV69GTEwM61JEUXWnm9oYzP2AcmCu3Va4klxkUyd8agq//zpXI8uSXgupnuuHhIRg5syZOH78OJ599lnB25dV6E+fPo2ysjLF3j57Jrv+jzc2mPvBUMZD/ldUDb8Ug6/VahEREYGdO3eKEnpZDe8NBgP0er1id8XhC3Fs69GeX0r0ej0MBoMobcsy9BwnlE7z4tBlvvTubejTpw+uXbuG33//XfC2ZRP6P/74A+fPn0f//v1ZlyIKpd46KgdmWv/5FFvRarUYPnw49uzZI3jbsgl9XFwchgwZAp1OV/c3y9Dd4nLWJajakct5kgt+aGioKEN82YReyUP7utbZc7Zx5HIeOs2Lk8xlvSFDhiA1NRV3794VtF1ZhL6oqAgHDhzAsGHDWJciiseZtefEJ5UJvmbNmqFv377Yu3evoO3KIvRJSUnw8fFBy5YtWZciCt7LS0/EmlRJ9PhizOLLIvRKHtoD9x8QwUmLVNbth4aGIj4+XtCNM2UR+kOHDiEwMJB1GaJxbtaAdQlcDVbsZzvUb9euHZydnXHmzBnB2pR86AsKCnD16lXFbm+9bF+mYjd9VAIzBfPJPR8fHxiNRsHak3zojx07Bm9vb9jZyWrFMKcwLCf3VBd6o9EIHx8f1mWIZrvxKusSuHpitXKPh57jGDFTNsHv0aMHTpw4IdhknixC7+tb7a4/isDP5+XFco5vS82bN0f79u1x9uxZQdqTdOgLCgpw5coVeHl5sS5FNPxynTzZemLPx8cHQj3eXdKhV8MkXgM7Sf8TcDWISrpg07X6Qp7XS/o3Tg3n86UmM+sSOBlQTehPnDghys4hUrFsX6ZiHzutBkcu59nslmjLZJ4QW9ZLOvRZWVmKfxotJ2/XbpfY5Py+efPmsLOzw+3bt61uS9Khz8nJgaurK+syOK5WaZds87wBV1dX5OTkWN2OpEOfnZ0NFxcX1mWIZkPKJdYlcALIyi+yyXFcXFyQnZ1tdTuSDX1JSQkKCwvh5KTcBzhM6qfMh2+qzbXbJTaZyXdxcVF2T5+Tk4O2bdsqdudbgC/BVZIjl/NEP7d3dXVVdk+fk5Oj6KE9pzxin9uroqfnk3icnIh9bq/4iTylT+IBwGif9qxL4ATUzrGxqO0rfiJPDcN7W13q4WxD7G20FT+8z8vLU/TMPcB3wVWa2UHu2PKaeM9ZdHJyQl5e3U81rotkQ28ymWBvb8+6DFE1b6Tsv5/aRCVdEHVZrr29vSD31AsSekJIMCHkPCHkIiFknhBtlpeXK/ruOk6a9JoUpOhm4VKDl5CimwW9JuWxfl7M/RHs7OxQXm79k5CsDj0hRAtgJYBhALwAjCOEWH0DvBp6+pw7fAMNKdFrUhBpvx7tNLegIUA7zS1E2q9/7OCLdb1eSj29H4CLlNJLlNIyAN8DCLO2UZPJpPie3qVFQ9YlcFXMtduKxqTsodcakzLMtdv6WO2INUGr1WphMpmsvtNOiNC7Aai6tCzrwWsPIYRMJYRkEEIybt68WWejQtxCKHViX+LhHo8ruVXD6/UPsZiTeUKtTrXZRB6ldC2l1JdS6uvs7Fzn99vZ2Qn6VA8p4rP30pJNW9Xwev2vIom5tNpkMkGr1VodfiFCfw1A1VUm7R68ZhWhzl+krLCMP8NOSj41haOIPvwo9CKqw6em8Hq3IeaCK6HmuYQIfToAd0JIZ0KIDsBYAFY/cU8NPT0/p5cWg7kf5pVPQZa5FcyUIMvcCvPKp8Bg7lfvNsTu6YWY57K6BUqpiRAyE8BeAFoAGyilp60uTKDLE1LWzrEx3wJbYgzmfjCU1T/kjxJznkaoy9iCTI9TSuMBxAvRlkWLFi0E2RpIymy1+QJnO2L+m96+fRstWrSwuh3JrsgTap2xlPHZe2Vxc2iIw/OCRGtfqPtRJBt6oW4j5DhbEXuTTKFuN5ds6IW6jVDK+PBeWdwcGuKNwV1Fa1+o280lHXre03NycrdY3IlnxQ/vLT29klfm8XN6ZRF7o9Ps7GxlD++bNm0Ke3t73LnDV61x0if20B5QQU8PKH8yb8trffhTaxWgWQOtqLP2FoqfyAPUMZnHn2Unf7Z6foHiJ/KA+z19VlYW6zJEY8tHHXPyVlxcjHv37gmyhZykQ//000/jxIkTrMsQjZj7qXG2MTvIXfRzeeD+E5y9vLyg0VgfWUmHXshnckuVmwO/6UaubDF5Z2E0GuHj4yNIW5IP/bFjx2A2m1mXIhq+ZZZ82fKSq2pC7+TkBEdHR1y8eJF1KaLht9fKk5tDQ5uenqkm9IA6hvicvOi0xCaX6CyKi4uRmZmJ7t27C9Ke5EPv6+uLjIwM1mVwHID7gc/8KMSmxzxx4gQ8PDzQsKEwo0LJh17pPb0tewzOOiwCDwg7tAdkEnolT+aJ/UxzTjjTn3+SyXFVF3o1TOZx0tasgRb/jRxus8tzj1Jd6AEgICAABw8eZF2GKFj9InH1d3JRMLNj5+bm4vLly4JN4gEyCX1oaChiY2NZlyGaZg20rEvgatC7c0umx4+Pj0dgYKBgk3iATEIfHByM5ORkFBUpc6cZL1frNzvkhKUhwH8jhzNfKm0wGKDX6wVtUxahd3R0hK+vLxITE1mXIgr/LtbfRMEJ6x+B7qxLQGlpKfbt24fhw4cL2q4sQg8Aer0eBoPVz9DguDr17txSEnMtycnJ8PLyQuvWrQVtVzaht5zXK/HS3RuDu0LD99KQBCkM6S3EGNoDMgr9E088AWdnZxw9epR1KaJoouOTeazNDmI/pLeglPLQA/eH+EqdxeeTeey4OTS02X3x9XXixAnodDp4enoK3rbsQq/U8/otr/Xh99bbmE5LKm+ekVLggf8N7YV6Jn1VgjzLzlb8/PyQm5uLM2fOwMvLi3U5ghN733Tuf1ito68PSim2bt2Kf//736K0L6ueXqPRYPLkyVi9ejXrUpjTa1KQopuFSw1eQopuFvSaFNYlycbsIHfJBh4ADh8+jPLycvTv31+U9mUVegCYOnUqvvnmG9y7d491KYKr73JPvSYFkfbr0U5zCxoCtNPcQqT9eh78epDauXt1oqOjMX36dFGG9oAMQ9++fXsMGDAAmzdvZl2KKOqz7HOu3VY0JmUPvdaYlGGu3VaxypK92UHusgj89evX8eOPP2LChAmiHUN2oQeAGTNmYOXKlYp85FV9rhG7kls1vJ4rdDmyV3VmXuqBB4Avv/wSo0ePhoODg2jHkNVEnkVQUBCKi4uRmpqKgIAA1uUIzs2hIa7drnnDzGzaCu2qCX425ct5q5JDz15VRUUFVq9ejV27dol6HFn29BqNBtOnT0d0dDTrUpj41BSOIqp76LUiqsOnpnBGFUlLswZa2QUeAOLi4uDm5oaePXuKehyrQk8IGUMIOU0IMRNCfIUqqj4mTpyIuLg43Lhxw5aHtYnD84JqXR1mMPfDvPIpyDK3gpkSZJlbYV75FBjM/WxYpbS4OTSsDPvJRcGyCzxwfwJvxowZoh+HWHNeTAjxBGAGsAbAW5TSeu1g6evrS4XY7HLy5Mlwd3fHvHnzrG5LiiLWpOLI5TzWZUhe784tJbNe/q+6ePEiAgICcOXKFUHunSeEGCml1XbEVvX0lNKzlNLz1rRhjZkzZ2LFihWKvc9+y2t9+I04tbDMyMs98ADw2WefYcqUKYJullETm53TE0KmEkIyCCEZN2/eFKTNHj16oE+fPlixYoUg7UlRr05sd26RIp2WyGpGvi6ZmZnYsWMH3nrrLZscr87hPSEkEUDbar60gFK6+8H3HASD4T0AnD9/Hv369UNmZiYcHR0FaVNqlu3LRFTSBdZlMNWsgRbNG9mjnWNjRfTsVYWHh6NHjx6YP3++YG3WNry36py+ygEOglHogfur9BwdHfHJJ58I1qYUdZoXx7oEJuQ4E19fGRkZCAsLw4ULF9C4sXDPxhPtnF4q3n//faxfvx7Xrl1jXYqoWG/SaEuzg9zRu3NLRQceAObPn4//+7//EzTwdbF29v4FACsAOAO4DeA4pXRoXT8ndE8PAO+88w7y8/Oxdu1aQduVmmX7MrHq4EWUVShvNaKbQ0OM9mkPQB1bgycmJmLGjBk4ffo07O3tBW1b9OH94xIj9Pn5+ejatStSUlLg4eEhaNtS1DcyqdZVe3Jh2f57Ur8uqgi6hdlshp+fH+bOnYvwcOEXVdUWelkuw62Oo6Mj3nzzTSxcuBDbtm1jXY7oDs8LwrJ9mUi7lItjV/Jl1/O7OTRU5KRcfe3YsQOUUowePdrmx1ZMTw8ARUVF6Nq1K7Zv3w5/f3/B25cy7/cTUFBawbqMP7HMQ5zJvgMv1xaV232rqVd/VGlpKby9vbFy5UoMHjxYlGOooqcHgMaNG+Pzzz/H5MmTYTQabbLQQSpOLgpGxJpU+Hdxwor9F2Cm969n23IE4ObQEDcLSjH9+Sex6uBF9OjgqNqevDaLFi1Ct27dRAt8XRTV0wP3txp68cUX4eHhgSVLlohyDDlYti8Tbwzuiog1qTh2JR/OzRrg2u0SNGugrdeIoLo3jN6dWyIrv+ihdqpeUeABr1t6ejpGjBiB//znP2jbtrrlL8JQxUReVdevX8czzzyD2NhY9OrVS7TjyE3VNwL/Lk5Iu5RbGWLLppx3i8vh5doCW17rg4g1qcjKv7/Eubrzb0t7XP2UlpaiZ8+eWLhwIcaNGyfqsVQXegD4/vvvsXjxYtUN8znpevfdd3Hu3Dns2LFDtK2wLBS/OKc6ERER8PDwwKJFi1iXwnFIT0/Hl19+iejoaNEDXxfFhp4QglWrVmHjxo1IT09nXQ6nYqWlpZg4cSKWL18u6nl8fSk29ADQpk0bLF++HBMnTkRJifwXsnDytGjRInh4eGDs2LGsSwGg8NAD/xvmv/fee6xL4VQoLS1NMsN6C8WHnhCCNWvWYOvWrdi6lW8RzdlOdnY2xowZg3Xr1kliWG+hqMU5NXF2dsauXbswePBguLu7o0ePHqxL4hSupKQEL7zwAqZNmybKk2etofie3uLZZ59FdHQ0Ro4cievXr7Muh1MwSimmTp2Kzp07491332Vdzp+ooqe3GDNmDE6ePIkXX3wR+/fvh06nq/uHOO4xLV26FKdOnUJKSopkzuOrUk1Pb/HBBx/A2dkZr7/+uiKfkMOxlZCQgC+++AK7du2y6cYYj0N1oddoNPj6669x5MgRrFy5knU5nIKcP38eEyZMwLZt29ChQwfW5dRIVcN7i6ZNm2L37t0ICAiAp6cngoKCWJfEydzt27eh1+uxZMkS9O3bl3U5tVJdT2/RuXNnfPfdd3jppZdw7Ngx1uVwMlZUVISRI0di2LBhmDRpEuty6qTa0APA888/j1WrViEkJASnTp1iXQ4nQ5ZLcx06dMDSpUtZl1MvqhzeVzVq1CiUlpZi6NCh2L9/vyr21+OEUVZWhjFjxsDBwQEbNmyARiOPPlT1oQeAcePGobS0FIMGDcKBAwfw5JNPsi6Jk7jy8nKMGzcOWq0W33zzDezs5BMl+VQqsokTJ6K8vBzPP/88EhMT8dRTT7EuiZOo0tJShIeHw2w2Y/v27YJvXy02Hvoq/v73v0On0yEwMBB79+6Ft7c365I4iSkuLsaoUaPQpEkTbN68WZYLvHjoHzFhwgQ0bNgQgwcPRnx8PHr27Mm6JE4iCgsLodfr0bZtW2zatElWQ/qq5DHzYGMRERFYvXo1goODERsby7ocTgKysrIwYMAAdO7cGV999ZVsAw/w0Ndo5MiRiI2NxfTp0xEZGcmX7KpYamoqevfujYiICKxbtw5arZZ1SVbhoa9F7969ceTIEezcuRPjx49HUVER65I4G9u4cSPCwsKwbt06vP3225K8geZx8dDXwc3NDcnJydBoNOjfvz+ysrJYl8TZgMlkwhtvvIElS5bg559/RkhICOuSBMNDXw+NGjXC119/jfDwcPTu3RupqamsS+JElJ+fj5CQEJw5cwZHjhxR3OVbHvp6IoRg7ty5WLduHcLCwrBx40bWJXEiOHv2LPz8/ODt7Y24uDg4OjqyLklwPPSPKSQkBMnJyViyZAkmT56MO3fusC6JEwClFF9++SX69++PBQsW4IsvvpD1DH1teOj/Ak9PTxiNRtjb28Pb2xsJCQmsS+KscPXqVQwbNgzR0dFISkrCxIkTWZckKh76v6hZs2ZYvXo1NmzYgGnTpvFeX4YsvXvPnj3Rr18/pKWloXv37qzLEh0PvZUGDRqEkydP8l5fZh7t3RcuXCi7NfR/lVWhJ4R8Rgg5Rwg5QQj5gRDiIFBdssJ7fflQa+9elbU9/T4AT1NKuwPIBDDf+pLk69Fef8eOHXwln4ScO3dOtb17VVaFnlL6E6XU9OCPaQDaWV+SvFl6/ZiYGHz44Yfw9/fHgQMHWJelallZWZgyZQqee+45DBo0SJW9e1VCntNPAvBjTV8khEwlhGQQQjJu3rwp4GGlKTAwEEajEf/85z8xZcoUDB06lO/FZ2N5eXl4++238cwzz8DZ2RmZmZl46623VNm7V1Vn6AkhiYSQU9V8hFX5ngUATAC+rakdSulaSqkvpdTX2dlZmOolTqPRYNy4cTh79izCwsIQEhKCcePG4eLFi6xLU7TCwkIsWbIEHh4eKCgowMmTJ7FkyRJFLrT5SyilVn0AmAggFUDj+v6Mj48PVaOCggK6ePFi6uTkRKdPn05zcnJYl6QoZWVldNWqVdTV1ZWOGTOGnj9/nnVJzADIoDXkz9rZ+2AAcwHoKaX8FrQ6NG3aFAsXLsS5c+fQqFEjeHp64tVXX0V6ejrr0mTt+vXr+Oijj/DEE09g+/bt2L17N7Zu3YquXbuyLk2SrD2n/zeAZgD2EUKOE0JWC1CT4rVq1QpffPEFMjMz4enpifDwcPTq1QsbN25EcXEx6/JkgVKKQ4cOYdy4cXjqqafw+++/Y/fu3UhMTISvry/r8qStpiGAmB9qHd7XxGQy0T179tCQkBDq5ORE58yZQy9cuMC6LEm6e/cujY6Opk8//TT18PCgUVFRND8/n3VZkgOxhvecMLRaLYYPH464uDgcPXoUdnZ2CAgIQHBwMHbu3Kn6zTvMZjPS09Px+uuvo2PHjkhKSkJUVBTOnj2LWbNmwcHBgXWJskIog8Ujvr6+NCMjw+bHlZOSkhJs27YNMTExSE9Px4ABA6DX6zFixAi4uLiwLk90xcXF2L9/PwwGA2JjY9GiRQuEh4dj6tSpcHNzY12e5BFCjJTSas9zeOhlID8/HwkJCTAYDEhISIC7uzv0ej30ej28vb0VsYUTcH9CLi4uDrGxsdi/fz+effZZ6PV6hIaG8km5x8RDryDl5eU4dOgQDAYDdu/eDUopRowYgYCAAPj4+MDd3V02j1e6ceMGjEYj0tPTkZCQgDNnzmDIkCEIDQ1FSEgInJycWJcoWzz0CkUpxenTpxEfH4+jR4/CaDQiNzcXPXr0gK+vL3x8fCTzRmAJuOUjIyMD9+7dQ8+ePeHj44NBgwZhwIABaNCgAdM6lYKHXkVyc3Px66+/Vgar6huBl5cXXF1d4eLiAhcXl8rPnZ2drdrWmVKKu3fvIjs7Gzk5OcjJyan8/NKlSzAajQ8F3PKG1KVLF8WcmkgND73KWd4IMjMzHwqk5b+3b9+Gs7Nz5ZtBkyZNYG9vDzs7u8oPs9mM8vJymEwmmEwmlJaW4saNG5Xt2NnZ/enNxMXFBR07dkTPnj15wG2Mh56rVVlZGa5fv17ZSxcVFVWG2xJ0rVZb+QZgb28Pe3t7tG7dujLgTZs2Zf3X4KqoLfTK3PmPeyw6nQ7t27dH+/btWZfC2YA8pnk5jhMMDz3HqQwPPcepDA89x6kMDz3HqQwPPcepDA89x6kMDz3HqQyTFXmEkJsAfheh6VYAbonQrljkVC+vVRxi1dqRUlrtttNMQi8WQkhGTUsPpUhO9fJaxcGiVj685ziV4aHnOJVRWujXsi7gMcmpXl6rOGxeq6LO6TmOq5vSenqO4+rAQ89xKqO40BNCPiOEnCOEnCCE/EAIcWBdU00IIWMIIacJIWZCiCQvMRFCggkh5wkhFwkh81jXUxtCyAZCyA1CyCnWtdSGENKeEHKAEHLmwb//bFseX3GhB7APwNOU0u4AMgHMZ1xPbU4BGAXgZ9aFVIcQogWwEsAwAF4AxhFCvNhWVasYAMGsi6gHE4A3KaVeAPwBvG7L/6+KCz2l9CdKqenBH9MAtGNZT20opWcppedZ11ELPwAXKaWXKKVlAL4HEMa4phpRSn8GkMe6jrpQSnMopb8++LwAwFkANntsj+JC/4hJAH5kXYSMuQG4WuXPWbDhL6caEEI6AegB4IitjinLjTEJIYkA2lbzpQWU0t0PvmcB7g+jvrVlbY+qT62cOhFCmgLYAeCflNK7tjquLENPKR1U29cJIRMBjAAQRBkvRKirVom7BqDqFrntHrzGWYkQYo/7gf+WUrrTlsdW3PCeEBIMYC4APaVU3c94tl46AHdCSGdCiA7AWAAGxjXJHrn/1I8vAZyllC619fEVF3oA/wbQDMA+QshxQshq1gXVhBDyAiEkC0AfAHGEkL2sa6rqwYToTAB7cX+yaSul9DTbqmpGCPkOQCoAD0JIFiFkMuuaatAXwCsAAh/8jh4nhITY6uB8GS7HqYwSe3qO42rBQ89xKsNDz3Eqw0PPcSrDQ89xKsNDz3Eqw0PPcSrz/wGllpoafsGzOwAAAABJRU5ErkJggg==\n",
"text/plain": [
"