From 7800272033933871d4b823d74209352cbf242be3 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 15 Feb 2024 14:42:47 -0500 Subject: [PATCH] add eigenproblem examples --- notes/Eigenproblems.ipynb | 1299 +++++++++++++++++++++++++++++++++++++ 1 file changed, 1299 insertions(+) create mode 100644 notes/Eigenproblems.ipynb diff --git a/notes/Eigenproblems.ipynb b/notes/Eigenproblems.ipynb new file mode 100644 index 0000000..af54aa5 --- /dev/null +++ b/notes/Eigenproblems.ipynb @@ -0,0 +1,1299 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4f446c8e", + "metadata": {}, + "source": [ + "# Eigenvalue problems\n", + "\n", + "For a square matrix $A$, an **eigenvector** $x$ is a solution to:\n", + "\n", + "$$\n", + "Ax = \\lambda x \\qquad \\Longleftrightarrow \\qquad (A - \\lambda I)x = 0\n", + "$$\n", + "for some $x\\ne 0$ and some scalar $\\lambda$, the **eigenvalue**.\n", + "\n", + "The important thing about eigenvectors is that they make a complicated matrix $A$ **act like a scalar** $\\lambda$, which makes it easier to understand.\n", + "\n", + "For a typical $m \\times m$ matrix $A$, you will almost always have $m$ eigenvalues $\\lambda_1, \\ldots, \\lambda_m$ and corresponding linearly independent eigenvectors $x_1, \\ldots, x_k$. That allows you to form a **basis of eigenvectors**: you can write *any* vector as a linear combination of eigenvectors, and then anything you want to do with $A$ becomes *easy (scalar) on each term*." + ] + }, + { + "cell_type": "markdown", + "id": "8112d161", + "metadata": {}, + "source": [ + "## Computing eigenvalues and eigenvectors\n", + "\n", + "From the second form of the eigen-equation above, we immediately see that for any eigenvalue $\\lambda$, the matrix $A - \\lambda I$ must be **singular**, which implies that its **determinant is zero**:\n", + "\n", + "$$\n", + "\\det(A - \\lambda I) = 0\n", + "$$\n", + "\n", + "If you think a little bit about determinant formulas, which are a sum of terms that multiply one entry from each row and column of the matrix, you'll quickly see that $\\det(A - \\lambda I)$ is a **polynomial of degree m in λ**, called the **characteristic polynomial of A**.\n", + "\n", + "So, one way to find eigenvalues, which you probably already learned, is to compute this determinant to find the polynomial, and then **calculate the roots of the polyomial**. Once you have an eigenvalue, finding an eigenvector is relatively easy: it comes from the nullspace of $A - \\lambda I$, and can be computed by Gaussian elimination.\n", + "\n", + "However, while this is an easy procedure for $2\\times 2$ matrices, it is not so obvious how to carry it out for large matrices, where you would have to find **roots of high-degree polynomials**. Efficient and accurate solution of eigenproblems is a famously tricky problem in numerical linear algebra, and the modern algorithms are *not obvious*.\n", + "\n", + "Here, we will cheat and use the function `eigen` in Julia's `LinearAlgebra` standard library:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "02daa5e0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n", + "values:\n", + "2-element Vector{Float64}:\n", + " -0.3722813232690143\n", + " 5.372281323269014\n", + "vectors:\n", + "2×2 Matrix{Float64}:\n", + " -0.824565 -0.415974\n", + " 0.565767 -0.909377" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using LinearAlgebra\n", + "\n", + "A = [1.0 2.0\n", + " 3.0 4.0]\n", + "λ, X = eigen(A)" + ] + }, + { + "cell_type": "markdown", + "id": "a16f128e", + "metadata": {}, + "source": [ + "Here, our 2x2 matrix has two eigenvalues:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "dd5fc376", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2-element Vector{Float64}:\n", + " -0.3722813232690143\n", + " 5.372281323269014" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "λ" + ] + }, + { + "cell_type": "markdown", + "id": "490048cd", + "metadata": {}, + "source": [ + "and corresponding eigenvalues given by the *columns* of the matrix $X$:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "2808625b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2×2 Matrix{Float64}:\n", + " -0.824565 -0.415974\n", + " 0.565767 -0.909377" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X" + ] + }, + { + "cell_type": "markdown", + "id": "5ff64dcd", + "metadata": {}, + "source": [ + "Let's check that this works fore the first column, denoted `X[:,1]` in Julia, by computing $Ax - \\lambda x$:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "d93ec2e4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2-element Vector{Float64}:\n", + " 0.0\n", + " 5.551115123125783e-17" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A * X[:,1] - λ[1] * X[:,1]" + ] + }, + { + "cell_type": "markdown", + "id": "4459722c", + "metadata": {}, + "source": [ + "Yup, it's zero up to roundoff errors! (The computer only does arithmetic to 15–16 digits by default.)" + ] + }, + { + "cell_type": "markdown", + "id": "8e617b0e", + "metadata": {}, + "source": [ + "These eigenvectors will form a basis for $\\mathbb{R}^2$, but they are not orthogonal:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "515895b8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2×2 Matrix{Float64}:\n", + " 1.0 -0.171499\n", + " -0.171499 1.0" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X' * X" + ] + }, + { + "cell_type": "markdown", + "id": "ae9efc5f", + "metadata": {}, + "source": [ + "The off-diagonal elements are the dot products of the eigenvectors with each other, and are not even close to zero: they are **not orthogonal**." + ] + }, + { + "cell_type": "markdown", + "id": "e8911a41", + "metadata": {}, + "source": [ + "## Using eigenvalues and eigenvectors: Matrix powers & exponentials\n", + "\n", + "The key point is the **matrix acts like a scalar λ** when acting on an eigenvector,\n", + "for **any operation you might want to do with the matrix**.\n", + "\n", + "For example:\n", + "\n", + "* Matrix *inverses*: $A^{-1} x = \\lambda^{-1} x$\n", + "* General matrix *powers*: $A^n x = \\lambda^n x$\n", + " - These show up in lots of problems, including linear *recurrence equations*.\n", + "* Matrix *exponentials*: $e^{A}x = \\left(I + A + A^2/2 + A^3/6 + \\cdots + A^n/n! + \\cdots\\right) x = e^{\\lambda} x$\n", + " - These are important when solving *systems of differential equations* $\\frac{dx}{dt} = Ax$: for any initial condition $x(0)$, the solution is $x(t) = e^{At} x(0)$.\n", + " \n", + "For *any* $x$, you expand in the basis of eigenvectors:\n", + "$$\n", + "x = c_1 x_1 + c_2 x_2 + \\cdots + c_m x_m = \\sum_k c_k x_k = X c\n", + "$$\n", + "for some coefficients $c = X^{-1} x$, and then act $A$ on *each term as a scalar*. For example,\n", + "with matrix powers:\n", + "$$\n", + "A^n x = \\sum_k c_k \\boxed{ \\lambda_k^m} x_k\n", + "$$\n", + "multiplies *each term* by $\\lambda_k^m$.\n", + "\n", + "This means, for example, that multiplying a random vector by $A$ many times tends towards the largest eigenvector. For example, with our 2x2 matrix above:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "f291f3be", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2-element Vector{Float64}:\n", + " 1.1064255702282437\n", + " 0.6604584519709495" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = randn(2) # a random 2-component vector" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "48a9a3e4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2-element Vector{Float64}:\n", + " 5.1264936878916495e72\n", + " 1.1207236302712575e73" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y = A^100 * x # should be parallel to x₂" + ] + }, + { + "cell_type": "markdown", + "id": "08235ffe", + "metadata": {}, + "source": [ + "This should be nearly parallel to $x_2$. Recall that the eigenvalues were $\\lambda_1 \\approx -0.37$ and $\\lambda_2 \\approx 5.4$, so:\n", + "$$\n", + "x = c_1 x_1 + c_2 x_2 \\implies A^{100} x = (-0.37)^{100} c_1 x_1 + (5.4)^{100} c_2 x_2 \\approx (5.4)^{100} c_2 x_2 \n", + "$$\n", + "since the second term is exponentially larger than the first.\n", + "\n", + "Here, we can check that it is parallel to $x_2$ by looking at the ratio of the components:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "250978f6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.186140661634507" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y[2]/y[1]" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "eb4c1420", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.1861406616345076" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x₂ = X[:,2] # second eigenvector\n", + "x₂[2]/x₂[1]" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "f7cac0ec", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "3-element Vector{Float64}:\n", + " 0.0\n", + " -0.0\n", + " -6.366701349055913e-17" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cross([y; 0], [x₂; 0]) / norm(y)" + ] + }, + { + "cell_type": "markdown", + "id": "3bc89010", + "metadata": {}, + "source": [ + "Yup, they are parallel!\n", + "\n", + "It should also be that $\\Vert A^n x \\Vert$ grows nearly exponentially, proportional to $\\lambda_2^n$. Let's plot it:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "6f321d2f", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGwCAYAAABy28W7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABaVUlEQVR4nO3deVhV9drG8e8GBMQUc8IJh8rKGUUzTSobKCzHLOeh1LLwOFCWU5nmUJZmBWrmSS1TKUWzsoy3SdJOKYI5VZomahhpBoIKulnvHz9FcShRYLE39+e69tVZa69388g5b/v2Nz0Oy7IsRERERIoID7sLEBERETmbwomIiIgUKQonIiIiUqQonIiIiEiRonAiIiIiRYrCiYiIiBQpCiciIiJSpHjZXUBeZWdn8/vvv1O6dGkcDofd5YiIiMglsCyLI0eOULVqVTw8/nlsxOXCye+//05gYKDdZYiIiMhl2Lt3L9WrV//HZ1wunJQuXRowf7gyZcrYXI2IiIhcirS0NAIDA3O+x/+Jy4STqKgooqKicDqdAJQpU0bhRERExMVcypIMh6v11klLS8Pf35/U1FSFExEREReRl+9v7dYRERGRIkXhRERERIoUl1lzkldOp5MTJ07YXYbb8fb2/tctYCIiIlfCZcLJuQtiL8ayLA4cOMDff/9dOIUVMx4eHtSuXRtvb2+7SxERETfldgtik5OT+fvvv6lUqRJ+fn46qC0fnT4Ar0SJEtSoUUO/WxERuWR5WRDrMiMnl8LpdOYEk/Lly9tdjluqWLEiv//+OydPnqREiRJ2lyMiIm7IrRYPnF5j4ufnZ3Ml7uv0dM6/Ta+JiIhcLrcKJ6dpuqHg6HcrIiIFzWWmdS51QayIiIhcHqcT4uIgORmqVIGQEPD0LPw6XGbkJDw8nG3btrF+/Xq7SxEREXE7MTFQqxa0aQM9eph/1qpl7hc2lwknIiIiUjBiYqBLF9i3L/f9/fvN/cIOKAonRcTXX39NrVq1ALOu47fffsvzZ4SFhTFu3DhuvvlmatasybZt2wr154uIiOtxOmHoULjQwSKn7w0bZp4rLAonF+B0wtdfw+LF5p9FZZlLQkICJUqUICQk5ILvb9myhdq1a/O///2PgQMH8tFHHxVyhSIi4mri4s4fMTmbZcHevea5wqJwco6iNOd2riFDhvDUU0+xadMmzj07LzU1lRIlStCvXz/AbPktW7Zs4RcpIiIuJTk5f5/LDwonZylqc25nW7RoEVdffTXh4eEcOXKEXbt25Xp/y5Yt3HTTTbmu69evD8DixYvx9fVl//79Oe8PGDCARo0akZqaWjh/ABERKXIsCzZuzH3Piwv3patSpRAKOsXtw4llQUbGv7/S0mDIkH+ecxs61Dx3KZ+Xn00BMjIyGD16NC+99BLVq1fH39+fxMTEXM9s2bKFhg0b5lxv3ryZBg0aANCtWzduuOEGpkyZAsD48eNZvXo1n376Kf7+/vlXqIiIuIy0NOjZE1555cy95vzANupxK9/k3HM4IDDQbCsuLG5/zsnRo3DVVVf+8y3LjKhc6nd5ejqUKnXlPxdg0qRJ3HvvvdStWxeAevXqkZiYyAMPPJDzzNatW7nrrrsAOHnyJOnp6TnTOg6Hg0mTJtGlSxeqVq3Ka6+9RlxcHNWqVcufAkVExKVs2ADdusGvv5pzTHp0yybgvelMZhQlOMlExnIra3IO3pwxo3DPO3GZcBIeHk54eHhO46DiYteuXcyZM4ctW7bk3GvQoMF5Iyevv/56zn/28vJix44dud6///77qVevHuPHj+fzzz/PmfIREZHiw7Lgtdfg6afhxAmoWRM+mPknzSP7Ap8C8D4P8ihzAAfVq5tg0rlz4dbpMuHkcvn5mVGMf7NmDbRt++/PrVoFt956aT83PwwfPpxDhw5RvXr1nHvZ2dl5HvVYvXo1P/30E06nk4CAgPwpTkREXMahQ9CvH3z8sbnu3Bnm9/ua0gN7wu+/g68v2dNnUOnGR5l1wGHrCbFuH04cjkubXgkNherVzeLXC60XcTjM+6GhhfdfVGxsLGvXriUhIQEvrzP/Va1fv55HHnmEQ4cOXVL35Y0bN/Lggw/y5ptvsmTJEp599lk++OCDgixdRESKkLg46N7dfMf5+MCrrzgZlDIBR4cXzJfejTfC++/j0bAht9tdLMUgnFwqT08z1NWliwkiZweU073uCnPO7eTJkwwdOpQRI0YQFBSU670yZcoAkJiYyJ133vmPn/Pbb79x3333MXLkSHr37k29evVo3rw58fHxBAcHF1T5IiJSBDidMHkyPP88ZGfD9ddDzBv7qT+ph5kyAHj4YXjjjfxbKJkP3H63Tl507gxLl8K5MybVq5v7hTnn9tFHH3Ho0CEGDx583nuBgYH4+fmdt+7kXH/99RdhYWG0b9+e0aNHAxAcHEy7du0YM2ZMQZQtIiJFRHKyGe1/7jkTTPr0gcRJn1C/R2MTTK66ChYuhLffLlLBBDRycp7OnaFDB/u7Mnbq1Inhw4df8D2Hw0FGRsa/fka5cuXYvn37efc//PDDK65PRESKrs8+M2Hkzz9N7pj9eha9to6CB6ebB5o2hSVLoE4dewu9CIWTC/D0hNtvt7sKERGRvDlxAsaOhalTzXXjxhDz8q9cM7qb2T8M5lCvqVPN4pMiSuFERETEDfz2mzm75PvvzXV4OExvEY33AwPhyBG4+mqYPx/at7ezzEtS6GtOfv75Z4KCgnJeJUuWZMWKFYVdRpFTq1Ythg0bBsC4ceMKvS+O3T9fREQu37JlEBRkgknZsrBi0VEisx7Fu083E0xuuQU2bXKJYALgsM7tIFeI0tPTqVWrFnv27KHUJS7GOX0IW2pqas6uldOOHz/O7t27qV27Nr6+vgVRcrGn37GISNFx/DhERMCsWea6ZUv44PmtVIvoClu3mu2mo0eb7Tpe9k6W/NP397lsrXTlypXceeedlxxMRERExPjpJ+jaFX780VyPfMbihdr/xavjEDh2DAICzG6cU61NXEmep3XWrFlDu3btqFq1Kg6H44JTMjNnzsz5m3VwcDBxcXEX/Kz333+frl275rloERGR4mzBAggONsGkUiX4v5g0puzpgdeggSaYhIaaaRwXDCZwGeEkIyODxo0bExkZecH3o6OjGTZsGGPGjCEhIYGQkBDCwsJISkrK9VxaWhpr166l7b+cGZ+ZmUlaWlqul4iISHF05IjZItyvn2lse+edsHXBBu4ccWprsKcnvPgifPqpGTlxUXme1gkLCyMsLOyi70+fPp3+/fszYMAAAGbMmMHq1auZNWsWU6ZMyXnuww8/5J577vnXdQtTpkxh/PjxeS1TRETErSQmwkMPwY4dJoNMGG8xsuRreLQ/1cWvRg0TUFq2tLvUK5avu3WysrKIj48nNDQ01/3Q0FDWrVuX696lTumMGjWK1NTUnNfevXvzs2QREZEizbIgMhJatDDBJDAQvv3wEKP/1x6PJ4ebYNKpk0kvbhBMIJ8XxB48ePCCXW8DAgI4cOBAznVqaio//PADy5Yt+9fP9PHxwcfHh6ioKKKionA6nflZsoiISJF1+DD07w/Ll5vr9u3h3UfjKPPYqS5+3t4wfTo88cSZRnBuoEDOOXGc8wuyLCvXPX9/f/744w+8vb0v+TPDw8PZtm0b69evz7c6RUREiqp168zZJcuXmwzy+qtOVgS/QJn2t5tgcv315mCT8HC3CiaQz+GkQoUKeHp65holAUhJSTlvNCWvoqKicjrqioiIuKvsbJgyBW69FZKS4LrrYMPK3/nPyrtxjDuri198vEkvbihfw4m3tzfBwcHExsbmuh8bG0urVq2u6LM1ciIiIu7ujz/g3nvNuWlOJ/ToAZte+oyGvYPgq69MF78FC8zrqqvsLrfA5DmcpKenk5iYSGJiIgC7d+8mMTExZ6twREQEc+fO5e2332b79u0MHz6cpKQkBg0alK+Fu6OxY8fi4+NDjx49Lvh+WFgY48aN4+abb6ZmzZps27atkCsUEZGCEhtrGvXFxoKfH8ybc4KFVZ/G74Ew0164cWMzWtKnj92lFjwrj7766isLOO/Vt2/fnGeioqKsmjVrWt7e3lbTpk2tb775Jq8/5jyRkZFW3bp1reuvv94CrNTU1POeOXbsmLVt2zbr2LFjV/zz7JCammpFRkZagLVjx47z3q9evbo1b948y7Is64UXXrBefPHFQq7Q9X/HIiJFzYkTljVqlGU5HJYFltWggWXt+HyXZbVoYW6AZYWHW5aL/3s3NTX1ot/f57K1t87lcPfeOseOHeOqq65i6dKldOrUKed+amoqTZo0YdeuXQBMnToVf39/HnvssUKtzx1+xyIiRUVSEnTvbha/Ajz2GLx261J8nhgAqammi99//wudO9taZ37IS2+dQu9KLP/s5MmT+Pn5sWXLllz3t2zZwk033ZTrun79+oVdnoiI5JMVK8x61nXroEwZWLrwGLMdj+PT80ETTG6+GRIS3CKY5JXLhJPisltn7NixpKenXzCcNGzYMOd68+bNNGjQoLDLExGRK5SZCUOGmHPTDh+G5s1hy9KfeGDqzTB7tnnomWdgzRqoVcvWWu1ib//kPAgPDyc8PDxnWOiSWZZpQFDY/PzyvO88Pj6e2bNnc999950XTrZu3cpdpxo4nTx5kvT0dMqWLZtf1YqISD5yOiEuDpKToUoVCAkxR87/8gt062YGRACeetJiyo0L8OoYbr6rKlaEd9+Fe+6x9w9gM5cJJ5ft6FF7tlulp5stX5coOzubxx57jMGDB9OiRQt69uxJVlZWzkF1r7/+es6zXl5e7NixI99LFhGRKxcTA0OHwr59Z+5Vr25mZ95+23w9VKgAC2cd4Z4Pn4BpC81Dd95pgkmVKvYUXoS4zLSOu3vjjTf4888/mTBhAg0bNuTkyZP8/PPPAOzdu5fbb7+devXq0ahRIz744AObqxURkQuJiYEuXXIHEzDXr79ugsntt8PWhQncMzoYFi4EDw+YOBFWr1YwOcVlRk4uu7eOn5/5X0Nh8/O75Ef379/Ps88+y+LFiylVqhR16tTBx8cnZ52Jl5cXM2bMICgoiJSUFJo2bUrbtm0plYeRGRERKVhOpxkx+ac9sGVKW8S2j8Sr/VOQlWWGVBYvhtatC69QF+Ay4eSy15w4HHmaXrHDkCFDCAsL47777gPMtE3dunVz1p1UqVKFKqfSdKVKlShXrhx//fWXwomISBESF3f+iMnZruYv/nukP14RK8yN9u3NPE/58oVSnytxmXDirj7++GO+/PJLtm/fnut+w4YNz1sUC7Bhwways7MJDAwsrBJFROQSJCdf/L1WrGURPahJEk4vbzynvQz/+Y/bNezLLwonNrv//vs5fPjwefffeeed8+4dOnSIPn36MHfu3MIoTURE8uBCy0UcZDOSF5nAc3jhZAfXkRoVTbNHmxZ+gS7EZRbEFpdzTi4mMzOTTp06MWrUqCtuoigiIvkrOxu+/z73vQAOsJp7mMwYvHDyHj1oX20jTformPwblwknxbkrsWVZ9OvXjzvuuIPevXvbXY6IiJwlJQXuuw9Gjjxz725i2URj7ub/yMCPR3ib3ixk0uul8fS0r1ZX4TLhpDhbu3Yt0dHRrFixgqCgIIKCgti8ebPdZYmIFHtffmmOoP/sM/D1hbdmnmB7p9F8xj0EkMJmGtCMDfxf4MMsXeYojifRXxatOXEBrVu3Jjs72+4yRETklJMnYcIEczyJZUG9ehDz6h5uGN8jp4vf7+0eY/sDrzKrZsmcE2Ll0iiciIiI5MG+fdCjh9k6DNC/P0TdvQKfrg/D33+bLn5z51L1wQd5yNZKXZemdURERC7Rxx+baZy4OChdGpbMP87ckv/Bp1snE0yaNzeNcx580O5SXZrLhJPivltHRETsk5UFw4dDu3Zw6BAEB8PmZb/QdUZLiIw0Dz31FHz7LVxzjb3FugGXmda57BNiRURErsCvv0LXrhAfb66HDYOpjRZSotMgyMgwXfwWLIC2bW2t0524zMiJiIhIYVuyBJo0McGkXDn4JDqdVw/3o8QjvU0wuf12SExUMMlnbhlOrH/quiRXRL9bESkOjh6FgQOhe3c4cgRCQmDbkh9p+1wzM0ri4QHjx8P//R9Uq2Z3uW7HZaZ1LkWJEiUAOHr0KCVLlrS5GveUlZUFgKf2xImIm9qyxUzjbNtmWt+MHWMxLmA2nu2GQ2YmVK0KixbBbbfZXarbcqtw4unpSdmyZUlJSQHAz88Ph5oq5Zvs7Gz+/PNP/Pz88PJyq//piIhgWTB3LgwZAsePQ+XKEP3m39z6zgCYuMw8dN99MH++WWciBcbtvmEqV64MkBNQJH95eHhQo0YNhT4RcSupqfDYYxAdba7vuQcWDfkf5Z7oBnv2QIkS8OKLZsuO/v1X4NwunDgcDqpUqUKlSpU4ceKE3eW4HW9vbzw83HKpkogUU+vXQ7dusGsXeHnB5InZPMkreHQYY46CveYaszJWR1kUGpcJJ1FRUURFReF0Oi/peU9PT62LEBGRi7IsePVV07DvxAmoVQs+iEqh2et9YPVq81DXrvDmm6AjLAqVw3Kx7RenzzlJTU2lTJkydpcjIiIu6OBB6NcPPvnEXD/wAMzv8yVXPdYTDhwwXfxefx0GDNA0Tj7Jy/e3y4yciIiI5Ic1a0xvnP37wccHZrxykscOjMfRcdKZLn7R0dCggd2lFltaPCAiIsWC02k6CbdpY4LJDTfAxpX7GPT+HTgmnWov3L+/WYSiYGIrjZyIiIjb+/136NkTvv7aXPfrB7PafoRv937w11+mi9+bb5pT18R2GjkRERG39umn0LixCSalSsF7b2cyr+xwfB9qb4JJcDBs3KhgUoQonIiIiFvKyjKNgtu2NQtgg4Jg8/Kd9Ii6BWbMMA8NGwZr18J119lYqZxL0zoiIuJ2du0yAyE//GCuBw+Gac0W4/3AY6ZZTrly5qTXdu1srVMuzJaRk927d9OmTRvq1atHw4YNycjIsKMMERFxQx98YDoJ//ADXH01rFxylDeODcC7Xw8TTFq3Np2EFUyKLFtGTvr168fEiRMJCQnhr7/+wsfHx44yRETEjRw7Zk6Xf/NNc92qFSx9fgtVhp3VxW/MGBg3zhwFK0VWof+3s3XrVkqUKEFISAgA5cqVK+wSRETEzWzfbg5z3bzZZJBRIy3G15iLV/uzuvi99x7ccYfdpcolyPO0zpo1a2jXrh1Vq1bF4XCwYsWK856ZOXMmtWvXxtfXl+DgYOLi4nLe27FjB1dddRXt27enadOmTJ48+Yr+ACIiUnxZFrz9NjRrZoJJQAB8EZPKpF3d8Xr8URNM7rkHNm1SMHEheQ4nGRkZNG7cmMjIyAu+Hx0dzbBhwxgzZgwJCQmEhIQQFhZGUlISACdOnCAuLo6oqCi+++47YmNjiY2NvbI/hYiIuC2n02wDXrzY/PN0i7UjR6BXL3Nu2tGjcPfdsGX+eto82dSc8OrlBS+9BKtWQaVKdv4RJI/yPK0TFhZGWFjYRd+fPn06/fv3Z8CAAQDMmDGD1atXM2vWLKZMmUL16tVp3rw5gYGBALRt25bExETuvvvuC35eZmYmmZmZOddpaWl5LVlERFxUTAwMHQr79p25V726uffmm7BzJ3h6wgsTLJ7xeRWP9qe6+NWsaToJ33yzfcXLZcvX3TpZWVnEx8cTGhqa635oaCjr1q0DoHnz5vzxxx8cPnyY7Oxs1qxZQ926dS/6mVOmTMHf3z/ndTrUiIiIe4uJgS5dcgcTMNcjRphgEhgIaz88yKh17fB46kkTTDp3hoQEBRMXlq/h5ODBgzidTgICAnLdDwgI4MCBAwB4eXkxefJkbr31Vho1akSdOnW4//77L/qZo0aNIjU1Nee1d+/e/CxZRESKIKfTjI5Y1sWfKVkSfoxcQ4vHgkx7YR8fmDkTli41e4jFZRXIbh3HOe2lLcvKde/fpobO5uPjg4+PD1FRUURFReE8PdkoIiJuKy7u/BGTs3ng5Kljk/DvNB6ys00Xv+hoc069uLx8HTmpUKECnp6eOaMkp6WkpJw3mpJX4eHhbNu2jfXr11/R54iISNGXnHzx96rwO//HXUxgHI7sbNPFLz5ewcSN5Gs48fb2Jjg4+LzdN7GxsbRq1eqKPjsqKop69erRvHnzK/ocEREp+qpUufD9e/mUTTSmDV+TTim2j34X5s0zHf3EbeR5Wic9PZ2dO3fmXO/evZvExETKlStHjRo1iIiIoHfv3jRr1oyWLVsyZ84ckpKSGDRo0BUVGh4eTnh4OGlpafj7+1/RZ4mISNF29Ch4eJgZG4ASZDGJMYzgFQASCGJY5Wi+nHC9jVVKQclzONmwYQNt2rTJuY6IiACgb9++zJ8/n65du3Lo0CEmTJhAcnIyDRo0YNWqVdSsWTP/qhYREbd04gQ8+6w5nuS0a9jFIrrTAtPFL5LBjOBl3ovyxdPTpkKlQDks65/WQhcdZy+I/eWXX0hNTaVMmTJ2lyUiIvlkzx7TSfi778z1449DjxIf0OiNAZSx0jhMWR7hbeIDOzFjhtkxLK7j9MzHpXx/u0w4OS0vfzgREXENMTHmpNe//wZ/f3g76hid48508Uut34qvH12Ef6OahISgERMXlJfvb7VlFBER2xw/Dk89BVFR5rpFC1j6wnaqP3lWF7+RI/EfP54OJUrYW6wUGpcJJzrnRETEvfz8s+kkvGmTuX56hMWkOvPw6vgfsyK2UiVYuNA0zZFiRdM6IiJS6N55B554AjIyoGJFeG/2Ee5eNggWLTIP3HUXvPsuVK5sb6GSb/Ly/Z2v55yIiIj8k/R06NvXvDIyoE0b2PruRu5+pqkJJp6eMGkSrF6tYFKMaVpHREQKRWKimcb55Rdzhsnz4yzGlHkDj/YjICvLdPFbvBhuucXuUsVmmtYREZECZVmmH9+TT0JmJlSrBh+8+Rct33oEPvzQPNSxI/z3v1CunK21SsHRbh0RESkSDh82W4SXLzfX998PCwd9i//jPWDvXvD2hmnTIDzc7MwRQeFEREQKyHffQbdukJQEJUrAyy86GXL0RRwdxoHTCXXqmE7CTZrYXaoUMVoQKyIi+So7G158EUJCTDC59lpYvzKZoavuwfHsWBNMevUynYQVTOQCXGbkRAtiRUSKvj/+gD594PPPzXW3bjD3wdWU6tsHUlLAz8+cuNa3r6Zx5KK0IFZERPLFF1+YAZEDB6BkSXhj+gke2f0sjqmnuvg1bGimcerWtbdQsYXOORERkUJz8iSMHWsOcj1wAOrXh8QP99D/ndvOBJPHH4fvv1cwkUviMtM6IiJS9OzdazoJr11rrgcOhDfaxODz0Fld/ObOhS5dbK1TXItGTkRE5LKsXAmNG5tgUro0vP/OceZ4D8anxwMmmLRoAQkJCiaSZy4TTqKioqhXrx7Nmze3uxQRkWItMxOGDoUOHcw5Js2awZZlP/PgtJvPtBceMQLi4qB2bXuLFZekBbEiInLJduwwO3A2bjTXw4fD1Abv4DXkVBe/ChVMV7+wMHsLlSJHC2JFRCTfLVoETZuaYFK+PKx6P53ph/ri1f+sLn6bNimYyBVTOBERkX+UkQGPPAI9e5quwrfeClsXbSLs2WZmlMTDAyZMgNhYqFrV7nLFDWi3joiIXNTmzaaT8Pbt5sy0Z8daPFdxFp7tI8508Vu82BwHK5JPNHIiIiLnsSx480246SYTTKpUgW9WHGb81i54Dgk3weT++800joKJ5DONnIiISC5//w2PPgoffGCu770X3hv8HeXCu8OePaaL39SpZsuOjqCXAqBwIiIiOX74wUzj/PYbeHnBlEnZRDin4tHhVMO+a6+FJUvM/mGRAuIy0zo650REpOBkZ8Mrr8Att5hgUqsW/O/DP3jqizA8Ro8yweT0HmIFEylgOudERKQYcTrN2WjJyWYdSUgI/PWXaRL86afmmS5dYF6vL7hq0Fld/F5/Hfr31zSOXLa8fH9rWkdEpJiIiTHLRPbtO3OvYkXTuO/wYfDxgdemneTR35/H0WmyWRVbv77pJFy/vn2FS7GjcCIiUgzExJgRkXPHyv/80/yzWjWIfXsvdSec08Vvxgzw8yvUWkUUTkRE3JzTaUZM/mkSP/T4Sm7s1s8MoZQuDW+9ZVbGithA4URExM3FxeWeyjmbN5lM5WmGHnrd3GjWzOzGufbawitQ5Bwus1tHREQuT3Lyhe9fxw7W0YqhmGDyU9sIM6WjYCI2UzgREXFzF5rO6c4iNtKUYDZykPLcx8ccGDENvL0Lv0CRcyiciIi4sehoGDTozLUfGcylP4voSWnS+YZbaUIimwPv0yn0UmTYEk68vLwICgoiKCiIAQMG2FGCiIhbO3rUHEHfrRscOQI33ggN2cx6mtOft8nGwfOM406+ZL+jOjNmgKen3VWLGLYsiC1btiyJiYl2/GgREbe3davZaLN1qzkzbdRIiwnV58DwYXhmHed3qtCDRXzD7QQGmt3CnTvbXbXIGdqtIyLiJiwL/vtfGDIEjh2DgABY/GYqbd4bCFNMFz/r3jB2P7aAx45V5PlTJ8RqxESKmjxP66xZs4Z27dpRtWpVHA4HK1asOO+ZmTNnUrt2bXx9fQkODiYuLi7X+2lpaQQHB9O6dWu++eabyy5eRESMtDTo0cOcm3bsGNx9N2yd9wNthjcx7YW9vODll3F88jG3dKxI9+5w++0KJlI05TmcZGRk0LhxYyIjIy/4fnR0NMOGDWPMmDEkJCQQEhJCWFgYSUlJOc/89ttvxMfHM3v2bPr06UNaWtpFf15mZiZpaWm5XiIicsaGDdC0qTmexNPTdBL+LHQa5dvfArt3my5+334LTz0FHtoHIS7AugKAtXz58lz3brrpJmvQoEG57t14443WyJEjL/gZ9957r7V+/fqL/oxx48ZZwHmv1NTUKyldRMTlZWdb1quvWlaJEpYFllWjhmV9/3GKZYWFmRtgWQ8+aFmHD9tdqoiVmpp6yd/f+Rqhs7KyiI+PJzQ0NNf90NBQ1q1bB8Dhw4fJzMwEYN++fWzbto1rrrnmop85atQoUlNTc1579+7Nz5JFRFzSoUPQvj0MHw4nTkDHjrD5ja+56dEg017Y1xdmzzZ7icuWtbdYkTzK1wWxBw8exOl0EhAQkOt+QEAABw4cAGD79u089thjeHh44HA4eO211yhXrtxFP9PHxwcfHx+ioqKIiorC6XTmZ8kiIi4nLg66d4f9+82ZadNfdvLEwQk4Or5gxktuvBHefx8aNrS7VJHLUiC7dRwOR65ry7Jy7rVq1YrNmzfn+TPDw8MJDw8nLS0Nf3//fKlTRMSVOJ0weTI8/zxkZ0OdOhDzxn4aTO4Ba9aYhx5+GN54A0qVsrVWkSuRr+GkQoUKeHp65oySnJaSknLeaIqIiFy65GTo1Qu+/NJc9+4Ns9t9gl/PvmaO56qrzDROz572FiqSD/J1zYm3tzfBwcHExsbmuh8bG0urVq2u6LOjoqKoV68ezZs3v6LPERFxNZ99Bo0bm2Di5wfvzM3inYpP4vfQ/SaYNGkCGzcqmIjbyPPISXp6Ojt37sy53r17N4mJiZQrV44aNWoQERFB7969adasGS1btmTOnDkkJSUx6OzmDpdB0zoiUtycOAFjx8LUqea6USOIeWUX147pBuvXm5tDhpgHfHzsK1Qkn+U5nGzYsIE2bdrkXEdERADQt29f5s+fT9euXTl06BATJkwgOTmZBg0asGrVKmrWrHlFhWpBrIgUJ7/9ZvrifP+9uX7iCXj15mi8uzxqTly7+mqYNw86dLC1TpGC4LCsCzXTLrpOj5ykpqZSpkwZu8sREcl3y5ZB//6Qmgr+/jB/5lE6fj0M3nrLPHDLLbB4MQQG2lqnSF7k5ftbRwWKiBQRx4+bEZIuXUwwuflm2PrBNjpOvskEE4cDxoyBr79WMBG35jKN/zStIyLu7KefTCfhH3801888bTHx2rfx6vCfM138Fi6Eu+6yt1CRQqBpHRERmy1YYEZMjh6FihVh0ew07lo6yEzdAISGwjvvmIAi4qLy8v3tMiMnIiLu5sgRCA+Hd98113fcAUuejqdieFf49VfTxW/SJBgxQg37pFhROBERsUFiIjz0EOzYYXLH+OctRpd6DY92T5s9xDVqmDbDLVvaXapIoXOZKK5D2ETEHVgWREZCixYmmFSvDt9+eIix6zvg8eSpLn6dOpn0omAixZTWnIiIFJLDh80W4eXLzXW7dvDuo3H4P94D9u071cVvulmAck6PMhFXpzUnIiJFzLp1ppNwUhKUKAGvvOTkP0cm4+jwvOnid/31EB0NQUF2lypiO5eZ1hERcUXZ2fDii3DrrSaYXHstrP8omSEfh+IY95x5oHdviI9XMBE5xWVGTnTOiYi4mj/+MLnjdC/U7t3hrS6fUap3H/jzTyhVCmbOhD597C1UpIjRmhMRkQLwf/8HvXqZgFKyJETNOEG/nWNxvHxWF7/334cbbrC3UJFCojUnIiKFwOmEuDhIToYqVSAkxOzGGTcOpkwx/7lBA4iZ/ht1nj2ni9+0aeDra+8fQKSIUjgREbkMMTEwdKjZZHNalSqmUd9PP5nrRx+F129fhs+DZ3Xx++9/4YEH7ClaxEUonIiI5FFMjGnOd+6keHKyeZUsCe/MOU6XdRHQY5Z58+abzXH0tWoVer0irkbhREQkD5xOM2LyT6v1ml31Ew+8fHYXv2fghRfMHmIR+Vcus5VYJ8SKSFEQF5d7KudcfVjAp38G4/jxR9PF77PPzF5iBRORS+Yy4SQ8PJxt27axfv16u0sRkWIsOfnC96/iCAvowwL6UYqjHKh/J2zaBPfcU7gFirgBlwknIiJFQdmy599rTCLxBNOHd3HiwRgm8tNrq80KWRHJM4UTEZFL9OOPMHz42Xcswonke1pwPTvYS3Vu5xveDRxDyO2edpUp4vIUTkRE/oVlwezZcNNN8PPPcPXVUJbDLOMBIvkPPmTxIe1pQiJrHa2ZMQM8lU1ELpvCiYjIP/j7b3joIXj8ccjMhLZt4dd317G/QhCdWU4WJRjKDDqyAr/A8ixdCp072121iGvTVmIRkYv4/nvo1g1++81stnlxcjbDTkzFo8NYcDqxrruOH59aws1lgul06oRYjZiIXDmFExGRc2Rnm9PlR4+Gkyehdm1YNvMPmkzP3cXPMXs2zcqUoZm95Yq4HZcJJ+pKLCKF4c8/oW9f+PRTc/3QQ/B291hK9et9potfZCQ8/DA4HPYWK+Km1JVYROSUr76Cnj3NWSa+vvD69JMMSHoOx0svnuniFx0N9erZXaqIy8nL97cWxIpIsXfypOkkfOedJpjUrQsJHyYxcOFtOF481V740Ufhhx8UTEQKgctM64iIFIR9+8xoyZo15rp/f4i8awW+3R6Bw4ehTBl46y0zvyMihUIjJyJSbH38MQQFmWBy1VWwZEEmc/2G4Nu9kwkmzZtDQoKCiUghUzgRkWInKwsiIqBdOzh0CJo2hS0xv9B1Rkt44w3z0FNPwbffwjXX2FusSDGkaR0RKVZ+/RW6doX4eHM9dCi83HghJTo/DunpUKECLFhgTlsTEVto5EREio0lS6BJExNMypWDj6MzmJH6MCUe6W2Cye23Q2KigomIzWwLJ0ePHqVmzZo89dRTdpUgIsXE0aMwcCB07w5HjkDr1rB18Y/cN64ZzJ8PHh7w/PPwf/8H1arZXa5IsWfbtM6kSZNo0aKFXT9eRIqJrVvNNM7WrebMtDGjLZ6v8iae7YeZZjlVq8KiRXDbbXaXKiKn2DJysmPHDn766SfaauhURAqIZZkdwM2bm2BSuTJ8veJvXvj5ITwHn9XFLzFRwUSkiMlzOFmzZg3t2rWjatWqOBwOVqxYcd4zM2fOpHbt2vj6+hIcHExcXFyu95966immTJly2UWLiPyTtDQzhfPoo3DsGNxzD2x9+3tuHdoEli41XfymTYOPPoKKFe0uV0TOkedwkpGRQePGjYmMjLzg+9HR0QwbNowxY8aQkJBASEgIYWFhJCUlAfDhhx9y/fXXc/31119Z5SIiF7Bhg1n0Gh0NXl7w0pRsVt3xMuXatzbthWvXhrVrzV5iD+0JECmKrqi3jsPhYPny5XTs2DHnXosWLWjatCmzZs3KuVe3bl06duzIlClTGDVqFAsXLsTT05P09HROnDjBk08+yXPPPXfBn5GZmUlmZmbOdVpaGoGBgeqtIyK5WBbMmAHPPAMnTkDNmrB0ZgrN3ugLn31mHnroIZgzB/z9ba1VpDiyrbdOVlYW8fHxhIaG5rofGhrKunXrAJgyZQp79+7lt99+45VXXmHgwIEXDSann/f39895BQYG5mfJIuIGDh6E9u3NYMiJE9C5M2x+7UuaDQgywcTX14SSJUsUTERcQL6Gk4MHD+J0OgkICMh1PyAggAMHDlzWZ44aNYrU1NSc1969e/OjVBFxE2vWmCPoP/4YfHxg1hsnWVrvOUp3ust08atXD9avN3uJHQ67yxWRS1AgW4kd5/wLwLKs8+4B9OvX718/y8fHBx8fH6KiooiKisLpdOZXmSLiwpxOmDQJxo+H7Gy44QaIeX0f9Sb2gNOL8Pv3h9dfBz8/e4sVkTzJ13BSoUIFPD09zxslSUlJOW80Ja/Cw8MJDw/PmbMSkeLr99+hVy/46itz3a8fzGz7MSV79DPNckqXhjffNFt2RMTl5Ou0jre3N8HBwcTGxua6HxsbS6tWra7os6OioqhXrx7Nmze/os8REdf26admGuerr6BUKXhvXhbzyg6n5EOnuvgFB8PGjQomIi4szyMn6enp7Ny5M+d69+7dJCYmUq5cOWrUqEFERAS9e/emWbNmtGzZkjlz5pCUlMSgQYOuqFCNnIgUb1lZMGYMvPKKuQ4KgpiXf6X2yLO6+A0bBi++aBafiIjLynM42bBhA23atMm5joiIAKBv377Mnz+frl27cujQISZMmEBycjINGjRg1apV1KxZM/+qFpFiZfdu6NYNfvjBXA8eDNOaL8G786OmWU65cqZHTrt2ttYpIvnjis45KUxnL4j95ZdfdM6JiJtxOs061uRkqFIFQkLA09Mc6DpgAKSmQtmysGDWUdp/MRTmzjX/h61bm944OmZApEjLyzknLhNOTsvLH05EXENMDAwdCvv2nblXrRrUrw+ff26uW7WCD57fStXhZ3fxGwPjxpmjYEWkSMvL97f+P1pEbBUTA126mBNez7Z/v3kBjBppMaHmXLw6DDXNcipXhvfegzvuKPyCRaTAuUxjCe3WEXE/TqcZMfmn8dva5dOYuKs7Xo+f1cVv0yYFExE3pmkdEbHN11/DWevrzxPMBqLpyrXsMlM3kyfDk0+qYZ+IC9K0joi4hOTki71jMYwZvMQzeHOC9Iq1uGrlYrj55sIsT0Rs4jJ//dC0joj7qVz5/HvlOchK2vMqEXhzgqU8QOLbCQomIsWIpnVExBZ//QUPPwwrV565F8IaFtGD6uznOD5E8CofVx/E7t8ceHraV6uIXDlN64hIkbZ2rTldfu9es5Qk+6STsUziOcbjSTY/cQPdiOZHR2OWvoaCiUgx4zLTOiLi+k53Er7tNhNM6tSBjZ/8zoH6dzGecXiSzTz6EUw8fwU2ZulS6NzZ7qpFpLBp5ERECsWBA6aT8BdfmOteveDNjp/i17MPHDyIVaoUPw2djW+DXnxy1gmxIlL8uEw4Ofv4ehFxLZ9/Dr17Q0oK+PnBrNey6P3TGBxdznTxc0RHU/f666lrb6kiUgRoQayIFJgTJ+C550yjYICGDSFm2m6uG3tOF7+XXwZfX/sKFZECpwWxImK7PXvMotfvvjPXjz8Or7Zeis+DZ3Xxe/tt6NTJ1jpFpOjRglgRyXfLl0NQkAkm/v6wbOExZlqP49PzQRNMWrWCxEQFExG5IIUTEck3x4+bWZrOneHvv6FFC9jywXY6v9QCZs82nYRHjTLn1tesaXe5IlJEaVpHRPLFL79A165mQARgxFMWk2+Yj1fHwXD0KFSqBAsXwt1321qniBR9LjNyouPrRYqud9+Fpk1NMKlQAVYvPcLU5N54DXzEBJO77jKdhBVMROQSaLeOiFy29HQzjbNggblu0waWjEyg0uCusGOHOahkwgQYOVKdhEWKOe3WEZECt2mTmcb5+WeTO54fZzHGPxKPdk9BVhYEBsLixXDLLXaXKiIuRn+VEZE8sSyYOdMsdv35Z6hWDb5d+RfPbuyEx7AhJph07GjmeBRMROQyaORERC7Z4cMwYADExJjr+++HhYO+xf/xHqZZjrc3TJsG4eFmZ46IyGVQOBGRS/Ldd+ZQtT17oEQJePlFJ0OOvoijwzjT0a9OHYiOhiZN7C5VRFycpnVE5B9lZ8NLL5lGfHv2wLXXwg8rDzB01T04nh1rgkmvXhAfr2AiIvlCIyciclEpKdCnD6xeba67dYO5D31Oqb5ndfGLioK+fTWNIyL5xmVGTnTOiUjh+uILaNzYBJOSJWHurBMsqjmKUp3vMcGkYUPYsAH69VMwEZF8pXNORCSXkyfh+edh8mSzM6d+fYh5dQ/Xjzuni9+0aSa1iIhcAp1zIiKXZe9e6NEDvv3WXA8cCG/csRyfhx4xzXL8/WHuXOjSxdY6RcS9ucy0jogUrJUrzTTOt99C6dLw/jvHmeM9GJ/uZ3XxS0hQMBGRAqdwIlLMZWbC0KHQoYM5x6RZM9gS8wsPTm9pFrsCjBgBcXFQu7a9xYpIsaBpHZFibMcOswNn40ZzPXw4TG34Ll4dH4eMDNPF7513ICzM3kJFpFjRyIlIMbVokekkvHEjlC8Pq95PZ/pf/fB6pI8JJm3amAY6CiYiUsgUTkSKmYwM6N8fevY0XYVvvRW2LtpE2LPNTHthDw/TSTg2FqpWtbtcESmGCn1a58iRI9xxxx2cOHECp9PJkCFDGDhwYGGXIeLWnE6zRCQ5GapUMae7enrC5s2mk/D27eZokueetXiu4iw82keYxSfVqpkhlVtvtfuPICLFWKGHEz8/P7755hv8/Pw4evQoDRo0oHPnzpQvX76wSxFxSzExZoHrvn1n7lWvDvfdZwZGjh83geX9Nw/TesFAWLbMPHT//TBvnllnIiJio0IPJ56envj5+QFw/PhxnE4nLnYOnEiRFRNjdvqe+/9S+/bBm2+a/xwWBu8N/o6rnziri99LL8GwYTrpVUSKhDyvOVmzZg3t2rWjatWqOBwOVqxYcd4zM2fOpHbt2vj6+hIcHExcXFyu9//++28aN25M9erVefrpp6mgv6mJXDGn04yY/FPWv9o/m49av8TV7U918bvmGli3zmzTUTARkSIiz+EkIyODxo0bExkZecH3o6OjGTZsGGPGjCEhIYGQkBDCwsJISkrKeaZs2bJs2rSJ3bt3s2jRIv7444+L/rzMzEzS0tJyvUTkfHFxuadyzlWRFBaltsVzzEiTZLp2NVt1mjUrvCJFRC5BnsNJWFgYEydOpHPnzhd8f/r06fTv358BAwZQt25dZsyYQWBgILNmzTrv2YCAABo1asSaNWsu+vOmTJmCv79/ziswMDCvJYsUC8nJF3/vDr5gE425l9Wc9C4Jb70Fixeb4+hFRIqYfN1KnJWVRXx8PKGhobnuh4aGsm7dOgD++OOPnNGPtLQ01qxZww033HDRzxw1ahSpqak5r7179+ZnySJuo0qV8+95cpIXGEssd1OFA2ylHhtnr4cBAzSNIyJFVr4uiD148CBOp5OAgIBc9wMCAjhw4AAA+/bto3///liWhWVZDB48mEaNGl30M318fPDx8SEqKoqoqCicTmd+liziFpxO+OKL3Peqs5dF9CAE08VvLgOYWu01tvfxs6FCEZFLVyC7dRzn/I3Msqyce8HBwSQmJub5M8PDwwkPD89puSwixv79ppPw2bOj7VnJ2zxMef4ijdI8xhyiHd1Y+ro570REpCjL12mdChUq4OnpmTNKclpKSsp5oykicuU++cR0El6zBq66ChbNy2THfcP4kA6U5y82EExTNrI2sBtLl8JFloqJiBQp+RpOvL29CQ4OJjY2Ntf92NhYWrVqdUWfHRUVRb169WjevPkVfY6IO8jKgiefNOemHToETZrA5uU76f5GK6775DUA9nYZzq/vrGPuV9exe7eCiYi4jjxP66Snp7Nz586c6927d5OYmEi5cuWoUaMGERER9O7dm2bNmtGyZUvmzJlDUlISgwYNuqJCNa0jYuzaZToJr19vrocMgVeCF1Oi82Nw5Ijp4jd/PoH3309Xe0sVEbkseQ4nGzZsoE2bNjnXERERAPTt25f58+fTtWtXDh06xIQJE0hOTqZBgwasWrWKmjVrXlGhWhArAtHR8OijkJYGV18N78zK4P7Ph0Dft80Dt94K771nzqsXEXFRDsvFzo4/PXKSmppKmTJl7C5HpFAcPWpOl3/rLXN9yy2w9PktVB7y0Jkufs8+a15ehd6VQkTkX+Xl+1v/FhMp4rZtg4cegq1bTQYZPcpifPW38Gw39EwXv/feg7NGNEVEXFm+LogtSFoQK8WNZcF//2tOl9+6FQIC4IvlqUzc0RXPJx4zweTeeyExUcFERNyKpnVEiqC0NBg0yJwwDxAaCouG/UD58G6we7eZupkyBSIiwMNl/o4hIsWYpnVEXNiGDWY3zq+/mgPTJr2QzYgSr+LRfiScPAm1asGSJdCihd2liogUCP2VS6SIsCyYMQNatTLBpEYN+G7lnzzzbTs8RjxlgkmXLpCQoGAiIm7NZUZOtJVY3NmhQ9CvH3z8sbnu1AkWPPw1pQf2hN9/Bx8fk1wee0wN+0TE7WnNiYjN4uKge3fTI8fbG159xcnjB1/AMfEFyM6GG280B5z8Q4NMEZGiLi/f35rWEbGJ0wkvvAC3326CyfXXw8aP9vPE0jtwTBhvgkm/fmYRioKJiBQjLjOtI+JOkpOhVy/48ktz3bs3vNn+E0r26GvmeK66CmbNMg+JiBQzLjNyonNOxF189pnpJPzll1CqFLz73yzeqfgkJR88q4tffLyCiYgUW1pzIlJITpyAsWNh6lRz3agRLJ+2i2tGn9PFb+pUswBWRMSN6JwTkSLmt9/M2SXff2+uw8Nh+s3v4/3AwDNd/ObNgw4dbK1TRKQocJlpHRFXtWwZBAWZYFK2LCxfdIzIE4/h3burCSa33GKOoFcwEREBFE5ECszx4/DEE+bctNRUuPlm2PL+NjpOvgnmzDnVxW80fP21OXFNREQATeuIFIiffoKuXeHHH831M09bTLxuHl4dBsOxY6aL38KFcNdd9hYqIlIEuczIiXbriKtYsACCg00wqVgRYpel8eLenng92t8Ek7vvhk2bFExERC5Cu3VE8smRI2ah67vvmus774TFI+KpGN71TBe/iRPh6afVSVhEih3t1hEpZAkJZhpnxw6TOyaMtxjl9xoe7Z42e4hr1IDFi01XPxER+Uf665vIFbAseOMNs9h1xw6oXh3WrjzEmB864PHkcBNMOnY06UXBRETkkmjkROQy/fUX9O8PK1aY6/bt4Z2Bcfg/dlYXv2nTzFyPOgmLiFwyhRORf+F0ms7ByclQpQqEhJgzS7p3h6QkKFECXnnJyX+OTMbR4XnTsK9OHdNJuEkTu8sXEXE5Cici/yAmBoYOhX37ztwrUwbS000Gue46WBqZTOOp53Txi4qC0qXtKVpExMUpnIhcREyMOUDt3P1saWnmnyEh8Omw1ZTq3Rv+/BP8/GDmTOjbt/CLFRFxIy6zIFbnnEhhcjrNiMnFNtp7cYLuic9Q6oF7TTBp1Mh0ElYwERG5YjrnROQCvv4a2rS58Hs1+Y0ldONmTnXxe+IJs/DV17fQ6hMRcTV5+f52mZETkcKUnHzh+51ZRiJB3Mz3/I0/ccOWmvUlCiYiIvlG4UTkAn79Nfe1D8eJ4gmW0YWypPIdNxNEIs4OD9hToIiIG1M4ETlLZiYMGQLPPnvm3g38xPe04AlmAfAiz3Aba8gOrEVIiE2Fioi4Me3WETnll1+gWzdzmCuYQ9WuXrmAKJ6gFEdJoSJ9eIfPHfcCMGOGaZcjIiL5SyMnIsDChaaTcEIClC8Pn75/hA/9+zCffpTiKF9wB43ZxGrupXp1WLoUOne2u2oREfekkRMp1jIyYPBgmD/fXN92G7w/OpFKgx/K6eKXPW48nreMYnqKZ84JsRoxEREpOIUeTvbu3Uvv3r1JSUnBy8uLZ599lgcffLCwyxDhxx9NJ+GffjKdhJ8da/Fc+Sg82j0JWVmmi9+iRXiEhHC73cWKiBQjhR5OvLy8mDFjBkFBQaSkpNC0aVPatm1LqVKlCrsUKaYsC958E4YNMwtgq1aF6NmHaT2vPyxfbh5q1w7mzTNzPCIiUqgKPZxUqVKFKlWqAFCpUiXKlSvHX3/9pXAiheLvv2HgQLNmBKBtW1j4xDqufuKsLn4vv2y27KiTsIiILfK8IHbNmjW0a9eOqlWr4nA4WHG6X/xZZs6cSe3atfH19SU4OJi4uLgLftaGDRvIzs4mMDAwz4WL5NX335smwUuXgpcXvDI1m49umcLVHW41weTaa+G778y59QomIiK2yXM4ycjIoHHjxkRGRl7w/ejoaIYNG8aYMWNISEggJCSEsLAwkpKScj136NAh+vTpw5w5cy6vcpFLlJ1tBkNat4bffoPateH7lX/wZOy9eIwZbRrpdO8OGzeaLTsiImKrK+qt43A4WL58OR07dsy516JFC5o2bcqsWbNy7tWtW5eOHTsyZcoUADIzM7n77rsZOHAgvXv3/sefkZmZSWZmZs51WloagYGB6q0jl+TPP00vvk8/NdcPPgjzesRSalBv+OMPKFkSIiPh4Yc1WiIiUoBs662TlZVFfHw8oaGhue6Hhoaybt06ACzLol+/ftxxxx3/GkwApkyZgr+/f85LU0Byqb7+GoKCTDDx9YU5M08Sfe1oSnW+xwSTBg1gwwZ45BEFExGRIiRfw8nBgwdxOp0EBATkuh8QEMCBAwcAWLt2LdHR0axYsYKgoCCCgoLYvHnzRT9z1KhRpKam5rz27t2bnyWLG3I6Ydw4uOMO+P13qFsXEj5MYuDC23C8OMVs13n0UfjhB6hXz+5yRUTkHAWyW8dxzt9CLcvKude6dWuys7Mv+bN8fHzw8fEhKiqKqKgonE5nvtYq7mX/fujRA9asMdePPAJRd6/At9sjcPgwlCkDb70FDz1kb6EiInJR+TpyUqFCBTw9PXNGSU5LSUk5bzQlr8LDw9m2bRvr16+/os8R9/XJJ9C4sQkmV10Fi+dn8t9SQ/Dt3skEk+bNzfn0CiYiIkVavoYTb29vgoODiY2NzXU/NjaWVq1aXdFnR0VFUa9ePZo3b35FnyPuJysLnnwS7r8fDh2Cpk1h87Jf6PZaS3jjDfPQk0/Ct9/CNdfYW6yIiPyrPE/rpKens3Pnzpzr3bt3k5iYSLly5ahRowYRERH07t2bZs2a0bJlS+bMmUNSUhKDBg26okLDw8MJDw/PWe0rArBrl+kkfHpAbcgQeCVoISUeeBzS080JrwsWwH332VuoiIhcsjyHkw0bNtCmTZuc64iICAD69u3L/Pnz6dq1K4cOHWLChAkkJyfToEEDVq1aRc2aNfOvahEgOtqsa01Lg6uvhndnZ3Dfp4Phkfnmgdtug/feg2rVbK1TRETy5orOOSlMZy+I/eWXX3TOSTF29Kjpi/PWW+b6lltg6XM/Unno2V38njUvtQ8WESkS8nLOicuEk9Py8ocT97Ntm1nPunWrOZpk9CiL8VXfxPPJYWe6+L33Htx+u92liojIWWw7hE2koFgW/Pe/0KyZCSYBAfBlzN9M/OUhPAc/boJJ27aQmKhgIiLi4lwmnGi3TvGVlgY9e8KAAXDsGNx9N2yd9z23Dz+7i98r8NFHULGi3eWKiMgV0rSOFGnx8dC1K/z6q1k+MnFCNk97TTMN+06eNF38liyBm26yu1QREfkHmtYRl2dZMGMGtGxpgkmNGrDuwz8Z+e39eDzztAkmDz5oDlVTMBERcSsFcny9yJU4dMg0Cf7oI3PdsSMs6PcVZQb2hORk08Xvtddg4EA17BMRcUMuE07UW6d4iIszvXH27QNvb5g+9SRPHHoBR6cXzHBK3brmgJOGDe0uVURECojWnEiR4HTClCmmm3B2NtSpAzGv76PBlJ5nuvj1729GTEqVsrdYERHJs7x8f7vMyIm4B6fTjI4kJ0OVKhASAikp0KsXfPmleaZ3b5h9/8f49epn5niuugrmzIHu3W2tXURECofCiRSamBgYOtRM2ZxWoYJp3JeWBn5+MPv1LHpveQa6zjAPNG1qpnGuu86WmkVEpPC5TDjRmhPXFhMDXbqYZSNnO3jQ/LNmTfhizq9cO7qr2T8MJsm89BL4+BRusSIiYiutOZEC53RCrVq5R0zO9US5JUSeeBTHkSNQrhzMmwft2xdajSIiUrC05kSKlLi4iweTkhzlNYYy8K+55kbr1rBoEQQGFl6BIiJSpOgQNilwyckXvl+PrfzATQxkLtk42NJxLHz1lYKJiEgxp3AiBS4r69w7FgN4i/U0pwFbSaYydxPLwaEvmD45IiJSrCmcSIE53Ul40KAz90qTxmK68xaP4scxVhNKExLZEXgnISH21SoiIkWHy4QTdSV2Lamp5qTXAQPg+HFo1AiasYEEmtCNaE7iyTO8SFs+JcURwIwZprGfiIiIdutIvlu/Hrp1g127TOCYNNFihPcMeOYZPE6eYA816MYS/kdLAgNNg7/One2uWkRECpJ264gtsrNh+nQYNco0Da5ZEz6YdZDmMx+Gjz8GwOrUmb0Pz2VI+tVMOXVCrEZMRETkbAonki9SUqBvX/jsM3P9wAMw7+E1lB7YA/bvNwepvfoqjkGDaK1OwiIi8g8UTuSKffml6Y2TnAy+vvDadCcDUybhaD/eDKfccIM5gr5xY7tLFRERF+AyC2Kl6Dl5EsaOhbvuMsGkbl3Y+PHvPPrB3TieP9VeuG9f2LBBwURERC6ZRk7ksiQlmd04a9ea64ED4fWwT/Ht3hf+/BNKlYJZs0yLYRERkTzQyInk2fLlZiBk7VooUwbeX5jFHP8R+HZua4JJ48ameZ+CiYiIXAaXCSc658R+x49DeLjZ9vv333DTTbB55W4efD0EXnnFPDR4MPzvf2adiYiIyGXQOSdySbZvN2eX/PijuR4xAiY1XUqJQQPMiWtly8Lbb0OnTrbWKSIiRVNevr9dZuRE7GFZMG8eNGtmgknFivD5h8eYeuRxSnR/0ASTli0hMVHBRERE8oXCiVxUWhr07AmPPAJHj8Kdd8LWpdu5e2wLmD3bPDRyJHzzjTlxTUREJB9ot45c0IYNZhrn11/NCa4vTLB4JmA+HmGDTVKpVAnefRdCQ+0uVURE3IxGTiSX00fQt2plgkmNGvDtp0cYta03HgNODaHcdRds2qRgIiIiBUIjJ5Ljzz+hXz9Ytcpcd+4M84YkUGZgV9ixwwyhTJhgpnI8lGtFRKRgKJwIAF99ZdaXJCefaoMz3WLQyUgcoU9BVhYEBsLixXDLLXaXKiIibs6Wv/526tSJq6++mi5dutjx4+UsJ0/Cc8+Zxa6nj6CPj/2Lxz/vhGPoEBNMOnQwu3EUTEREpBDYEk6GDBnCO++8Y8ePlrMkJUGbNvDCC2bLcP/+EP/6Wur3DIIPPwRvb3j9dXMkbLlydpcrIiLFhC3hpE2bNpQuXdqOHy2nrFgBQUHw7bdQujQsec/J3GsmU/Le22DvXrjuOvjuO/jPf8DhsLtcEREpRvIcTtasWUO7du2oWrUqDoeDFStWnPfMzJkzqV27Nr6+vgQHBxMXF5cftUo+OH7c5I1OneDwYWjeHH78/ABd374HxowBp9MsPtm4EZo2tbtcEREphvIcTjIyMmjcuDGRkZEXfD86Opphw4YxZswYEhISCAkJISwsjKSkpMsqMDMzk7S0tFwvuTw//QQ33wyn/6t76ilYO+5zanVoDF98AX5+5gj6d981wykiIiI2yHM4CQsLY+LEiXTu3PmC70+fPp3+/fszYMAA6taty4wZMwgMDGTWrFmXVeCUKVPw9/fPeQUGBl7W5xRnlgXz50NwsDmepGJF+HTlCV72GkWJ+++BlBRo2NCcvPbww5rGERERW+XrmpOsrCzi4+MJPedwrtDQUNatW3dZnzlq1ChSU1NzXnv37s2PUouNI0egd2+TOY4ehTvugM0f7+HeKbfBiy+ahwYNgu+/N1t1REREbJav55wcPHgQp9NJQEBArvsBAQEcOHAg5/qee+5h48aNZGRkUL16dZYvX07z5s0v+Jk+Pj74+PgQFRVFVFQUTqczP0t2a/Hx5gj6nTvN+Wnjx8PIG5fjec8j8Pff4O8Pb70FDz5od6kiIiI5CuQQNsc50wKWZeW6t3r16jx/Znh4OOHh4Tktl+XiLAteew2efhpOnDDnpy2Zf5xWMU/B2CjzUIsW5lC12rXtLVZEROQc+RpOKlSogKenZ65REoCUlJTzRlPkyjmdEBdnDk+rUgVCQswOnIcfho8/Ns906gRvj/yFso91NQepAYwYAZMmQYkSttUuIiJyMfkaTry9vQkODiY2NpZOnTrl3I+NjaVDhw5X9Nma1sktJgaGDoV9+87cq1jRnPh6+LA5gn76dHj8qndx3PE4ZGRAhQrwzjsQFmZf4SIiIv8iz+EkPT2dnTt35lzv3r2bxMREypUrR40aNYiIiKB37940a9aMli1bMmfOHJKSkhg0aNAVFappnTNiYqBLFzN9c7Y//zT/rFoVPluaTsM3B8OCBeZmmzawcKF5U0REpAjLczjZsGEDbdq0ybmOiIgAoG/fvsyfP5+uXbty6NAhJkyYQHJyMg0aNGDVqlXUrFnzigrVyInhdJoRk3ODydnqn9xEg4e7ws8/m+7Bzz8Po0ebVbEiIiJFnMOy/ulrrug5PXKSmppKmTJl7C6n0H39tRkEuTCLx5nFdCLwJROqVYNFi+DWWwuxQhERkfPl5fu7QHbrSMFJTr7wfX/+Zi4D6MIyAPY3uZ9qn88z60xERERciC2N/y5HVFQU9erVu+h5KMXF4cPn32vB/0gkiC4sI4sSDGc6O6atVDARERGXpGkdF3HyJEycCC+8ANnZ5p6DbEbwMpMYgxdOfuUauhHNH4HN2L1bS0xERKTo0LSOm9m5E3r1MifMgznP5Oe4FBbQh3sxB9otoSuDeJM0hz9LZyiYiIiI63KZaZ3iyLJMk+CgIBNM/P3hvfdgzbgv2FO2MfeymqOUZABv0Z3FlAn0Z+lSuEhPRhEREZfgMiMnxW0r8aFD8Oij5kwTMBtu3nn7JDXnPQ+9JuNrWVj167PlyWju9K1Pr1MnxGrEREREXJ3WnBRBsbHQt6/ZmePlZdaaPNV1L569e8C335qHBg6EGTPAz8/WWkVERC6F1py4qOPHYdQokzkAbrjBHFPSdN9KCH4Y/voLSpc2nYS7drW1VhERkYKiNSdFxObN0Lz5mWDyxBOw8btMmr4zDDp0MMGkWTNISFAwERERt+Yy4cRdzznJzoZXXzW5Y8sWqFTJdBSOGr4Tv7tawWuvmQeHD4e1a+Haa+0tWEREpIBpzYmN9u+Hfv3g//7PXN9/P/z3v1Dpi8Xw2GNw5AiULw/z55s3RUREXFRevr9dZuTE3SxbBo0amWBSsiTMmgUrF2dQaVR/6NHDBJNbb4XERAUTEREpVrQgtpAdOQJDhpjBEICmTc3ZJTee3AI3PQTbt4PDAc8+a15e+q9IRESKF42cFKLvvjMHqs2fb/LHqFHw3TqLG9fMMatht2+HKlXgiy9g/HgFExERKZb07VcITvfFmTgRnE6oWRPefRdCGqVCn0fh/ffNg/feCwsWmFWxIiIixZTLjJy46m6dnTuhdWszEOJ0mh45mzZBiO96M6fz/vtmhOTll+GTTxRMRESk2NNunQJyui/O0KGQkWH64syeDd0eOrV3eORIM6RSqxYsWQItWthdsoiISIHRCbE2O3TInC6/fLm5vu02eOcdqOF3ENr1hVWrzBtdupjTXsuWta1WERGRosZlpnVcxeefQ8OGJpiUKAEvvWTWt9bY/Q00bmyCia+vGUZ5/30FExERkXNo5CSfHD9uZmpOH+h6441mi3DTxk544QXzys42b7z/vkkwIiIich6Fk3zw44/m3LStW811eDhMnQp+h/fDnT3hm2/MGw8/DG+8AaVK2VesiIhIEadpnSuQnQ3Tp5sjSrZuNRttPvkEIiPB7+tV5lCTb76Bq66ChQvNClkFExERkX+kkZPLtH8/9O1r1pMAtGsHc+dCpbJZ8NRomDbNvNGkCURHQ5069hUrIiLiQlxm5KQonXOydKlZMvLFF6YvzuzZ8OGHUCl9lznU5HQwGTLEHAurYCIiInLJdM5Jnn62ObfkdF+cZs3MbM0NN2AWuQ4caB66+mqYNw86dCjU+kRERIoqdSUuAOvWnemL4+EBY8aYezfUOAaPPQZdu5pgcsstppOwgomIiMhlUTj5FydOwHPPQUgI7N5t+uJ8/bXpk1Nixza46SaYM8d08hs92rxZo4bdZYuIiLgsLYg9xemEuDhITjaNgU+HkV694PvvzTO9epmdOP5lLHh7HgweDMeOQUCAmd+56y57/xAiIiJuQOEEiIkxa0n27Ttz7+qr4ehRyMw0h7jOmgXdumGmbno9DosWmQfvvtu0GA4IsKN0ERERt1Psw0lMjGlxc+6y4MOHzT/r14dPP4XAQCA+3iSUnTvB09PM7Tz9tFmEIiIiIvmiWIcTp9OMmPzTfqW0NKhaxYLXXocRI8wilBo1YPFiaNWq8IoVEREpJop1OImLyz2VcyEZew9x+LZHqLBupbnRqZM5ba1cuYIvUEREpBiyZT7i448/5oYbbqBOnTrMnTvXjhIAs/j1n9zCtyQSZIKJt7dZDbtsmYKJiIhIASr0kZOTJ08SERHBV199RZkyZWjatCmdO3emnA1f+FWqXPi+B05G8iLjGYcXTo4GXo/fymhz0ImIiIgUqEIfOfnhhx+oX78+1apVo3Tp0rRt25bVq1cXdhmA2S5cvbo5ouS0yiSzmnuYxFi8cLLMrzc+m+MVTERERApJnsPJmjVraNeuHVWrVsXhcLBixYrznpk5cya1a9fG19eX4OBg4uLict77/fffqVatWs519erV2b9//+VVf4U8PeG118x/djgglNUkEsRdfEE6pejLAhzvvoOn/1W21CciIlIc5TmcZGRk0LhxYyIjIy/4fnR0NMOGDWPMmDEkJCQQEhJCWFgYSUlJAFyolY/j7KGLc2RmZpKWlpbrlZ86dzaN/CaUeZnV3EsAKWyiEfdXjqfDsj507pyvP05ERET+RZ7DSVhYGBMnTqTzRb61p0+fTv/+/RkwYAB169ZlxowZBAYGMmvWLACqVauWa6Rk3759VLnY4g9gypQp+Pv757wCAwPzWvK/6twZRi1vgeXhwS93P0Hq6u/5Yt8NCiYiIiI2yNc1J1lZWcTHxxMaGprrfmhoKOvWrQPgpptuYsuWLezfv58jR46watUq7rnnnot+5qhRo0hNTc157d27Nz9LzuHZ5lYc27Zx/edR3Brqi6dngfwYERER+Rf5ulvn4MGDOJ1OAs45yj0gIIADBw6YH+jlxbRp02jTpg3Z2dk8/fTTlC9f/qKf6ePjg4+PD1FRUURFReF0OvOz5NxuuKHgPltEREQuSYFsJT53DYllWbnutW/fnvbt2+fpM8PDwwkPDyctLQ1/f/98qVNERESKnnyd1qlQoQKenp45oySnpaSknDeakldRUVHUq1eP5s2bX9HniIiISNGWr+HE29ub4OBgYmNjc92PjY2l1RX2oQkPD2fbtm2sX7/+ij5HREREirY8T+ukp6ezc+fOnOvdu3eTmJhIuXLlqFGjBhEREfTu3ZtmzZrRsmVL5syZQ1JSEoMGDcrXwkVERMQ95TmcbNiwgTZt2uRcR0REANC3b1/mz59P165dOXToEBMmTCA5OZkGDRqwatUqataseUWFFsqCWBEREbGdw7rQqWhF2OkFsampqZQpU8buckREROQS5OX725auxCIiIiIX4zLhRLt1REREigdN64iIiEiB07SOiIiIuCyFExERESlSXCacaM2JiIhI8aA1JyIiIlLg8vL9XSCN/wrS6SyVlpZmcyUiIiJyqU5/b1/KmIjLhZMjR44AEBgYaHMlIiIikldHjhzB39//H59xuWmd7Oxsfv/9d0qXLo3D4cjXz05LSyMwMJC9e/dqyqgA6fdcOPR7Lhz6PRcO/Z4LT0H9ri3L4siRI1StWhUPj39e8upyIyceHh5Ur169QH9GmTJl9D/+QqDfc+HQ77lw6PdcOPR7LjwF8bv+txGT01xmt46IiIgUDwonIiIiUqQonJzFx8eHcePG4ePjY3cpbk2/58Kh33Ph0O+5cOj3XHiKwu/a5RbEioiIiHvTyImIiIgUKQonIiIiUqQonIiIiEiRonAiIiIiRYrCySkzZ86kdu3a+Pr6EhwcTFxcnN0luZUpU6bQvHlzSpcuTaVKlejYsSM///yz3WW5vSlTpuBwOBg2bJjdpbil/fv306tXL8qXL4+fnx9BQUHEx8fbXZZbOXnyJGPHjqV27dqULFmSa665hgkTJpCdnW13aS5tzZo1tGvXjqpVq+JwOFixYkWu9y3L4vnnn6dq1aqULFmS22+/na1btxZafQonQHR0NMOGDWPMmDEkJCQQEhJCWFgYSUlJdpfmNr755hvCw8P53//+R2xsLCdPniQ0NJSMjAy7S3Nb69evZ86cOTRq1MjuUtzS4cOHueWWWyhRogSffvop27ZtY9q0aZQtW9bu0tzKSy+9xOzZs4mMjGT79u1MnTqVl19+mTfeeMPu0lxaRkYGjRs3JjIy8oLvT506lenTpxMZGcn69eupXLkyd999d05/uwJniXXTTTdZgwYNynXvxhtvtEaOHGlTRe4vJSXFAqxvvvnG7lLc0pEjR6w6depYsbGx1m233WYNHTrU7pLczjPPPGO1bt3a7jLc3n333Wc98sgjue517tzZ6tWrl00VuR/AWr58ec51dna2VblyZevFF1/MuXf8+HHL39/fmj17dqHUVOxHTrKysoiPjyc0NDTX/dDQUNatW2dTVe4vNTUVgHLlytlciXsKDw/nvvvu46677rK7FLe1cuVKmjVrxoMPPkilSpVo0qQJb731lt1luZ3WrVvzxRdf8MsvvwCwadMmvv32W9q2bWtzZe5r9+7dHDhwINf3oo+PD7fddluhfS+6XOO//Hbw4EGcTicBAQG57gcEBHDgwAGbqnJvlmURERFB69atadCggd3luJ0lS5awceNG1q9fb3cpbm3Xrl3MmjWLiIgIRo8ezQ8//MCQIUPw8fGhT58+dpfnNp555hlSU1O58cYb8fT0xOl0MmnSJLp37253aW7r9Hffhb4X9+zZUyg1FPtwcprD4ch1bVnWefckfwwePJgff/yRb7/91u5S3M7evXsZOnQon3/+Ob6+vnaX49ays7Np1qwZkydPBqBJkyZs3bqVWbNmKZzko+joaBYuXMiiRYuoX78+iYmJDBs2jKpVq9K3b1+7y3Nrdn4vFvtwUqFCBTw9Pc8bJUlJSTkvNcqV+89//sPKlStZs2YN1atXt7sctxMfH09KSgrBwcE595xOJ2vWrCEyMpLMzEw8PT1trNB9VKlShXr16uW6V7duXZYtW2ZTRe5pxIgRjBw5km7dugHQsGFD9uzZw5QpUxROCkjlypUBM4JSpUqVnPuF+b1Y7NeceHt7ExwcTGxsbK77sbGxtGrVyqaq3I9lWQwePJiYmBi+/PJLateubXdJbunOO+9k8+bNJCYm5ryaNWtGz549SUxMVDDJR7fccst52+F/+eUXatasaVNF7uno0aN4eOT+qvL09NRW4gJUu3ZtKleunOt7MSsri2+++abQvheL/cgJQEREBL1796ZZs2a0bNmSOXPmkJSUxKBBg+wuzW2Eh4ezaNEiPvzwQ0qXLp0zUuXv70/JkiVtrs59lC5d+rx1PKVKlaJ8+fJa35PPhg8fTqtWrZg8eTIPPfQQP/zwA3PmzGHOnDl2l+ZW2rVrx6RJk6hRowb169cnISGB6dOn88gjj9hdmktLT09n586dOde7d+8mMTGRcuXKUaNGDYYNG8bkyZOpU6cOderUYfLkyfj5+dGjR4/CKbBQ9gS5gKioKKtmzZqWt7e31bRpU21xzWfABV/z5s2zuzS3p63EBeejjz6yGjRoYPn4+Fg33nijNWfOHLtLcjtpaWnW0KFDrRo1ali+vr7WNddcY40ZM8bKzMy0uzSX9tVXX13w38l9+/a1LMtsJx43bpxVuXJly8fHx7r11lutzZs3F1p9DsuyrMKJQSIiIiL/rtivOREREZGiReFEREREihSFExERESlSFE5ERESkSFE4ERERkSJF4URERESKFIUTERERKVIUTkRERKRIUTgRERGRIkXhRERERIoUhRMREREpUhRORMR2O3fuxOFw8Mknn3DnnXfi5+fHDTfcwPfff293aSJiA4UTEbHdpk2bcDgcTJs2jbFjx7Jp0yZq1KjByJEj7S5NRGygcCIittu0aRP+/v5ER0fTpk0b6tSpQ8eOHfnzzz/tLk1EbKBwIiK227RpE+3ataNixYo593bt2sV1111nY1UiYheFExGx3aZNm2jZsmWuewkJCQQFBdlTkIjYSuFERGyVmprKnj17aNKkSa77iYmJCicixZTCiYjYatOmTXh6etK4ceOce3v27OHw4cMKJyLFlMKJiNhq06ZN3HjjjZQsWTLnXkJCAmXLlqVWrVr2FSYitnFYlmXZXYSIiIjIaRo5ERERkSJF4URERESKFIUTERERKVIUTkRERKRIUTgRERGRIkXhRERERIoUhRMREREpUhROREREpEhROBEREZEiReFEREREihSFExERESlS/h9uOtibl5TREQAAAABJRU5ErkJggg==", + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "PyObject Text(0.5, 24.0, '$n$')" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using PyPlot\n", + "semilogy(0:10, [norm(A^n * x) for n=0:10], \"bo-\")\n", + "semilogy(0:10, λ[2] .^ (0:10), \"r-\")\n", + "legend([L\"\\Vert A^n x \\Vert\", L\"\\lambda_2^n\"])\n", + "xlabel(L\"n\")" + ] + }, + { + "cell_type": "markdown", + "id": "8605fdc7", + "metadata": {}, + "source": [ + "### Example: Cyclic permutations:\n", + "\n", + "Another example is the $4 \\times 4$ cyclic permutation matrix $P$" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "2b07afee", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4×4 Matrix{Int64}:\n", + " 0 1 0 0\n", + " 0 0 1 0\n", + " 0 0 0 1\n", + " 1 0 0 0" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P = [ 0 1 0 0\n", + " 0 0 1 0\n", + " 0 0 0 1\n", + " 1 0 0 0 ]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "91f2be85", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4-element Vector{Int64}:\n", + " 2\n", + " 3\n", + " 4\n", + " 1" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P * [1,2,3,4]" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "950df8e5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4-element Vector{Int64}:\n", + " 1\n", + " 2\n", + " 3\n", + " 4" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P^4 * [1,2,3,4]" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "3150598c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4×4 Matrix{Int64}:\n", + " 1 0 0 0\n", + " 0 1 0 0\n", + " 0 0 1 0\n", + " 0 0 0 1" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P^4" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "a6d10f08", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4-element Vector{ComplexF64}:\n", + " -0.9999999999999997 + 0.0im\n", + " 1.2902104469069482e-16 - 1.0000000000000002im\n", + " 1.2902104469069482e-16 + 1.0000000000000002im\n", + " 1.0 + 0.0im" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigvals(P) # just computes eigenvalues" + ] + }, + { + "cell_type": "markdown", + "id": "9ad66bd0", + "metadata": {}, + "source": [ + "Up to roundoff errors, the eigenvalues are $-1, -i, +i, +1$.\n", + "\n", + "This is a consequence of the fact that $P^4 = I$, which means that $\\lambda^4 = 1$:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "fd3d134d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4×4 Matrix{Int64}:\n", + " 1 0 0 0\n", + " 0 1 0 0\n", + " 0 0 1 0\n", + " 0 0 0 1" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P^4" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "266caef6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "4-element Vector{ComplexF64}:\n", + " 0.9999999999999987 - 0.0im\n", + " 1.0000000000000009 + 5.160841787627796e-16im\n", + " 1.0000000000000009 - 5.160841787627796e-16im\n", + " 1.0 + 0.0im" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigvals(P) .^ 4 # each element to 4th power" + ] + }, + { + "cell_type": "markdown", + "id": "8bbc42ce", + "metadata": {}, + "source": [ + "### Complex eigenvalues of real matrices\n", + "\n", + "Notice that, even though $P$ is a *real* matrix, its **eigenvalues (and eigenvectors)** can be **complex**.\n", + "\n", + "This isn't too surprising if you think about the connection to roots of polynomials — even for the quadratic equation, you know that *real polynomials can have complex roots*.\n", + "\n", + "If you make a random real matrix, it almost always has some complex roots if the matrix is big enough. For example:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "22d5edcc", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "20×20 Matrix{Float64}:\n", + " 0.829171 0.138794 0.678026 … -0.91498 0.948123 -0.803583\n", + " 2.19259 0.748492 0.414546 0.171528 0.68522 0.223575\n", + " 0.816218 0.153332 -0.914678 -1.20477 -0.153085 0.521608\n", + " -1.47155 -1.49438 -0.274829 0.218383 0.0826769 -0.469391\n", + " 0.743134 -0.177163 0.313937 -1.24328 0.613885 0.271068\n", + " -0.186778 -0.171095 -1.03668 … 2.46496 2.15091 0.0320823\n", + " 0.301177 -0.115305 -0.339294 -0.488035 0.46601 -0.365528\n", + " 1.73961 0.323871 -1.09787 -0.503415 -1.14928 1.29621\n", + " -0.792071 0.683768 -0.908675 -1.20089 -0.672091 0.44782\n", + " -0.351495 -0.536343 0.2996 -2.24491 1.40864 -1.52049\n", + " -0.370341 1.21907 -0.896898 … -0.180509 -0.429355 -0.141925\n", + " 0.989002 0.140579 0.300553 -0.914419 0.0226391 0.141268\n", + " -0.793007 -0.101244 0.21641 -1.11298 0.467109 0.575488\n", + " 0.281743 -0.593686 -0.147228 -0.0557582 -0.537067 -0.317039\n", + " -0.613932 0.285307 -1.4214 0.281757 -2.05607 -0.0201645\n", + " 0.56785 0.164016 0.0313567 … -1.22978 1.06791 -0.235113\n", + " 0.309827 -0.558453 -1.46202 -0.826678 -1.0413 -0.769869\n", + " 0.69229 1.00292 -0.879679 -0.797199 1.12466 -0.839403\n", + " -1.20952 0.898374 -1.81996 0.503379 0.918547 0.297399\n", + " 0.255088 0.617499 0.280616 0.55378 0.434486 0.493207" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "R = randn(20,20) # a random 20x20 matrix:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "69dc90f2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "20-element Vector{ComplexF64}:\n", + " -4.029758270745915 + 0.0im\n", + " -2.6633757669626856 - 1.3053111840992218im\n", + " -2.6633757669626856 + 1.3053111840992218im\n", + " -2.1982876157312465 + 0.0im\n", + " -2.1438324706632796 - 2.3199797501785926im\n", + " -2.1438324706632796 + 2.3199797501785926im\n", + " -0.052050103895014396 - 1.7072423620543062im\n", + " -0.052050103895014396 + 1.7072423620543062im\n", + " 0.17696831746269315 - 3.9684099248490288im\n", + " 0.17696831746269315 + 3.9684099248490288im\n", + " 0.9903843645924151 - 0.4121854874725685im\n", + " 0.9903843645924151 + 0.4121854874725685im\n", + " 1.469730788006183 - 2.31275995485107im\n", + " 1.469730788006183 + 2.31275995485107im\n", + " 1.7756728065541905 + 0.0im\n", + " 2.195870725270419 - 3.853100900798987im\n", + " 2.195870725270419 + 3.853100900798987im\n", + " 2.997204454240207 - 2.006535728143899im\n", + " 2.997204454240207 + 2.006535728143899im\n", + " 3.76663409591321 + 0.0im" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "λ = eigvals(R) # its complex eigenvalues" + ] + }, + { + "cell_type": "markdown", + "id": "c7365336", + "metadata": {}, + "source": [ + "One thing to notice is that the complex numbers come in complex-conjugate pairs. This is *always* the case for real $A$, since\n", + "$$\n", + "\\overline{Ax = \\lambda x} = \\bar{A} \\bar{x} = \\bar{\\lambda}\\bar{x} = A\\bar{x}\n", + "$$\n", + "i.e. $\\bar{x}$ is also an eigenvector with eigenvalue $\\bar{\\lambda}$." + ] + }, + { + "cell_type": "markdown", + "id": "0a6f7499", + "metadata": {}, + "source": [ + "Let's plot the eigenvalues in the complex plane:" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "5c6d984a", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGZCAYAAAAUzjLvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIJ0lEQVR4nO3deVhUZf8/8PdhZN9kEQVBUTAVFfdMEoUUUXFPS+2rgmbpo4VlZrapPfVQ6pOauVUGLila4poapKBY+rhnmVpugYApUoCoIOP9+2N+TI7DzgxnhvN+XddcOvfcc87nzMxhPnMv55aEEAJERESkWBZyB0BERETyYjJARESkcEwGiIiIFI7JABERkcIxGSAiIlI4JgNEREQKx2SAiIhI4ZgMEBERKRyTASIiIoVjMmAGrl69CkmSEBcXJ3coRiNJEubOnSt3GDW2adMmtGnTBra2tpAkCadPn5Y7JKNJSUmBJElISUmROxQAQFZWFt5++210794d7u7ucHJyQufOnfHZZ59BrVbr1b99+zamT58OLy8v2NjYoEOHDoiPjzdILGq1Gh9//DH69esHb29v2NnZoXXr1njjjTfw999/l/qcpUuXolWrVrC2tkazZs0wb9483L9/3yDxGMry5cur/HfI1D4nVAZBJu/evXvi8OHD4saNG3KHYjQAxJw5c+QOo0Zu3LghLC0txaBBg0RKSoo4fPiwKCgokDsso0lOThYARHJystyhCCGE2Llzp/Dx8RFvvfWW+Pbbb0ViYqJ45ZVXhIWFhYiKitKrHxYWJurXry9Wrlwp9u/fL55//nkBQHz11Vc1jiU/P184OjqKF154QXz99dciOTlZ/Pe//xUuLi4iICBA3LlzR6f++++/LyRJErNnzxbJycli/vz5wsrKSkyaNKnGsRhSmzZtRK9evar0nNzcXHH48GGRm5trnKDIIJgMkEmoC8nAoUOHBACxadMmWfb/6BeMsZlaMpCTkyOKior0yqdOnSoAiLS0NG3Zt99+KwCIDRs26NQNCwsTXl5eori4WKc8NTVVXLt2rdT93rt3T2zbtk2nrLi4WGRnZ+vV/frrrwUAsW7dOm1Zdna2sLGxES+88IJO3Q8++EBIkiTOnj1bxhHXvqokA0VFReL+/fvGDYgMht0EMvv9998xZswYeHh4wNraGq1bt8ayZct06pTVTbB9+3YEBgbC2toazZs3x5IlSzB37lxIkqRTTwiB5cuXo0OHDrC1tYWLiwtGjBiBy5cv69QLCQlB27ZtcezYMQQHB8POzg7NmzfHhx9+iAcPHgAAbt68CSsrK7zzzjt6x3L+/HlIkoRPPvlEW/df//oXAgIC4ODgAA8PDzz11FNITU2t8HUp7TgAIC4uDpIk4erVqzrlmzZtQvfu3WFvbw8HBweEh4fj1KlTOnUuX76MUaNGwcvLC9bW1mjYsCF69+5dqab8HTt2oHv37rCzs4OjoyPCwsJw+PBh7eORkZHo0aMHAODZZ5+FJEkICQkpc3s1eW0AwNfXFwMHDkRCQgI6duwIGxsbzJs3DwCwbNky9OzZEx4eHrC3t0e7du0wf/58vSbnyrzfJc6fP49+/frBzs4O7u7umDx5MvLz80uN7csvv0T79u1hY2MDV1dXDBs2DOfOndOpExkZCQcHB5w/fx7h4eGwt7eHp6cnPvzwQwDAkSNH0KNHD9jb2+Oxxx7DmjVrKnxNXFxcYGlpqVf++OOPAwCuXbumLdu6dSscHBwwcuRInbpRUVHIzMzE//73P23ZgwcPMG3aNPTs2RN//PGHTv07d+5g4MCBGDduHK5fv64tV6lUcHNzKzOW9PR0bdnevXtx7949REVF6cUihMC2bdvKPe6Sc2L//v2YNGkS3Nzc4OTkhHHjxqGgoADXr1/HM888g/r168PT0xOvvfaa3mdh3rx56NatG1xdXeHk5IROnTph9erVEA+tY+fr64uzZ8/iwIEDkCQJkiTB19cXwD9dAevWrcOMGTPQuHFjWFtb4+LFi3rdBNnZ2fDx8UFQUJBOHL/++ivs7e0xduzYco+XjETmZETRzp49K5ydnUW7du3E2rVrRWJiopgxY4awsLAQc+fO1da7cuWKACBiY2O1ZXv27BEWFhYiJCREbN26VXz99deiW7duwtfXVzz6tk6aNElYWlqKGTNmiL1794oNGzaIVq1aiYYNG4rr169r6/Xq1Uu4ubmJFi1aiJUrV4qkpCTxr3/9SwAQa9as0dYbNmyY8PHxEWq1Wmc/r7/+urCystL+Ijp//ryYMmWKiI+PFykpKWLXrl1i4sSJwsLCQu/XJB5pGZgzZ47ecQghRGxsrAAgrly5oi0r+QU1YcIEsWvXLpGQkCC6d+8u7O3tdX5VtWzZUvj7+4t169aJAwcOiC1btogZM2ZU+Mv2q6++EgBE3759xbZt28SmTZtE586dhZWVlUhNTRVCCHHx4kWxbNkyAUD85z//EYcPHy73F11VXpvSNG3aVHh6eormzZuLL7/8UiQnJ4ujR48KIYR45ZVXxIoVK8TevXvF/v37xaJFi4S7u7teU3ll3+/r168LDw8P0bhxYxEbGyt2794tnnvuOdGkSRO9loH//Oc/AoAYPXq0+Pbbb8XatWtF8+bNhbOzs/jtt9+09caPHy+srKxE69atxZIlS0RSUpKIiooSAMTs2bPFY489JlavXi2+++47MXDgQAFAHD9+vMLXpTTjx48X9erV0/ml/sQTT4iuXbvq1f3ll18EALFq1Sqd8qysLBEQECCaNGkifv/9dyGEpiugZ8+eon79+uJ///tfpWIp+fxu375dW/bGG28IAOL27dt69d3d3cXo0aMrtc1mzZqJGTNmiMTERPHRRx8JlUolRo8eLTp16iTef/99kZSUJGbNmiUAiP/+978624iMjBSrV68WSUlJIikpSfz73/8Wtra2Yt68edo6J0+eFM2bNxcdO3YUhw8fFocPHxYnT54UQvzTStS4cWMxYsQIsWPHDrFr1y5x69atUluQDh06JOrVqydeeeUVIYQQBQUFIiAgQLRq1arU14GMj8mAjMLDw4W3t7deX9q0adOEjY2NyMnJEUKUngx07dpV+Pj4iMLCQm1Zfn6+cHNz0/kSPXz4cKknf3p6urC1tRWvv/66tqxXr14CgN4ftoCAABEeHq69v2PHDgFAJCYmasuKi4uFl5eXePrpp8s83uLiYnH//n3Ru3dvMWzYMJ3HqpsMpKWliXr16omXXnpJp15+fr5o1KiReOaZZ4QQmqZYAGLx4sVlxlcatVotvLy8RLt27XSSn/z8fOHh4SGCgoK0ZSV/9L7++usq7UOI8l+b0jRt2lSoVCpx4cKFCuO/f/++WLt2rVCpVNrPlBCVf79nzZolJEkSp0+f1qkXFham80f+r7/+Era2tmLAgAE69dLS0oS1tbUYM2aMtmz8+PECgNiyZYu27P79+6JBgwYCgPZLRgghbt26JVQqlXj11VcreFX0fffdd8LCwkL7pVOiRYsWOsdYIjMzU5vQPerGjRuiffv2wsvLSxw+fFg88cQTwt3dXSfW8ly7dk00bNhQdOnSReezNGnSJGFtbV3qcx577DHRt2/fcrdbck48eg4MHTpUABAff/yxTnmHDh1Ep06dytxeyWfmvffeE25ubuLBgwfax8rqJij57Pfs2bPMxx5Ncj/66CMBQGzdulWMHz9e2NraijNnzpR7rGQ87CaQyb1797Bv3z4MGzYMdnZ2KC4u1t4GDBiAe/fu4ciRI6U+t6CgAMePH8fQoUNhZWWlLXdwcMCgQYN06u7atQuSJOH//u//dPbRqFEjtG/fXm+Eb6NGjbRNmSUCAwN1mkf79++PRo0aITY2Vlv23XffITMzExMmTNB57sqVK9GpUyfY2NigXr16sLS0xL59+/Sajavru+++Q3FxMcaNG6dzfDY2NujVq5f2+FxdXeHn54cFCxbg448/xqlTp/Sawktz4cIFZGZmYuzYsbCw+Od0cXBwwNNPP40jR47gzp071Yq9pq9NYGAgHnvsMb3yU6dOYfDgwXBzc4NKpYKlpSXGjRsHtVqN3377TaduZd7v5ORktGnTBu3bt9epN2bMGJ37hw8fxt27dxEZGalT7uPjg6eeegr79u3TKZckCQMGDNDer1evHvz9/eHp6YmOHTtqy11dXeHh4aHXRF+RkydP4plnnsETTzyBmJgYvcdL64Yq77EGDRogOTkZXl5e6N69O65cuYKUlBSdWMuSk5ODAQMGQAiBTZs26XyWqhNLaQYOHKhzv3Xr1gCAiIgIvfJHX8v9+/ejT58+cHZ21n5m3n33Xdy6dQs3btyo1P4B4Omnn6503ZkzZyIiIgKjR4/GmjVrsHTpUrRr167SzyfDYjIgk1u3bqG4uBhLly6FpaWlzq3kD2R2dnapz/3rr78ghEDDhg31Hnu07M8//9TWfXQ/R44c0dtHaf2c1tbWuHv3rvZ+vXr1MHbsWGzdulU7TSouLg6enp4IDw/X1vv4448xZcoUdOvWDVu2bMGRI0dw7Ngx9OvXT2d7NfHnn38CALp27ap3fJs2bdIenyRJ2LdvH8LDwzF//nx06tQJDRo0wMsvv1xm3zegeZ8AwNPTU+8xLy8vPHjwAH/99VeV4zbEa1NaTGlpaQgODkZGRgaWLFmC1NRUHDt2TDsO5dFtV+b9vnXrFho1aqRX79Gyil6rksdL2NnZwcbGRqfMysoKrq6ues+3srLCvXv39MrLcurUKYSFhaFFixbYvXs3rK2tdR53c3PTiwfQfGkDKDUGACgsLMTdu3dhYWGBoqKiSiWCf/31F8LCwpCRkYGkpCQ0b95cL5Z79+6Vuq2cnJwyY3nUo/VKfiiUVv7wa3n06FH07dsXAPD555/jhx9+wLFjx/DWW28B0P/MlKe0974skiQhMjIS9+7dQ6NGjThWQGb15A5AqVxcXKBSqTB27FhMnTq11DrNmjUr87mSJGm/CB/28CAmAHB3d4ckSUhNTdX7gwig1LLKiIqKwoIFCxAfH49nn30WO3bswPTp06FSqbR11q9fj5CQEKxYsULnueV9+ZYo+ZIoLCzUifHR5MXd3R0A8M0336Bp06blbrNp06ZYvXo1AOC3337D5s2bMXfuXBQVFWHlypWlPqfkyzIrK0vvsczMTFhYWMDFxaXC43lUTV6bEqX9Yty2bRsKCgqQkJCg83rU5HoHbm5uep8rQP+zVtFrVfJeGdupU6fQp08fNG3aFImJiXB2dtar065dO2zcuBHFxcWoV++fP4M///wzAKBt27Z6z7l27Rp69+6N4uJinDlzBtHR0ejTpw92796NJ598stRY/vrrL/Tp0wdXrlzBvn37EBgYWGosJfvu1q2btvz69evIzs4uNRZDio+Ph6WlJXbt2qWTnFU0cLE0lW3FADSfk6lTp6JDhw44e/YsXnvtNe3gY6p9bBmQiZ2dHUJDQ3Hq1CkEBgaiS5cuerfSfrUBgL29Pbp06YJt27ahqKhIW3779m3s2rVLp+7AgQMhhEBGRkap+6hus1zr1q3RrVs3xMbGYsOGDSgsLNQbDS1Jkl6ycebMGZ1R+GUpGaV85swZnfKdO3fq3A8PD0e9evVw6dKlUo+vS5cupW7/sccew9tvv4127drh5MmTZcbRsmVLNG7cGBs2bNAZWV1QUIAtW7ZoZxhUVU1em4q2C+gmeUIIfP7559XeZmhoKM6ePYuffvpJp3zDhg0697t37w5bW1usX79ep/zatWvYv38/evfuXe0YKuv06dPo06cPvL29kZSUVGaiNmzYMNy+fRtbtmzRKV+zZg28vLx0vpQBzYyenj17QpIkHDx4EG3atMGuXbsQHByM8PBwJCcn6+2jJBG4fPkyEhMTy+xO6NevH2xsbPRmC5XMEhg6dGjlX4BqkCQJ9erV00nk7969i3Xr1unVfbTVqLrUajVGjx4NSZKwZ88exMTEYOnSpUhISKjxtql62DIgoyVLlqBHjx4IDg7GlClT4Ovri/z8fFy8eBE7d+7E/v37y3zue++9h4iICISHhyM6OhpqtRoLFiyAg4ODtqkTAJ588km88MILiIqKwvHjx9GzZ0/Y29sjKysLhw4dQrt27TBlypRqxT9hwgS8+OKLyMzMRFBQEFq2bKnz+MCBA/Hvf/8bc+bMQa9evXDhwgW89957aNasGYqLi8vd9oABA+Dq6oqJEyfivffeQ7169RAXF6czJQvQJA3vvfce3nrrLVy+fBn9+vWDi4sL/vzzTxw9ehT29vaYN28ezpw5g2nTpmHkyJFo0aIFrKyssH//fpw5cwZvvPFGmXFYWFhg/vz5eO655zBw4EC8+OKLKCwsxIIFC/D3339rp8JVVU1em/KEhYXBysoKo0ePxuuvv4579+5hxYoV1erKKDF9+nR8+eWXiIiIwPvvv4+GDRviq6++wvnz53Xq1a9fH++88w7efPNNjBs3DqNHj8atW7cwb9482NjYYM6cOdWOoTIuXLiAPn36AAA++OAD/P777/j999+1j/v5+aFBgwYANONewsLCMGXKFOTl5cHf3x8bN27E3r17sX79ep0vxgcPHmDw4MFwcHBAUlKStivOxsYGW7duxejRozFo0CBcvHhR23Vy9+5d7fTWxYsXo7i4WGcMUIMGDeDn5wdA04z/9ttv45133oGrqyv69u2LY8eOYe7cuXj++ecREBBg1NctIiICH3/8McaMGYMXXngBt27dwsKFC0ttNWzXrh3i4+OxadMmNG/eHDY2NtX6QTFnzhykpqYiMTERjRo1wowZM3DgwAFMnDgRHTt2LLNVlIxIztGLpJkpMGHCBNG4cWNhaWkpGjRoIIKCgsT777+vUwePzCYQQoitW7eKdu3aCSsrK9GkSRPx4Ycfipdfflm4uLjo7efLL78U3bp1E/b29sLW1lb4+fmJcePG6UzX6tWrl2jTpo3ec8ePHy+aNm2qV56bmytsbW0FAPH555/rPV5YWChee+010bhxY2FjYyM6deoktm3bVur2UMpFh44ePSqCgoKEvb29aNy4sZgzZ4744osv9KYWCiHEtm3bRGhoqHBychLW1taiadOmYsSIEeL7778XQgjx559/isjISNGqVSthb28vHBwcRGBgoFi0aJHeBWZKs23bNtGtWzdhY2Mj7O3tRe/evcUPP/ygU6cqswmq8tqUpmnTpiIiIqLUx3bu3Cnat28vbGxsROPGjcXMmTPFnj179EZ0V+X9/vXXX0VYWJiwsbERrq6uYuLEiWL79u2ljhL/4osvRGBgoLCyshLOzs5iyJAhetMsx48fL+zt7fX2XVZM5R1viZJR9WXdHj1/8vPzxcsvvywaNWokrKysRGBgoNi4cWOp2z59+nSpFxESQjMT5MCBAzplJedsWbfx48frbWfJkiXiscce057Pc+bMKfUiSmUd97Fjx3TKS2bk3Lx5U6e8tNf+yy+/FC1bthTW1taiefPmIiYmRqxevVrvXLt69aro27evcHR0FAC0n5PyPvuPziZITEwUFhYWeuf7rVu3RJMmTUTXrl11ZklR7ZCEeKjtk8za/fv30aFDBzRu3BiJiYlyh0NERGaC3QRmbOLEiQgLC4OnpyeuX7+OlStX4ty5c1iyZIncoRERkRlhMmDG8vPz8dprr+HmzZuwtLREp06dsHv3bm2/KRERUWWwm4CIiEjhOLWQiIhI4ZgMEBERKRyTASKFE0IgLy8P7DEkUi4mA0QKl5+fD2dn5ypdCpmI6hYmA0RERArHZICIiEjhmAwQEREpHJMBIiIihWMyQEREpHBMBoiIiBSOyQAREZHCMRkgIiJSOCYDRERECsdkgIiISOGYDBARESkckwGiOiQmJgaSJGH69Om1ul+1GkhJATZu1PyrVtfq7omohurJHQARGcaxY8fw2WefITAwsFb3m5AAREcD1679U+btDSxZAgwfXquhEFE1sWWAqA64ffs2nnvuOXz++edwcXGptf0mJAAjRugmAgCQkaEpT0iotVCIqAaYDBDVAVOnTkVERAT69OlTYd3CwkLk5eXp3KpDrda0CAih/1hJ2fTp7DIgMgdMBojMXHx8PE6ePImYmJhK1Y+JiYGzs7P25uPjU639pqbqtwg8TAggPV1Tj4hMG5MBIjOWnp6O6OhorF+/HjY2NpV6zuzZs5Gbm6u9paenV2vfWVmGrUemhYNClYUDCInM2IkTJ3Djxg107txZW6ZWq3Hw4EF8+umnKCwshEql0nmOtbU1rK2ta7xvT0/D1iPTwUGhyiMJUVqPHxGZg/z8fPzxxx86ZVFRUWjVqhVmzZqFtm3bVriNvLw8ODs7Izc3F05OTpXet1oN+PpqBguW9ldEkjRfIFeuAI/kI2TCSgaFPvqeSpLm32++YUJQF7FlgMiMOTo66n3h29vbw83NrVKJQE2oVJpfiiNGaL4oHv7yKPniWLyYiYA5qWhQqCRpBoUOGcL3ta7hmAEiqrbhwzW/FBs31i339uYvSHPEQaHKxZYBojomJSWlVvc3fLjml2JqqmawoKcnEBzMX47miINClYvJABHVmEoFhITIHQXVFAeFKhe7CYiICICmRcfb+58xH4+SJMDHR1OP6hYmA0REBOCfQaGAfkLAQaF1G5MBIiLS4qBQZeJ1BogUrrrXGaC6Ta3moFAl4QBCIiLSw0GhysJuAiIiIoVjMkBERKRwTAaIjCgzMxMLFy6UOwwionIxGSAysNu3byMuLg59+vRBkyZNsGrVKrlDIiIqF5MBIgNQq9XYtWsXRo0aBQ8PD7z22mto0aIFDh48iN9//13u8IiIysXZBEQ1cOTIEaxbtw5ff/01bt26hd69e2PDhg2IiIiApaWl3OEREVUKWwaIqik6Oho9evTAhQsXMHHiRHh5eeF///sfDhw4gMzMTLnDIyKqNCYDRNW0evVqbNq0Cd9//z1iYmJw9epVrFy5EkeOHIG/vz9GjhyJI0eOyB0mEVGFmAwQVdPIkSPRt29f7X2VSoXRo0fj8OHDOHjwIFQqFXr27ImgoCAZoyQiqhgvR0xkRNeuXcOyZcsQExMjdyhl4uWIiYjJAJHCMRkgInYTEBERKRyTASIiIoVjMkBERKRwTAaIqmn48OHIy8sDAKxduxaFhYUyR0REVD0cQEhUTVZWVvjjjz/g6ekJlUqFrKwseHh4yB1WlXEAIRHxcsRE1dSqVSvMnj0boaGhEEJg8+bNZX6Zjhs3rpajIyKqPLYMEFXTjz/+iFdffRWXLl1CTk4OHB0dIUmSXj1JkpCTkyNDhJXDlgEiYjJAZAAWFha4fv06uwmIyCxxACGRAVy5cgUNGjSQOwwiomrhmAEiA2jatCn+/vtvrF69GufOnYMkSWjdujUmTpwIZ2dnucMjIioXuwmIDOD48eMIDw+Hra0tHn/8cQghcPz4cdy9exeJiYno1KmT3CGWid0ERMRkgMgAgoOD4e/vj88//xz16mka3IqLi/H888/j8uXLOHjwoMwRlo3JABExGSAyAFtbW5w6dQqtWrXSKf/111/RpUsX3LlzR6bIKsZkgIg4gJDIAJycnJCWlqZXnp6eDkdHRxkiIiKqPCYDRAbw7LPPYuLEidi0aRPS09Nx7do1xMfH4/nnn8fo0aPlDo+IqFycTUBkAAsXLoQkSRg3bhyKi4sBAJaWlpgyZQo+/PBDo+13xYoVWLFiBa5evQoAaNOmDd59913079/faPskorqHYwaIDOjOnTu4dOkShBDw9/eHnZ2dUfe3c+dOqFQq+Pv7AwDWrFmDBQsW4NSpU2jTpk2ltsExA0TEZICojnF1dcWCBQswceLEUh8vLCzUWWExLy8PPj4+TAaIFIxjBojqCLVajfj4eBQUFKB79+5l1ouJiYGzs7P25uPjU4tREpEpYssAkZn7+eef0b17d9y7dw8ODg7YsGEDBgwYUGZ9tgwQ0aOYDBCZuaKiIqSlpeHvv//Gli1b8MUXX+DAgQMICAio1PM5ZoCI2E1AZAAHDx7UziJ4WHFxsdGvPmhlZQV/f3906dIFMTExaN++PZYsWWLUfRJR3cJkgMgAQkNDkZOTo1eem5uL0NDQWo1FCKHTDUBEVBFeZ4DIAIQQkCRJr/zWrVuwt7c32n7ffPNN9O/fHz4+PsjPz0d8fDxSUlKwd+9eo+2TiOoeJgNENTB8+HAAgCRJiIyMhLW1tfYxtVqNM2fOICgoyGj7//PPPzF27FhkZWXB2dkZgYGB2Lt3L8LCwoy2TyKqe5gMENWAs7MzAE3LgKOjI2xtbbWPWVlZ4YknnsCkSZOMtv/Vq1cbbdtEpBxMBohqIDY2FkIICCGwdOlSLkpERGaJAwiJakgIgQ0bNuD69etyh0JEVC1MBohqyMLCAi1atMCtW7fkDoWIqFqYDBAZwPz58zFz5kz88ssvcodCRFRlvAIhkQG4uLjgzp07KC4uhpWVlc5AQgClXoPAVPAKhETEAYREBrB48WK5QyAiqja2DBApHFsGiIgtA0QGdvfuXdy/f1+njF+yRGTKOICQyAAKCgowbdo0eHh4wMHBAS4uLjo3IiJTxmSAyABef/117N+/H8uXL4e1tTW++OILzJs3D15eXli7dq3c4RERlYtjBogMoEmTJli7di1CQkLg5OSEkydPwt/fH+vWrcPGjRuxe/duuUMsE8cMEBFbBogMICcnB82aNQOgGR9QMpWwR48eOHjwoJyhERFViMkAkQE0b94cV69eBQAEBARg8+bNAICdO3eifv368gVGRFQJTAaIDCAqKgo//fQTAGD27NnasQOvvPIKZs6cKXN0RETl45gBIiP4448/cOLECfj5+aF9+/Zyh1MujhkgIiYDRArHZICI2E1AVaJWAykpwMaNmn/VarkjMh379u3DwIED4efnB39/fwwcOBDff/+93GGRieE5RKaIyQBVWkIC4OsLhIYCY8Zo/vX11ZQr3aeffop+/frB0dER0dHRePnll+Hk5IQBAwbg008/lTs8MhE8h8hUsZuAKiUhARgxAnj00yJJmn+/+QYYPrz24zIVjRs3xuzZszFt2jSd8mXLluGDDz5AZmamTJFVjN0EtYPnEJkyJgNUIbVa8+vl2rXSH5ckwNsbuHIFUKlqNTST4ejoiFOnTsHf31+n/Pfff0fHjh1x+/ZtmSKrGJMB4+M5RKaO3QRUodTUsv+IAZpfOunpmnpKNXjwYGzdulWvfPv27Rg0aJAMEZEp4TlEpo6rFlKFsrIMW68uat26NT744AOkpKSge/fuAIAjR47ghx9+wIwZM/DJJ59o67788styhUky4TlEpo7dBFShlBTNQKeKJCcDISHGjsY0lVyKuCKSJOHy5ctGjqZq2E1gfDyHyNQxGaAKlfR3ZmToD34C2N9p7pgMGB/PITJ1HDNAFVKpgCVLNP8vGflcouT+4sX8I0ZUFp5DZOqYDFClDB+umfrUuLFuubc3p0QRVQbPITJl7CagKlGrNSOes7IAT08gOJi/ZswduwlqF88hMkVMBogUjskAEbGbgIiISOGYDBAZgK+vL9577z2kpaXJHQqRonEhqOphMkBkADNmzMD27dvRvHlzhIWFIT4+HoWFhXKHRaQoXAiq+pgMEBnASy+9hBMnTuDEiRMICAjAyy+/DE9PT0ybNg0nT5402n5jYmLQtWtXODo6wsPDA0OHDsWFCxeMtj8iU1WyENSjl33OyNCUMyEoHwcQEhnB/fv3sXz5csyaNQv3799H27ZtER0djaioKEiPTjSvgX79+mHUqFHo2rUriouL8dZbb+Hnn3/Gr7/+Cnt7+0ptgwMIydxxIaiaYzJAZED379/H1q1bERsbi6SkJDzxxBOYOHEiMjMz8emnnyI0NBQbNmww2v5v3rwJDw8PHDhwAD179qzUc5gMkLnj5Z5rjgsVERnAyZMnERsbi40bN0KlUmHs2LFYtGgRWrVqpa3Tt2/fSn9BV1dubi4AwNXVtcw6hYWFOuMZ8vLyjBoTkbFxIaiaYzJAZABdu3ZFWFgYVqxYgaFDh8LS0lKvTkBAAEaNGmW0GIQQePXVV9GjRw+0bdu2zHoxMTGYN2+e0eIgqm2enoatp0TsJiCqIbVajXXr1mHw4MHl/iI3tqlTp+Lbb7/FoUOH4O3tXWa90loGfHx82E1AZosLQdUcZxMQ1ZBKpcLkyZO1TfRyeOmll7Bjxw4kJyeXmwgAgLW1NZycnHRuROaMC0HVHJMBIgNo164dLl++XOv7FUJg2rRpSEhIwP79+9GsWbNaj4HIFHAhqJphNwGRASQmJmLWrFn497//jc6dO+tN6zPWr+9//etf2LBhA7Zv346WLVtqy52dnWFra1upbXA2AdUlXAiqepgMEBmAhcU/jWwPX0dACAFJkqA20jVRy7pmQWxsLCIjIyu1DSYDRMTZBEQGkJycLMt+mcsTkSGwZYBI4R5tGWAzK5HysGWAyIDu3LmDtLQ0FBUV6ZQHBgbKFFHVJCQA0dG6l3X19taM1OYALKK6iy0DRAZw8+ZNREVFYc+ePaU+bqwxA4ZQ0jKwbl0uxo1z0punXTIsgSOyieouTi0kMoDp06fjr7/+wpEjR2Bra4u9e/dizZo1aNGiBXbs2CF3eJUya1bpF2wpKZs+nWvDE9VV7CYgMoD9+/dj+/bt6Nq1KywsLNC0aVOEhYXByckJMTExiIiIkDvECmVmlv2YEEB6umYsARd6Iap72DJAZAAFBQXw8PAAoFkk6ObNmwA0FyM6efKknKEZFBd6IaqbmAwQGUDLli1x4cIFAECHDh2watUqZGRkYOXKlfCsQ6uj1KFDIaKHsJuAyACmT5+OrP//s3nOnDkIDw/HV199BSsrK8TFxckbXCV5eWl++Ze30EtwcO3HRUTGx9kEREZw584dnD9/Hk2aNIG7u7vc4ZTr0dkEgG5CwNkERHUfuwmIjMDOzg6dOnUy+UTgYYMHc6EXIqViywCRAajVasTFxWHfvn24ceMGHjx4oPP4/v37ZYqsYrwCIRFxzACRAURHRyMuLg4RERFo27ZtmQsImQOVitMHiZSGLQNEBuDu7o61a9diwIABcodSZVy1kIg4ZoDIAKysrODv7y93GERE1cJkgMgAZsyYgSVLlnBJYSIySxwzQGQAhw4dQnJyMvbs2YM2bdrA0tJS5/GEhASZIiMiqhiTASIDqF+/PoYNGyZ3GERE1cIBhGaK07/IUDiAkIjYMmCGEhKA6Gjg2rV/yry9gSVLeGEYIiKqOrYMmJmEBGDECP3rx/OSsbWvU6dO2LdvH1xcXNCxY8dyry1gyisXsmWAiNgyYEbUak2LQGnpmxCahGD6dGDIEHYZ1IYhQ4bA2toaADB06FB5gyEiqgG2DJiRlBQgNLTiesnJvIIcVR5bBoiI1xkwI/9/hVyD1SMiIgLYTWBWPD0NW48Mx8XFpdQxA5IkwcbGBv7+/oiMjERUVJQM0RERlY/JgBkJDtbMGsjIKH3cgCRpHg8Orv3YlO7dd9/FBx98gP79++Pxxx+HEALHjh3D3r17MXXqVFy5cgVTpkxBcXExJk2aJHe4REQ6mAyYEZVKM31wxAjNF//DCUHJj9LFizl4UA6HDh3C+++/j8mTJ+uUr1q1ComJidiyZQsCAwPxySefMBkgIpPDAYRmqLTrDPj4aBIBTiuUh4ODA06fPq23WNHFixfRoUMH3L59G5cuXUJgYCAKCgpkirJ0HEBIRBxAaIaGDweuXtXMGtiwQfPvlStMBOTk6uqKnTt36pXv3LkTrq6uAICCggI4OjrWdmhERBViN4GZUqk4fdCUvPPOO5gyZQqSk5Px+OOPQ5IkHD16FLt378bKlSsBAElJSejVq5fMkRIR6WM3AZGB/PDDD/j0009x4cIFCCHQqlUrvPTSSwgKCpI7tHKxm4CImAwQKRyTASLimAEiA7t79y7y8vJ0bsZ08OBBDBo0CF5eXpAkCdu2bTPq/oio7mEyQGQAd+7cwbRp0+Dh4QEHBwe4uLjo3IypoKAA7du3x6effmrU/RBR3cUBhEQGMHPmTCQnJ2P58uUYN24cli1bhoyMDKxatQoffvihUffdv39/9O/f36j7IKK6jckAkQHs3LkTa9euRUhICCZMmIDg4GD4+/ujadOm+Oqrr/Dcc8/JHaJWYWEhCgsLtfeN3Y1BRKaP3QREBpCTk4NmzZoBAJycnJCTkwMA6NGjBw4ePChnaHpiYmLg7Oysvfn4+MgdEhHJjMkAkQE0b94cV69eBQAEBARg8+bNADQtBvXr15cvsFLMnj0bubm52lt6errcIRGRzNhNQGQAUVFR+Omnn9CrVy/Mnj0bERERWLp0KYqLi/Hxxx/LHZ4Oa2trWFtbyx0GEZkQJgNEBvDKK69o/x8aGorz58/j+PHj8PPzQ/v27WWMjIioYkwGiIygSZMmaNKkSa3s6/bt27h48aL2/pUrV3D69Gm4urrWWgxEZN54BUIiAzl69ChSUlJw48YNPHjwQOcxY3YVpKSkIDQ0VK98/PjxiIuLq/D5vAIhEbFlgMgA/vOf/+Dtt99Gy5Yt0bBhQ0iSpH3s4f8bQ0hICJjTE1FNsGWAyAAaNmyIjz76CJGRkXKHUmVsGSAiTi0kMgALCws8+eSTcodBRFQtTAaIDOCVV17BsmXL5A6DiKha2E1AZAAPHjxAREQEfvvtNwQEBMDS0lLn8YSEBJkiqxi7CYiIAwiJDOCll15CcnIyQkND4ebmZvRBg0REhsSWASIDcHR0RHx8PCIiIuQOpcrYMkBEHDNAZACurq7w8/OTOwwiomphMkBkAHPnzsWcOXNw584duUMhIqoydhMQGUDHjh1x6dIlCCHg6+urN4Dw5MmTMkVWMXYTEBEHEBIZwNChQ+UOgYio2tgyQKRwbBkgIo4ZICIiUjh2ExBVk6urK3777Te4u7vDxcWl3GsL5OTk1GJkRERVw2SAqJoWLVoER0dHAMDixYvlDYaqTa0GUlOBrCzA0xMIDgZUKrmjIqpdHDNApHBKHjOQkABERwPXrv1T5u0NLFkCDB8uX1xEtY1jBohIkRISgBEjdBMBAMjI0JSb8HISRAbHlgEihVNiy4BaDfj66icCJSRJ00Jw5Qq7DEgZ2DJARIqTmlp2IgAAQgDp6Zp6RErAZICIFCcry7D1iMwdkwEiUhxPT8PWIzJ3TAaISHGCgzVjAsq6NIQkAT4+mnpESsBkgIgUR6XSTB8E9BOCkvuLF3PwICkHkwEiUqThw4FvvgEaN9Yt9/bWlPM6A6QknFpIpHBKnFr4MF6BkIiXIyYihVOpgJAQuaMgkhe7CYiIiBSOyQAREZHCMRkgIiJSOCYDRERECsdkgIiISOFknU3AKT2mj+8RkenheVl3mMp7KVsykJAAREfrrhzm7a25Khgv9mEa+B4RmR6el3WHKb2XsnQTJCQAI0boLyGakaEpT0iQIyp6GN8j87J8+XI0a9YMNjY26Ny5M1K59m6dxPOy7jC197LWr0CoVgO+vmWvJS5JmszoyhU2e8mF75F52bRpE8aOHYvly5fjySefxKpVq/DFF1/g119/RZMmTSp8vtKvQGgueF7WHab4XlYqGRBCID8/3yA7TE0FBg6suN6uXVwxTC58j8zLU089hfbt22PRokXasi5dumDgwIGYO3euXv3CwkIUFhZq7+fn5yMgIADp6elMBkwYz8u6Q4730tHREVJZy3SikslAyS8HIiIiMj8VtfyxZaAW5eXlwcfHx+R/gdX0PTKX46wpUzjOrKwstGrVComJiejWrZu2fOHChdiwYQNOnjyp95yqtgyYwnEamzkcoyH+dprDcRqCqR+nob4Hq3KcFbUMVGo2gSRJBntB+/XT9IVkZAClpSElfSX9+tXdfi8nJyeT/ICWMNR7ZOrHaShyHuft27cBAA4ODjoxWFtbo169epWKKy8vD0DFx6GE99OUj9GQfztN+TgNyVSP09Dfg4Y4zlqfTaBSaaZNAJoDfpgkaV6VxYvrbiJgDsp/jzT/8j0yDe7u7lCpVLh+/bpO+Y0bN9CwYUOZoiJj4HlZd5jieynL1MLhw4FvvgEaN9Yt9/IS+OYbzpU1BWW9R97e4HtkQqysrNC5c2ckJSXplCclJSEoKEimqMhYeF7WHab2Xsp20aHhw4EhQzR9J2lp97F//1dYvnw07Oys5QrJ6KytrTFnzhxYW5vHMT78HlXl6ljmdpzVZSrH+eqrr2Ls2LHo0qULunfvjs8++wxpaWmYPHmyQbZvKsdpTOZ0jNU9LwHzOs6aMJfjrMl7CRj2OGv9OgNEZHjLly/H/PnzkZWVhbZt22LRokXo2bNnpZ7L6wwQEZMBIoVjMkBEXLWQiIhI4ZgMEBERKRyTASIiIoVjMkBERKRwJp0MFBYWokOHDpAkCadPn5Y7HIMaPHgwmjRpAhsbG3h6emLs2LHIzMyUOyyDunr1KiZOnIhmzZrB1tYWfn5+mDNnDoqKiuQOzeA++OADBAUFwc7ODvXr15c7HINRwtLIBw8exKBBg+Dl5QVJkrBt2za5QzK4mJgYdO3aFY6OjvDw8MDQoUNx4cIFucMyuBUrViAwMFB7Rb7u3btjz549codlVDExMZAkCdOnT6/Rdkw6GXj99dfh5eUldxhGERoais2bN+PChQvYsmULLl26hBEjRsgdlkGdP38eDx48wKpVq3D27FksWrQIK1euxJtvvil3aAZXVFSEkSNHYsqUKXKHYjCbNm3C9OnT8dZbb+HUqVMIDg5G//79kZaWJndoBpWXV4D69Yfi2We3AegFtVruiAzvwIEDmDp1Ko4cOYKkpCQUFxejb9++KCgokDs0g/L29saHH36I48eP4/jx43jqqacwZMgQnD17Vu7QjOLYsWP47LPPEBgYWPONCRO1e/du0apVK3H27FkBQJw6dUrukIxq+/btQpIkUVRUJHcoRjV//nzRrFkzucMwmtjYWOHs7Cx3GFWSm5srAIjc3Fyd8scff1xMnjxZp6xVq1bijTfeqM3wjGrLFiG8vYXQXCFec3NzuyO2bJE7MuO6ceOGACAOHDggdyhG5+LiIr744gu5wzC4/Px80aJFC5GUlCR69eoloqOja7Q9k2wZ+PPPPzFp0iSsW7cOdnZ2codjdDk5Ofjqq68QFBQES0tLucMxqtzcXLi6usodBlWgqKgIJ06cQN++fXXK+/btix9//FGmqAwrIQEYMQK4dk23/NYtG4wYoXm8rsrNzQWAOn0uqtVqxMfHo6CgAN27d5c7HIObOnUqIiIi0KdPH4Nsz+SSASEEIiMjMXnyZHTp0kXucIxq1qxZsLe3h5ubG9LS0rB9+3a5QzKqS5cuYenSpQa7TC4ZT3Z2NtRqtd5iRw0bNtRbFMkcqdVAdHTpK8YBmpVipk9HnewyEELg1VdfRY8ePdC2bVu5wzG4n3/+GQ4ODrC2tsbkyZOxdetWBAQEyB2WQcXHx+PkyZOIiYkx2DZrLRmYO3cuJEkq93b8+HEsXboUeXl5mD17dm2FZjCVPcYSM2fOxKlTp5CYmAiVSoVx48ZBmMEFIat6nACQmZmJfv36YeTIkXj++edlirxqqnOcdc2j658LIcpdE91cpKbqtwg8TAggPV1Tr66ZNm0azpw5g40bN8odilG0bNkSp0+fxpEjRzBlyhSMHz8ev/76q9xhGUx6ejqio6Oxfv162NjYGGy7tXY54uzsbGRnZ5dbx9fXF6NGjcLOnTt1/uCo1WqoVCo899xzWLNmjbFDrbbKHmNpb+C1a9fg4+ODH3/80eSbtKp6nJmZmQgNDUW3bt0QFxcHCwuTa5AqVXXez7i4OEyfPh1///23kaMznNIuR1xUVAQ7Ozt8/fXXGDZsmLZudHQ0Tp8+jQMHDsgVrkFs3AiMGVNxvQ0bgNGjjR9PbXnppZewbds2HDx4EM2aNZM7nFrRp08f+Pn5YdWqVXKHYhDbtm3DsGHDoHpoNSO1Wg1JkmBhYYHCwkKdxyqr1lYtdHd3h7u7e4X1PvnkE7z//vva+5mZmQgPD8emTZvQrVs3Y4ZYY5U9xtKU5GSFhYWGDMkoqnKcGRkZCA0NRefOnREbG2s2iQBQs/fT3D28NPLDyUBSUhKGDBkiY2SG4elp2HqmTgiBl156CVu3bkVKSopiEgFAc+zm8He1snr37o2ff/5ZpywqKgqtWrXCrFmzqpUIADIuYVyWJk2a6Nx3cHAAAPj5+cHb21uOkAzu6NGjOHr0KHr06AEXFxdcvnwZ7777Lvz8/Ey+VaAqMjMzERISgiZNmmDhwoW4efOm9rFGjRrJGJnhpaWlIScnB2lpaVCr1drrYvj7+2s/w+bG2Esjyyk4WLNufEZG6eMGJEmgUSM1goNN7k9ktUydOhUbNmzA9u3b4ejoqB334ezsDFtbW5mjM5w333wT/fv3h4+PD/Lz8xEfH4+UlBTs3btX7tAMxtHRUW+sR8nYsxqNAanRXIRacOXKlTo3tfDMmTMiNDRUuLq6Cmtra+Hr6ysmT54srl27JndoBhUbGysAlHqra8aPH1/qcSYnJ8sdWoXKmloohBDLli0TTZs2FVZWVqJTp051airali1CSJIQwAOdqYWAWgBqERLyidwhGkxZ52FsbKzcoRnUhAkTtJ/XBg0aiN69e4vExES5wzI6Q0wt5BLGRAqn5CWMExI0swoeHkzo4wMsXgwMHy5bWES1jskAkcIpORkANNMHU1OBrCzNGIHgYKCa3a5EZqtudIgREVWTSgWEhMgdBZG8zGdoNxERERkFkwEiIiKFYzJARESkcEwGiIiIFI7JABERkcIxGSAiqgFfX18sXrxY7jBq7Pr16wgLC4O9vT3q169fap24uDhYWlqid+/e2mWQqW5gMkBEZKbi4uLK/OKuqkWLFiErKwunT5/Gb7/9VmqdZ599FkePHsUvv/yCZcuWGWS/ZBqYDBCR4ty/f1/uEGrM0Mdw6dIldO7cGS1atICHh0epdWxtbdGxY0dMmzYN69evN+j+SV5MBoiozps7dy46dOiAL7/8Es2bN4e1tTWEEMjNzcULL7wADw8PODk54amnnsJPP/2kfd6lS5cwZMgQNGzYEA4ODujatSu+//77Ku07MjISQ4cOxbx587T7efHFF1FUVKSts3fvXvTo0QP169eHm5sbBg4ciEuXLmkfv3r1KiRJwubNmxESEgIbGxusX78eUVFRyM3NhSRJkCQJc+fOLTOOFStWwM/PD1ZWVmjZsiXWrVunfczX1xdbtmzB2rVrIUkSIiMjyz2mrl274ty5czh58mSVXgsyXUwGiEgRLl68iM2bN2PLli3aVSUjIiJw/fp17N69GydOnECnTp3Qu3dv5OTkAABu376NAQMG4Pvvv8epU6cQHh6OQYMGIS0trUr73rdvH86dO4fk5GRs3LgRW7duxbx587SPFxQU4NVXX8WxY8ewb98+WFhYYNiwYXjw4IHOdmbNmoWXX34Z586dQ+/evbF48WI4OTkhKysLWVlZeO2110rd/9atWxEdHY0ZM2bgl19+wYsvvoioqCgkJycDAI4dO4Z+/frhmWeeQVZWFpYsWVLu8cTGxgIAWwfqkhovl0REZq28VQvrijlz5ghLS0tx48YNbdm+ffuEk5OTuHfvnk5dPz8/sWrVqjK3FRAQIJYuXaq937RpU7Fo0aIy648fP164urqKgoICbdmKFSuEg4ODUKvVpT7nxo0bAoD4+eefhRD/rN66ePFinXqxsbHC2dm5zH2XCAoKEpMmTdIpGzlypBgwYID2/pAhQ8T48eMr3FZ6erqoV6+eGDJkiGjUqJEoLi6u8Dlk+tgyQESK0LRpUzRo0EB7/8SJE7h9+zbc3Nzg4OCgvV25ckXbRF9QUIDXX38dAQEBqF+/PhwcHHD+/Pkqtwy0b98ednZ22vvdu3fH7du3kZ6eDkDTHTFmzBg0b94cTk5OaNasGQDo7adLly7VOvZz587hySef1Cl78sknce7cuSpva+nSpQgICMBnn32G7Oxs7Nu3r1oxkWnhQkVEpAj29vY69x88eABPT0+kpKTo1S0ZoT9z5kx89913WLhwIfz9/WFra4sRI0bo9PfXhCRJAIBBgwbBx8cHn3/+Oby8vPDgwQO0bdtWbz+PHkN19lVCCKFXVpE7d+7giy++wH//+194eHggPDwc69evR9++fasdF5kGJgNEpEidOnXC9evXUa9ePfj6+pZaJzU1FZGRkRg2bBgAzRiCq1evVnlfP/30E+7evQtbW1sAwJEjR+Dg4ABvb2/cunUL586dw6pVqxAcHAwAOHToUKW2a2VlBbVaXWG91q1b49ChQxg3bpy27Mcff0Tr1q2rdBxr1qyBlZUVxowZAwD4v//7P0yaNAl37tzRafkg88NuAiIz9sEHHyAoKAh2dnYGm2+uFH369EH37t0xdOhQfPfdd7h69Sp+/PFHvP322zh+/DgAwN/fHwkJCTh9+jR++uknjBkzRm9QX2UUFRVh4sSJ+PXXX7Fnzx7MmTMH06ZNg4WFBVxcXODm5obPPvsMFy9exP79+/Hqq69Waru+vr64ffs29u3bh+zsbNy5c6fUejNnzkRcXBxWrlyJ33//HR9//DESEhLKHHBYGiEEPvnkE0ybNg1WVlYAgCFDhsDCwgLbtm2r9HbINDEZIDJjRUVFGDlyJKZMmSJ3KGZHkiTs3r0bPXv2xIQJE/DYY49h1KhRuHr1Kho2bAhAcyEeFxcXBAUFYdCgQQgPD0enTp2qvK/evXujRYsW6NmzJ5555hkMGjRIOw3QwsIC8fHxOHHiBNq2bYtXXnkFCxYsqNR2g4KCMHnyZDz77LNo0KAB5s+fX2q9oUOHYsmSJViwYAHatGmDVatWITY2FiEhIZU+hu+++w5paWmYPHmytszW1hZPP/00ZxXUAZIQQsgdBBHVTFxcHKZPn46///67ys/Ny8uDs7MzcnNz4eTkZPjgFC4yMhJ///03fz2TSeOYASKFKSwsRGFhofZ+Xl6ejNEQkSlgNwGRwsTExMDZ2Vl78/HxkTskIpIZkwEiEzN37lzt5WXLupUMcKuO2bNnIzc3V3srmetOxhEXF8cuAjJ57CYgMjHTpk3DqFGjyq1T1lS4yrC2toa1tXW1n09EdQ+TASIT4+7uDnd3d7nDICIFYTJAZMbS0tKQk5ODtLQ0qNVq7QI8/v7+cHBwkDc4IjIbnFpIZMYiIyOxZs0avfLk5ORKzyHn1EIiYjJApHBMBoiI3QRmSq0GUlOBrCzA0xMIDgZUKrmjIiIic8RkwAwlJADR0cC1a/+UeXsDS5YAw4fLFxcREZknXmfAzCQkACNG6CYCAJCRoSlPSJAnLiIiMl8cM2BG1GrA11c/ESghSZoWgitX2GVAlccxA0TElgEzkppadiIAAEIA6emaekRERJXFZMCMZGUZth4RERHAZMCseHoath4RERHAZMCsBAdrxgRIUumPSxLg46OpR0REVFlMBsyISqWZPgjoJwQl9xcv5uBBIiKqGiYDZmb4cOCbb4DGjXXLvb015bzOABERVRWnFpopXoGQDIVTC4mIVyA0UyoVUMl1aIiIiMrFbgIiIiKFYzJARESkcEwGiIiIFI7JABERkcIxGSAiIlI4ziYgIh2ctkqkPEwGiEgrIQGIjtZdHdPbW3PlS17QiqjuYjcBEQEAduwARozQXyY7I0NTnpAgT1xEZHy8AiGRwpVcgdDLKxeZmaVfgVCSNC0EV66wy4CoLmLLABEBADIzy35MCCA9XTOWgIjqHiYDRFRpWVlyR0BExsBkgIgqzdNT7giIyBiYDBARAMDLSzM2oDSSBPj4aKYZElHdw2SAiAAAH32k+ffRhKDk/uLFHDxIVFcxGSAiAMDgwcA33wCNG+uWe3trynmdAaK6i1MLiRSuZGphbm4unJyceAVCIgXiFQiJSIdKBYSEyB0FEdUmdhMQEREpHFsGiIiozmA3V/WwZYDITF29ehUTJ05Es2bNYGtrCz8/P8yZMwdFRUVyh0Yki4QEwNcXCA0FxozR/Ovry3U1KoMtA0Rm6vz583jw4AFWrVoFf39//PLLL5g0aRIKCgqwcOFCucMjqlUJCZoFtR4dEl+y0BZnxJSPswmI6pAFCxZgxYoVuHz5cqWf8+hsAiJzo1ZrWgAeXXGzBBfaqhi7CYjqkNzcXLi6upZbp7CwEHl5eTo3InOWmlp2IgBwoa3KYDJAVEdcunQJS5cuxeTJk8utFxMTA2dnZ+3Nx8enliIkMo7KLqDFhbbKxmSAyMTMnTsXkiSVezt+/LjOczIzM9GvXz+MHDkSzz//fLnbnz17NnJzc7W39PR0Yx4OkdFVdgEtLrRVNo4ZIDIx2dnZyM7OLreOr68vbGxsAGgSgdDQUHTr1g1xcXGwsKhajs8xA2TuSsYMZGToDyAEOGagMjibgMjEuLu7w93dvVJ1MzIyEBoais6dOyM2NrbKiQBRXaBSAUuWaGYNSJJuQsCFtiqHfzmIzFRmZiZCQkLg4+ODhQsX4ubNm7h+/TquX78ud2hEtW74cC60VRPsJiAyU3FxcYiKiir1saqc1uwmoLqEVyCsHiYDRArHZICI2E1ARESkcBxASFXCJjiimuE5RKaIyQBVWkICEB2te6Uvb2/NKF4OziGqGM8hMlUcM0CVUtYiICXTdjha13xxzEDt4DlEpozJAFWIi4DUbUwGjI/nEJk6DiCkCnEREKKa4TlEpo7JAFWIi4AQ1QzPITJ1TAaoQlwEhKhmeA6RqeOYAaoQFwGp2zhmwPh4DpGpY8sAVahkERDgn5HPJbgICFHFeA6RqWMyQJXCRUCIaobnEJkydhNQlfDqaXUPuwlqF88hMkVMBogUjskAEbGbgIiISOGYDBARESkckwEiIiKFYzJARESkcEwGiIiIFI7JABERkcIxGSAiIlI4JgNEREQKx2SAiIhI4ZgMEBERKRyTASIiIoVjMkBERKRwTAaIiIgUjskAERGRwjEZICIiUjgmA0RERArHZICIiEjhmAwQEREpHJMBIiIihWMyQGTGBg8ejCZNmsDGxgaenp4YO3YsMjMz5Q6LiMwMkwEiMxYaGorNmzfjwoUL2LJlCy5duoQRI0bIHRYRmRlJCCHkDoKIDGPHjh0YOnQoCgsLYWlpWWqdwsJCFBYWau/n5eXBx8cHubm5cHJyqq1QiciEsGWAqI7IycnBV199haCgoDITAQCIiYmBs7Oz9ubj41OLURKRKWIyQGTmZs2aBXt7e7i5uSEtLQ3bt28vt/7s2bORm5urvaWnp9dSpERkqpgMEJmYuXPnQpKkcm/Hjx/X1p85cyZOnTqFxMREqFQqjBs3DuX1/llbW8PJyUnnRkTKxjEDRCYmOzsb2dnZ5dbx9fWFjY2NXvm1a9fg4+ODH3/8Ed27d6/U/vLy8uDs7MwxA0QKVk/uAIhIl7u7O9zd3av13JLc/uEBgkREFWEyQGSmjh49iqNHj6JHjx5wcXHB5cuX8e6778LPz6/SrQJERADHDBCZLVtbWyQkJKB3795o2bIlJkyYgLZt2+LAgQOwtraWOzwiMiMcM0CkcBwzQERsGSAiIlI4JgNEREQKx2SAiIhI4ZgMEBERKRyTASIiIoVjMkBERKRwTAaIiIgUjskAERGRwjEZICIiUjgmA0RERArHZICIiEjhmAwQEREpHJMBIiIihWMyQEREpHBMBoiIiBSOyQAREZHCMRkgIiJSOCYDRERECsdkgIiISOGYDBARESkckwEiIiKFYzJARESkcEwGiIiIFK6e3AEQEZHpUauB1FQgKwvw9ASCgwGVSu6oyFiYDBARkY6EBCA6Grh27Z8yb29gyRJg+HD54iLjYTcBERFpJSQAI0boJgIAkJGhKU9IkCcuMi5JCCHkDoKI5JOXlwdnZ2fk5ubCyclJ7nBIRmo14OurnwiUkCRNC8GVK+wyqGvYMkBERAA0YwTKSgQAQAggPV1Tj+oWjhkgohrjYLO6ISvLsPXIfDAZIKIa4WCzusPT07D1yHxwzACRwtVkzEDJYLNH/4pIkubfb75hQmBOSsYMZGTov6cAxwzUZRwzQFQHFBYWokOHDpAkCadPn66VfarVmhaB0r40SsqmT9fUI/OgUmladIB/EroSJfcXL2YiUBcxGSCqA15//XV4eXnV6j452KxuGj5c06LTuLFuubc3W3rqMo4ZIDJze/bsQWJiIrZs2YI9e/bU2n452KzuGj4cGDKEg0KVhMkAkRn7888/MWnSJGzbtg12dnaVek5hYSEKCwu19/Py8qq1bw42q9tUKiAkRO4oqLawm4DITAkhEBkZicmTJ6NLly6Vfl5MTAycnZ21Nx8fn2rtPzhY03T8aN9yCUkCfHw09YjItDEZIDIxc+fOhSRJ5d6OHz+OpUuXIi8vD7Nnz67S9mfPno3c3FztLT09vVpxcrAZUd3BqYVEJiY7OxvZ2dnl1vH19cWoUaOwc+dOSA99E6vVaqhUKjz33HNYs2ZNpfZX08sRl3adAR8fTSLAwWZE5oHJAJGZSktL0+nvz8zMRHh4OL755ht069YN3t7eldqOIdYm4BUIicwbBxASmakmTZro3HdwcAAA+Pn5VToRMBQONiMybxwzQEREpHBsGSCqI3x9fcFePyKqDrYMEBERKRyTASIiIoVjMkBERKRwTAaIiIgUjskAERGRwjEZICIiUjgmA0RERArHZICIiEjhmAwQEREpHBcqIlI4IQTy8/Ph6OioswIiESkHkwEiIiKFYzcBERGRwjEZICIiUjgmA0RERArHZICIiEjhmAwQEREpHJMBIiIihWMyQEREpHD/D7yUETx5eYO6AAAAAElFTkSuQmCC", + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "using PyPlot\n", + "plot(real(λ), imag(λ), \"bo\")\n", + "xlabel(lpad(\"real part of λ\",60))\n", + "ylabel(lpad(\"imaginary part of λ\", 60))\n", + "title(\"eigenvalues of a random 20×20 matrix\")\n", + "\n", + "\n", + "function center_spines(ax)\n", + " ax.spines[\"left\"].set_position(\"center\")\n", + " ax.spines[\"right\"].set_color(\"none\")\n", + " ax.spines[\"bottom\"].set_position(\"center\")\n", + " ax.spines[\"top\"].set_color(\"none\")\n", + " #ax.spines[\"left\"].set_smart_bounds(true)\n", + " #ax.spines[\"bottom\"].set_smart_bounds(true)\n", + " ax.xaxis.set_ticks_position(\"bottom\")\n", + " ax.yaxis.set_ticks_position(\"left\")\n", + "end\n", + "center_spines(gca())" + ] + }, + { + "cell_type": "markdown", + "id": "3815694b", + "metadata": {}, + "source": [ + "### Hermitian/real-symmetric matrics and the discrete Laplacian\n", + "\n", + "As another example, let's construct the famous \"discrete Laplacian\" matrix $L$. (You may have seen $-L$ instead, which can be viewed as a finite-difference approximation for $d^2/dx^2$ on a uniform grid. Here, it turns out to be slightly nicer to flip the sign, which flips the sign of all the eigenvalues.)\n", + "\n", + "$$\n", + "L = \\begin{pmatrix} 2 & -1 & & & & \\\\\n", + "-1 & 2 & -1 & & & \\\\\n", + "& \\ddots & \\ddots & \\ddots & & \\\\\n", + "& & -1 & 2 & -1 & \\\\\n", + "& & & -1 & 2 & -1 \\\\\n", + "& & & & -1 & 2 \\end{pmatrix}\n", + "$$\n", + "which can be constructed in Julia in a variety of ways, e.g.:" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "778ceffd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×6 Matrix{Int64}:\n", + " 2 -1 0 0 0 0\n", + " -1 2 -1 0 0 0\n", + " 0 -1 2 -1 0 0\n", + " 0 0 -1 2 -1 0\n", + " 0 0 0 -1 2 -1\n", + " 0 0 0 0 -1 2" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = 6\n", + "L = diagm(-1 => fill(-1,m-1), 0 => fill(2,m), +1 => fill(-1,m-1))" + ] + }, + { + "cell_type": "markdown", + "id": "6329285f", + "metadata": {}, + "source": [ + "This is a **real-symmetric tridiagonal matrix** $L = L^T$. Julia actually has a special type for this sort of matrix, that allows it do computations extra efficiently (it doesn't even store the zero off-diagonal elements):" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "28a228c2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×6 SymTridiagonal{Int64, Vector{Int64}}:\n", + " 2 -1 ⋅ ⋅ ⋅ ⋅\n", + " -1 2 -1 ⋅ ⋅ ⋅\n", + " ⋅ -1 2 -1 ⋅ ⋅\n", + " ⋅ ⋅ -1 2 -1 ⋅\n", + " ⋅ ⋅ ⋅ -1 2 -1\n", + " ⋅ ⋅ ⋅ ⋅ -1 2" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "L = SymTridiagonal(fill(2,m), fill(-1,m-1))" + ] + }, + { + "cell_type": "markdown", + "id": "f436ff44", + "metadata": {}, + "source": [ + "The eigenvalues of *any* real-symmetric matrix are very special: they are **always real**: " + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "2cca5c08", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6-element Vector{Float64}:\n", + " 0.1980622641951617\n", + " 0.7530203962825329\n", + " 1.5549581320873713\n", + " 2.4450418679126287\n", + " 3.2469796037174667\n", + " 3.8019377358048385" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigvals(L)" + ] + }, + { + "cell_type": "markdown", + "id": "9147adfa", + "metadata": {}, + "source": [ + "Moreover, the eigenvectors are **orthogonal** (and can be normalized to **orthonormal**):" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "0816eb9d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×6 Matrix{Float64}:\n", + " 0.231921 -0.417907 0.521121 -0.521121 -0.417907 0.231921\n", + " 0.417907 -0.521121 0.231921 0.231921 0.521121 -0.417907\n", + " 0.521121 -0.231921 -0.417907 0.417907 -0.231921 0.521121\n", + " 0.521121 0.231921 -0.417907 -0.417907 -0.231921 -0.521121\n", + " 0.417907 0.521121 0.231921 -0.231921 0.521121 0.417907\n", + " 0.231921 0.417907 0.521121 0.521121 -0.417907 -0.231921" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Q = eigvecs(L)" + ] + }, + { + "cell_type": "markdown", + "id": "ee9f8087", + "metadata": {}, + "source": [ + "Let's check that they are orthogonal: $Q^T Q \\approx I$:" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "ae413c92", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×6 Matrix{Float64}:\n", + " 1.0 3.2e-16 -3.23485e-16 … -6.06815e-17 -1.26907e-16\n", + " 3.2e-16 1.0 1.63628e-16 1.09586e-16 8.61728e-17\n", + " -3.23485e-16 1.63628e-16 1.0 -4.28431e-17 -4.34141e-17\n", + " -1.92391e-16 -2.27856e-17 -1.05864e-16 -1.0623e-16 -7.7364e-17\n", + " -6.06815e-17 1.09586e-16 -4.28431e-17 1.0 -2.34469e-16\n", + " -1.26907e-16 8.61728e-17 -4.34141e-17 … -2.34469e-16 1.0" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Q' * Q" + ] + }, + { + "cell_type": "markdown", + "id": "69e47f3d", + "metadata": {}, + "source": [ + "It's a little hard to see because of the roundoff errors. Let's round to 3 decimal places:" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "8673e9fd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×6 Matrix{Float64}:\n", + " 1.0 0.0 -0.0 -0.0 -0.0 -0.0\n", + " 0.0 1.0 0.0 -0.0 0.0 0.0\n", + " -0.0 0.0 1.0 -0.0 -0.0 -0.0\n", + " -0.0 -0.0 -0.0 1.0 -0.0 -0.0\n", + " -0.0 0.0 -0.0 -0.0 1.0 -0.0\n", + " -0.0 0.0 -0.0 -0.0 -0.0 1.0" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "round.(Q' * Q, digits=3)" + ] + }, + { + "cell_type": "markdown", + "id": "d0350e16", + "metadata": {}, + "source": [ + "Let's try it for a random real-symmetric matrix:\n", + "\n", + "$$\n", + "B = R + R^T\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "20d048d5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "20×20 Matrix{Float64}:\n", + " 1.65834 2.33138 1.49424 … -0.22269 -0.261395 -0.548495\n", + " 2.33138 1.49698 0.567878 1.17445 1.58359 0.841074\n", + " 1.49424 0.567878 -1.82936 -2.08445 -1.97305 0.802223\n", + " -0.873223 -2.5189 -0.234729 -0.466553 1.98014 -0.757998\n", + " 1.19066 -1.5483 1.26955 -3.18594 1.4984 -0.443326\n", + " -1.28653 -0.70039 0.740466 … 1.43708 3.38325 -0.679315\n", + " 1.70158 -0.598096 -0.780981 1.26355 0.222589 -0.101874\n", + " 1.0822 1.0677 -0.795777 -1.73173 -1.37277 0.291422\n", + " 0.100375 1.97324 -0.942124 -1.77077 -1.75298 0.469773\n", + " -2.21561 -0.0290779 0.133419 -1.01377 1.71847 -1.6003\n", + " 0.950668 3.11415 -0.0590488 … -0.0276943 -0.692113 0.181601\n", + " -0.0765708 1.70557 0.651119 0.942239 -0.81868 0.0588557\n", + " 1.5453 -0.70579 0.0573533 -1.9466 1.0038 1.38634\n", + " 0.963923 0.535501 0.715344 1.47753 -1.9614 -0.0941366\n", + " 0.23409 1.8659 -2.18804 0.570987 -1.32442 1.34531\n", + " -1.83025 -0.434536 -1.73489 … -2.13741 0.532592 -1.74014\n", + " 2.5416 -0.0067432 -0.657476 -1.14267 -1.57061 0.745832\n", + " -0.22269 1.17445 -2.08445 -1.5944 1.62804 -0.285623\n", + " -0.261395 1.58359 -1.97305 1.62804 1.83709 0.731884\n", + " -0.548495 0.841074 0.802223 -0.285623 0.731884 0.986415" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "B = R + R'" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "aae3554e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "20-element Vector{Float64}:\n", + " -9.885815419012882\n", + " -9.084302010408306\n", + " -6.843285604985337\n", + " -5.942971471336738\n", + " -4.7145012740161105\n", + " -3.6586794713705375\n", + " -2.5987759869356073\n", + " -2.125904297844215\n", + " -0.8883340910207597\n", + " 0.5950921105698613\n", + " 1.4203119332808172\n", + " 1.8872687001818484\n", + " 2.3827569043089762\n", + " 3.2628973805220824\n", + " 4.577517156826758\n", + " 6.138927327909758\n", + " 7.186237946968724\n", + " 7.649004117817263\n", + " 8.873524934813709\n", + " 12.281154377914909" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigvals(B)" + ] + }, + { + "cell_type": "markdown", + "id": "031d6f94", + "metadata": {}, + "source": [ + "They are all real numbers, but they are not all positive!\n", + "\n", + "$L$ must be special in *another* way!" + ] + }, + { + "cell_type": "markdown", + "id": "65d950e0", + "metadata": {}, + "source": [ + "### Positive-definite matrices and the discrete Laplacian\n", + "\n", + "In fact, not only were the eigenvalues real, but they were actually **all positive**.\n", + "\n", + "A real-symmetric (or Hermitian) matrix with *positive* real eigenvalues is called **positive definite**.\n", + "\n", + "(The matrix $-L$ has all-negative eigenvalues, which is called **negative definite**.)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "b227a1eb", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "true" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "isposdef(L) # check positive-definiteness" + ] + }, + { + "cell_type": "markdown", + "id": "383da99f", + "metadata": {}, + "source": [ + "Positive-definiteness of $L$ is no accident, because it can be factored as a product of two \"difference matrices\" in a *very special* form:\n", + "\n", + "$$\n", + "L = \\underbrace{\\begin{pmatrix} 2 & -1 & & & \\\\\n", + "-1 & 2 & -1 & & \\\\\n", + "& \\ddots & \\ddots & \\ddots & \\\\\n", + "& & -1 & 2 & -1 \\\\\n", + "& & & -1 & 2 \\end{pmatrix}}_{m \\times m}\n", + "= D^T D = \\\\\n", + "-\n", + "\\underbrace{\\begin{pmatrix} -1 & 1 & & & \\\\\n", + " & -1 & 1 & & \\\\\n", + "& & \\ddots & \\ddots & & \\\\\n", + "& & & -1 & 1 & \\\\\n", + "& & & & -1 & 1 \\end{pmatrix}}_{m \\times (m+1) \\mbox{ matrix } -D^T}\n", + "\\underbrace{\\begin{pmatrix} 1 & & & & & \\\\\n", + "-1 & 1 & & & & \\\\\n", + "& \\ddots & \\ddots & & & \\\\\n", + "& & -1 & 1 & & \\\\\n", + "& & & -1 & 1 \\\\\n", + "& & & & -1 \\end{pmatrix}}_{(m+1) \\times m \\mbox{ matrix } D}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "89bf7016", + "metadata": {}, + "source": [ + "Let's check:" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "9d4b0a42", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "7×6 Matrix{Int64}:\n", + " 1 0 0 0 0 0\n", + " -1 1 0 0 0 0\n", + " 0 -1 1 0 0 0\n", + " 0 0 -1 1 0 0\n", + " 0 0 0 -1 1 0\n", + " 0 0 0 0 -1 1\n", + " 0 0 0 0 0 -1" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "D = diagm(m+1, m, -1 => fill(-1, m), 0 => fill(1, m))" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "acb7f38e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×6 Matrix{Int64}:\n", + " 2 -1 0 0 0 0\n", + " -1 2 -1 0 0 0\n", + " 0 -1 2 -1 0 0\n", + " 0 0 -1 2 -1 0\n", + " 0 0 0 -1 2 -1\n", + " 0 0 0 0 -1 2" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "D'D" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "51eafbaa", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "true" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "D'D == L" + ] + }, + { + "cell_type": "markdown", + "id": "05ecfaf0", + "metadata": {}, + "source": [ + "It is easy to see that *any* matrix $A = B^T B$ for *any* real $B$ must have eigenvalues $\n", + "\\lambda \\ge 0$.\n", + "\n", + "It is clearly real-symmetric so its eigenvalues and eigenvalues are real. If $Ax = \\lambda x$, then\n", + "$$\n", + "x^T A x = x^T B^T B x = (Bx)^T (Bx) = \\Vert Bx \\Vert^2 = \\lambda x^T x = \\lambda \\Vert x \\Vert^2\n", + "$$\n", + "Hence\n", + "$$\n", + "\\lambda = \\frac{x^T A x}{x^T x} = \\frac{\\Vert Bx \\Vert^2}{ \\Vert x \\Vert^2} \\ge 0\n", + "$$\n", + "which is obviously $\\ge 0$. Any such matrix $A$ is **positive semidefinite**.\n", + "\n", + "For it to be **positive definite**, i.e. $\\lambda > 0$ (not just $\\ge 0$), we must have $Bx \\ne 0$ for any $x \\ne 0$: **B must be full column rank (independent columns)**.\n", + "\n", + "In the case above, it is obvious that $D$ has rank $m$ ($m$ independent columns). For example, you can see this because $D^T$ is upper-triangular with $m$ nonzero pivots on the diagonal." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96b6d134", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "@webio": { + "lastCommId": null, + "lastKernelId": null + }, + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}