{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Arnoldi vs Power Iteration"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "\n",
        "import matplotlib.pyplot as pt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let us make a matrix with a defined set of eigenvalues and eigenvectors, given by `eigvals` and `eigvecs`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16. 17. 18.\n",
            " 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36.\n",
            " 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50.]\n",
            "[50.  1. 49.  2. 48. 47.  3. 46. 45.  4.  5. 44.  6. 43.  7.  8.  9. 42.\n",
            " 41. 40. 10. 11. 12. 13. 14. 39. 15. 38. 37. 36. 16. 35. 34. 17. 33. 32.\n",
            " 18. 31. 30. 29. 19. 28. 20. 21. 27. 22. 26. 25. 24. 23.]\n"
          ]
        }
      ],
      "source": [
        "np.random.seed(40)\n",
        "\n",
        "# Generate matrix with eigenvalues 1...50\n",
        "n = 50\n",
        "eigvals = np.linspace(1., n, n)\n",
        "eigvecs = np.random.randn(n, n)\n",
        "#To work with symmetric matrix, orthogonalize eigvecs\n",
        "eigvecs, R = la.qr(eigvecs)\n",
        "\n",
        "print(eigvals)\n",
        "\n",
        "A = la.solve(eigvecs, np.dot(np.diag(eigvals), eigvecs))\n",
        "print(la.eig(A)[0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Initialization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set up $Q$ and $H$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [],
      "source": [
        "Q = np.zeros((n, n))\n",
        "H = np.zeros((n, n))\n",
        "\n",
        "k = 0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Pick a starting vector, normalize it"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [],
      "source": [
        "x0 = np.random.randn(n)\n",
        "x0 = x0/la.norm(x0)\n",
        "\n",
        "# Set the first column of Q to be the normalized starting vector\n",
        "Q[:, k] = x0.copy()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a list to save arrays of Ritz values:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [],
      "source": [
        "ritz_values = []\n",
        "ritz_max = []"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Algorithm"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Carry out one iteration of Arnoldi iteration.\n",
        "\n",
        "Run this cell in-place (Ctrl-Enter) until H is filled."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [],
      "source": [
        "def Arnoldi_step(A,Q,H,k):\n",
        "    \n",
        "    u = A @ Q[:, k]\n",
        "\n",
        "    # Carry out Gram-Schmidt on u against Q\n",
        "    # to do Lanczos change range start to k-1\n",
        "    for j in range(0,k+1):\n",
        "        qj = Q[:, j]\n",
        "        H[j,k] = qj @ u\n",
        "        u = u - H[j,k]*qj\n",
        "\n",
        "    if k+1 < n:\n",
        "        H[k+1, k] = la.norm(u)\n",
        "        Q[:, k+1] = u/H[k+1, k]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "14\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD8CAYAAACxd9IeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAALJklEQVR4nO3dUahkhX3H8e+vu1oDoajJZVlcqRalYR8aZS9WsQ/FVGpNiD5ISQhhHxb2JQFDA6m2UBooNL7E5KEvS5TsQ4imJqBIoNjNhlAo2ptoUnVpdyOGKqt7JUqSl7Sb/Pswx3S93u2dnTszd2b/3w8M95wzc/f8Xe7XM+fcmdlUFZIufr+10wNImg9jl5owdqkJY5eaMHapCWOXmphr7EnuSPIfSU4luW+e+x5HkoeTnEny/DnbrkzyVJKTw9crdnLGtyW5OsnxJC8meSHJvcP2RZ33siTPJPnhMO/nh+3XJnl6+Jl4NMmlOz3r25LsSvJskieH9YWddRxziz3JLuAfgD8D9gMfT7J/Xvsf01eBOzZsuw84VlXXA8eG9UVwFvhsVe0HbgY+Nfx9Luq8vwRuq6oPAjcAdyS5GXgAeLCqrgPeBA7t3Ijvci9w4pz1RZ51S/M8st8EnKqql6rqv4FHgLvmuP8tVdX3gJ9u2HwXcHRYPgrcPc+ZzqeqTlfVD4blnzP6obyKxZ23quoXw+olw62A24DHhu0LM2+SfcCHga8M62FBZx3XPGO/Cvivc9ZfGbYtuj1VdXpYfg3Ys5PDbCbJNcCNwNMs8LzD0+LngDPAU8CPgbeq6uzwkEX6mfgS8Dng18P6+1jcWcfiBboLUKPXFi/U64uTvBf4JvCZqvrZufct2rxV9auqugHYx+iZ3gd2dqLNJfkIcKaqvr/Ts0zT7jnu61Xg6nPW9w3bFt3rSfZW1ekkexkdlRZCkksYhf61qvrWsHlh531bVb2V5DhwC3B5kt3DEXNRfiZuBT6a5E7gMuB3gC+zmLOObZ5H9n8Drh+uaF4KfAx4Yo77n9QTwMFh+SDw+A7O8hvDOeRDwImq+uI5dy3qvCtJLh+W3wPczug6w3HgnuFhCzFvVd1fVfuq6hpGP6ffqapPsICzXpCqmtsNuBP4T0bnan89z32POd/XgdPA/zA6JzvE6FztGHAS+Gfgyp2ec5j1jxg9Rf8R8Nxwu3OB5/0D4Nlh3ueBvxm2/x7wDHAK+Efgt3d61g1z/zHw5DLMutUtw3+EpIucF+ikJoxdasLYpSaMXWrC2KUmdiT2JId3Yr+TWKZZYbnmXaZZYfnm3WhbsW/jLavL9Je2TLPCcs27TLPC8s37DhPHviRvWZU0mPhFNUluAf62qv50WL8foKr+/v/5nvav4Dlw4MBM//z19XVWVlZmuo9pWaZZYTnmffnll3njjTey2X3beSPMZm9Z/cNt/HktrK2t7fQIuoitrq6e976Zv+ttuKix1Oc60sVgO7GP9ZbVqjoCHAGfxks7aTtX45f1LatSSxMf2avqbJJPA/8E7AIerqoXpjaZpKna1jl7VX0b+PaUZpE0Q75cVmpirrEfOHBg46eASJoTj+xSE8YuNWHsUhPz/Nz4d9nsvH30CcmSps0ju9SEsUtNGLvUxI6es29m43m85/DSdHhkl5owdqkJY5eaMHapiYW7QLeRF+yk6fDILjVh7FITxi41sfDn7Bv55hlpMh7ZpSaMXWrC2KUmjF1qYuku0G3GF95IW/PILjVh7FITxi41cVGcs2/kC2+kd/PILjVh7FITxi41cVGes2/G38WrO4/sUhPGLjVh7FITW8ae5OEkZ5I8f862K5M8leTk8PWK2Y4pabvGObJ/Fbhjw7b7gGNVdT1wbFhfKlX1jpt0sdsy9qr6HvDTDZvvAo4Oy0eBu6c7lqRpm/ScfU9VnR6WXwP2TGkeSTOy7Qt0NXoOfN7nwUkOJ1lLsra+vr7d3Uma0KSxv55kL8Dw9cz5HlhVR6pqtapWV1ZWJtzd7G08h/c8XhebSWN/Ajg4LB8EHp/OOJJmZZxfvX0d+Ffg95O8kuQQ8AXg9iQngT8Z1iUtsC1fG19VHz/PXR+a8iySZqjNG2Em4ZtndDHx5bJSE8YuNWHsUhPGLjXhBboL4AU7LTOP7FITxi41YexSE56zb4P/8oyWiUd2qQljl5owdqkJz9mnzN/Fa1F5ZJeaMHapCWOXmjB2qQkv0M2Yn1KrReGRXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qYsvYk1yd5HiSF5O8kOTeYfuVSZ5KcnL4esXsx5U0qXGO7GeBz1bVfuBm4FNJ9gP3Aceq6nrg2LAuaUFtGXtVna6qHwzLPwdOAFcBdwFHh4cdBe6e0YySpuCCztmTXAPcCDwN7Kmq08NdrwF7pjuapGkaO/Yk7wW+CXymqn527n01+qC1TT9sLcnhJGtJ1tbX17c1rKTJjRV7kksYhf61qvrWsPn1JHuH+/cCZzb73qo6UlWrVbW6srIyjZklTWCcq/EBHgJOVNUXz7nrCeDgsHwQeHz640malnE+SvpW4JPAvyd5btj2V8AXgG8kOQT8BPjzmUwoaSq2jL2q/gU4379O+KHpjiNpVnwFndSEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71ISxS00Yu9SEsUtNGLvUhLFLTRi71MSWsSe5LMkzSX6Y5IUknx+2X5vk6SSnkjya5NLZjytpUuMc2X8J3FZVHwRuAO5IcjPwAPBgVV0HvAkcmtmUkrZty9hr5BfD6iXDrYDbgMeG7UeBu2cxoKTpGOucPcmuJM8BZ4CngB8Db1XV2eEhrwBXned7DydZS7K2vr4+hZElTWKs2KvqV1V1A7APuAn4wLg7qKojVbVaVasrKyuTTSlp2y7oanxVvQUcB24BLk+ye7hrH/DqdEeTNE3jXI1fSXL5sPwe4HbgBKPo7xkedhB4fEYzSpqC3Vs/hL3A0SS7GP3P4RtV9WSSF4FHkvwd8Czw0AznlLRNW8ZeVT8Cbtxk+0uMzt8lLQFfQSc1YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhPGLjVh7FITxi41YexSE8YuNWHsUhNjx55kV5Jnkzw5rF+b5Okkp5I8muTS2Y0pabsu5Mh+L3DinPUHgAer6jrgTeDQNAeTNF1jxZ5kH/Bh4CvDeoDbgMeGhxwF7p7BfJKmZNwj+5eAzwG/HtbfB7xVVWeH9VeAq6Y7mqRp2jL2JB8BzlTV9yfZQZLDSdaSrK2vr0/yR0iagnGO7LcCH03yMvAIo6fvXwYuT7J7eMw+4NXNvrmqjlTValWtrqysTGFkSZPYMvaqur+q9lXVNcDHgO9U1SeA48A9w8MOAo/PbEpJ27ad37P/JfAXSU4xOod/aDojSZqF3Vs/5P9U1XeB7w7LLwE3TX8kSbPgK+ikJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qwtilJoxdasLYpSaMXWrC2KUmjF1qIlU1v50l68BPgPcDb8xtx9uzTLPCcs27TLPCcsz7u1W1stkdc439NztN1qpqde47nsAyzQrLNe8yzQrLN+9GPo2XmjB2qYmdiv3IDu13Ess0KyzXvMs0KyzfvO+wI+fskubPp/FSE8YuNWHsUhPGLjVh7FIT/wuTf1QOsXqiLAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "print(k)\n",
        "\n",
        "Arnoldi_step(A,Q,H,k)\n",
        "\n",
        "k += 1\n",
        "\n",
        "pt.spy(H)\n",
        "\n",
        "\n",
        "if k>1:\n",
        "    D = la.eig(H)[0]\n",
        "    max_ritz = D[np.argmax(np.abs(D))]\n",
        "    ritz_vals = np.zeros(k)\n",
        "    for i in range(k):\n",
        "        ritz_vals[i] = D[np.argmax(np.abs(D))]\n",
        "        D[np.argmax(np.abs(D))] = 0\n",
        "    ritz_max.append(max_ritz)\n",
        "    ritz_values.append(ritz_vals)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that $Q^T A Q =H$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "7.336036368452853e-15"
            ]
          },
          "execution_count": 22,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(Q[:,:k-1].T @ A @ Q[:,:k-1] - H[:k-1,:k-1])/ la.norm(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that $AQ-QH$ is fairly small"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.06030251657738512"
            ]
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(A @ Q[:,:k-1] - Q[:,:k-1]@H[:k-1,:k-1])/ la.norm(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that Q is orthogonal:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "1.349310944636171e-13"
            ]
          },
          "execution_count": 24,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm((Q.T.conj() @ Q)[:k-1,:k-1] - np.eye(k-1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compare Max Ritz Value to Power Iteration"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "48.04487764639638 49.98334764109434\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f8d39556ee0>"
            ]
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAerElEQVR4nO3dfXRU9b3v8ffXyIHUglhMXZDoJbYlPvE8ggqJSLVYQUoVjnQplUIPl/b4gHpAWJ51pXddb63pKmivV5dXCtSD1IoIrT0tbcUU8LEJQVExHCuxJVBNqWCpoSThe/+YScjDQGbI7JnZzOe1Fmtm/7L3nm94+PDLb//2b5u7IyIi4XNKpgsQEZETowAXEQkpBbiISEgpwEVEQkoBLiISUqem88POPPNMHzhwYDo/UkQk9Kqqqv7i7gUd29Ma4AMHDqSysjKdHykiEnpm9n68dg2hiIiElAJcRCSkFOAiIiGlABcRCSkFuIhISCU0C8XMaoG/Ac1Ak7tHzOwzwFPAQKAW+Gd3/yiYMkVEgrGuuo7yDTXs2d/AgL75zJ9QwpThhVl/bkhuGuEV7v6XNtsLgefd/X4zWxjbvjtllYlIaAQdVEGdf111HYvWbqehsRmAuv0NLFq7HaDb5w/y3C26M4TyFWBl7P1KYEq3qxGRwKyrrmPM/RspXvgLxty/kXXVdSk776K126nb34BzNKjCcP7yDTWtAduiobGZ8g01WX3uFokGuAO/NrMqM5sTazvL3ffG3v8ZOCvegWY2x8wqzayyvr6+m+WKnLyCCtiWc4cxBIM+/579DUm1Z8u5WyQa4GPdfQTwZeBfzays7Rc9+lSIuE+GcPfH3D3i7pGCgk53goqESlh7sWENwaDPP6BvflLt2XLuFgkFuLvXxV4/BJ4FRgEfmFl/gNjrhymrSiQLhbkXG9YQDPr88yeUkN8jr11bfo885k8oyepzt+gywM3sNDPr3fIe+BLwJvAz4ObYbjcD61NWlUgWCnMvNqwhGPT5pwwv5LvXDaawbz4GFPbN57vXDU7JRcYgz90ikVkoZwHPmlnL/k+6+6/M7PfAT81sNvA+8M8pq0qkG4KasRB0L7YuznlS1YudP6Gk3YwISG0IAoHNQknH+VMZquk6NyQQ4O7+HjA0Tvs+4ItBFCVyooKcuhVkyAYZsBDuEEzH+cMqrcvJigTteMMc3Q2AMPdiWz5DIXhyUYDLSSXIYY6w92Ll5KMAl5NK0GPJClnJJlrMSk4q6Zi6JZIt1AOXtAty3Yx0jCWLZAsFuKRVOhb40TCH5AoNoUhapWOBH5FcoQCXtErHAj8iuUIBLmmVjgV+RHKFAlzSSrNERFJHFzElrTRLRCR1FOCSdpolIpIaGkIREQkpBbiISEgpwEVEQkoBLiISUrqIKXEFuV6JiKSGAlw6Scd6JSLSfRpCkU60XolIOCjApROtVyISDgpw6UTrlYiEgwJcOtF6JSLhoIuY0onWKxEJBwW4xKX1SkS6actSKBwBxWVH23ZtgrqtMHZeSj4i4SEUM8szs2ozey62/UUz22pm28xsi5l9PiUViYicDApHwNMzo6EN0denZ0bbUySZMfDbgR1tth8BbnT3YcCTwL+nrCoRkbArLoNpK6KhvfG+6Ou0Fe175N2UUICbWREwEXi8TbMDfWLvTwf2pKwqEZGTQXEZRGbDpgeirykMb0h8DHwpsADo3abtm8B/mlkD8DFwSbwDzWwOMAfgnHPOOeFCRURCZ9cmqFwGZQuir8Wl6e2Bm9kk4EN3r+rwpTuAa9y9CFgO/CDe8e7+mLtH3D1SUFDQ7YJFRFJmy9KjY9Qtdm2KtndXy5j3tBUw/p6jwykdP68bEhlCGQNMNrNa4CfAeDP7BTDU3V+N7fMUcFnKqhIRSYcgLzTWbW0/5t0yJl63tfvnjjF3T3xns3HAvwFTgD8Dl7n7TjObTbQ3fv3xjo9EIl5ZWXnCxYqIpFxLaEdmR4c5UnyhMRXMrMrdIx3bT2geuLs3mdm/AM+Y2RHgI2BWN2sUEUm/thcayxZkXXgfT1K30rt7hbtPir1/1t0Hu/tQdx/n7u8FU6KI5LQgx6lbztX2QmMKx6iDprVQRCS7BTlOnYYLjUFSgItIdgvyhpg0XGgMktZCEZHsF9Q4dbw1SYrLQjMOrh64iGS/EI9TB0kBLiLdk46LjCEepw6SAlxEuifoVfdCPk4dpKRu5Oku3ciTOuuq6/TABckeIbgZJsxSeiOPZNa66joWrd3e+uT4uv0NLFq7HUAhLpkR4pthwkxDKCFUvqGmNbxbNDQ2U76hJkMVSc7TRcaMUICH0J79DUm1iwRKFxkzRgEeQgP65ifVLhLoTBFdZMwYBXgIzZ9QQn6PvHZt+T3ymD+hJEMVSdYLcqbI2Hmdx7yLy1L24F45Nl3EDKGWC5WahSIJa3s7umaKnDQU4CE1ZXihAluSo5kiJx0NoYjkCs0UOekowEVygWaKnJQU4CLZIOj1RDRT5KSkABfJBkGvJ6KZIiclXcQUyQaaJSInQD1wkWzRdpZIZLbCW7qkABfJFpolIklSgItkA80SkROgABdJlNYTkSyjABdJlNYTkSyTcICbWZ6ZVZvZc7FtM7P7zGynme0ws9uCK1MkC7SdKbLxvqNDHrrYKBmSzDTC24EdQJ/Y9kzgbOA8dz9iZp9NcW0i2UfriUgWSagHbmZFwETg8TbN3wL+p7sfAXD3D1NfnkiW0UwRySKJDqEsBRYAR9q0fQ64wcwqzeyXZvaFeAea2ZzYPpX19fXdq1YkkzRTRLJMlwFuZpOAD929qsOXegKHYk9K/n/Aj+Id7+6PuXvE3SMFBQXdLlgkYzRTRLKMufvxdzD7LjADaAJ6ER0DXwtEgC+7+y4zM2C/u59+vHNFIhGvrKxMSeEicW1ZGp0V0nZsetemaMhqRoeElJlVxTrL7XTZA3f3Re5e5O4DgenARne/CVgHXBHb7XJgZ+rKFTlBQS8KJZJFurOY1f3AKjO7AzgIfDM1JYl0gxaFkhySVIC7ewVQEXu/n+jMFDmGddV1em5lJmiqn+QI3YkZkHXVdSxau526/Q04ULe/gUVrt7Ouui7TpZ38NNVPcoQCPCDlG2poaGxu19bQ2Ez5hpoMVZQjNNVPcogCPCB79jck1S4poql+kkMU4AEZ0Dc/qfacEfSzH7UolOQQBXhA5k8oIb9HXru2/B55zJ9QkqGKsoSm+YmkjJ6JGZCW2SaahdKBpvmJpIwCPEBThhcqsOPRND+RlNAQiqSfpvmJpIQCXNJL0/xEUkYBLumlaX4iKaMxcEmveNP5iss0Di5yAtQDFxEJKQW4dBb0zTYikhIKcOlMN9uIhILGwKUz3WwjEgrqgUt8bW+2icxWeItkIQW4xKebbUSyngJcOtPNNiKhoACXznSzjUgo6CKmdKabbURCQT1wEZGQUoCLiISUAlxEJKQU4CIiIZVwgJtZnplVm9lzHdofMrODqS9NjklrlYgIyfXAbwd2tG0wswhwRkorkq5prRIRIcEAN7MiYCLweJu2PKAcWBBMaXJMbdcq2Xjf0ZtuNM1PJKck2gNfSjSoj7RpuwX4mbvvPd6BZjbHzCrNrLK+vv7EqpTOtFaJSM7rMsDNbBLwobtXtWkbAEwDftjV8e7+mLtH3D1SUFDQrWKlDa1VIpLzErkTcwww2cyuAXoBfYC3gH8A75oZwKfM7F13/3xglcpRbdcqKS6D4lINo4jkoC4D3N0XAYsAzGwc8G/uPqntPmZ2MIzhva66jvINNezZ38CAvvnMn1DClOGFmS6ra8dbq0QBLpIzcnYtlHXVdSxau52GxmYA6vY3sGjtdoDsD3GtVSIiJHkjj7tXdOx9x9o/nbqS0qN8Q01reLdoaGymfENNhioSEUlOzt6JuWd/Q1LtIiLZJmcDfEDf/KTaRUSyTc4G+PwJJeT3yGvXlt8jj/kTSjJUkYhIcnL2ImbLhcpQzkIRESGHAxyiIa7AFpGwytkhFBGRsFOAi4iElAI8KFqzW0QCpgAPitbsFpGA5fRFzEC1XbM7Mju6YqAWmxKRFFIPPEhas1tEAqQAD5LW7BaRACnAg9J2ze7x9xwdTlGIi0iKKMCDcrw1u0VEUkAXMYOiNbtFJGDqgYuIhJQCXEQkpBTgIiIhpQAXEQkpBbiISEgpwEVEQkoBLiISUgpwEZGQUoCLiIRUwgFuZnlmVm1mz8W2V5lZjZm9aWY/MrMewZUpIiIdJdMDvx3Y0WZ7FXAeMBjIB76ZwrpERKQLCQW4mRUBE4HHW9rc/T89BngNKAqmRBERiSfRHvhSYAFwpOMXYkMnM4BfxTvQzOaYWaWZVdbX159onSIi0kGXAW5mk4AP3b3qGLv8X2CTu2+O90V3f8zdI+4eKSgo6EapIiLSViLLyY4BJpvZNUAvoI+Z/Ye732Rm9wIFwH8PskgREemsyx64uy9y9yJ3HwhMBzbGwvubwATga+7eaWhFRESC1Z154I8CZwEvm9k2M/sfKapJREQSkNQTedy9AqiIvdfTfEREMkh3YoqIhFTuBviWpZ2fEL9rU7RdRCQEcjfAC0fA0zOPhviuTdHtwhGZrEpEJGG5O45dXAbTVkRDOzIbKpdFt/XUeBEJidztgUM0rCOzYdMD0VeFt4iESG4H+K5N0Z532YLoa8cxcRGRLJa7Ad4y5j1tBYy/5+hwikJcREIidwO8bmv7Me+WMfG6rZmsSkQkYbl7EXPsvM5txWUaBxeR0MjdHriISMgpwEVEQkoBLiISUgpwEZGQUoCLiISUAlxEJKQU4CIiIaUAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkFKAi4iElAJcRCSkFOAiIiGVcICbWZ6ZVZvZc7HtYjN71czeNbOnzOyfgitTREQ6SqYHfjuwo83294Al7v554CNgdioLExGR40sowM2sCJgIPB7bNmA8sCa2y0pgSgD1sa66jjH3b6R44S8Yc/9G1lXXBfExIiKhk+gj1ZYCC4Dese1+wH53b4pt7wYK4x1oZnOAOQDnnHNOUsWtq65j0drtNDQ2A1C3v4FFa7cDMGV43I8TEckZXfbAzWwS8KG7V53IB7j7Y+4ecfdIQUFBUseWb6hpDe8WDY3NlG+oOZFSREROKon0wMcAk83sGqAX0Ad4EOhrZqfGeuFFQMrHNvbsb0iqXUQkl3TZA3f3Re5e5O4DgenARne/EXgBmBrb7WZgfaqLG9A3P6l2EZFc0p154HcDd5rZu0THxJelpqSj5k8oIb9HXru2/B55zJ9QkuqPEhEJnUQvYgLg7hVARez9e8Co1Jd0VMuFyvINNezZ38CAvvnMn1CiC5giMY2NjezevZtDhw5luhRJgV69elFUVESPHj0S2j+pAM+EKcMLFdgix7B792569+7NwIEDic7ulbByd/bt28fu3bspLi5O6BjdSi8SYocOHaJfv34K75OAmdGvX7+kfppSgIuEnML75JHsn6UCXEQkpBTgIjkkqKUp1q1bh5nxzjvvpOR8LcaNG0dlZWWn9hUrVnDLLbcA8Oijj/LjH/+40z5t21esWMGePXtSVldFRQUvvfRS3M9Kp6y/iCkiqRHk0hSrV69m7NixrF69mu985zudvt7U1MSppwYTN3Pnzu2yfcWKFVx00UUMGDAg4fMer+aKigo+/elPc9lllx23hqCpBy6SI4JamuLgwYNs2bKFZcuW8ZOf/KS1vaKigtLSUiZPnswFF1xARUUF48aNY+rUqZx33nnceOONuDsAzz//PMOHD2fw4MHMmjWLf/zjH50+Z/ny5QwaNIhRo0bx4osvtrYvXryY73//+532b2lfs2YNlZWV3HjjjQwbNoyGhgaqqqq4/PLLGTlyJBMmTGDv3r1AtMc/b948IpEIDz74ID//+c8ZPXo0w4cP58orr+SDDz6gtraWRx99lCVLljBs2DA2b97croZt27ZxySWXMGTIEL761a/y0UcftZ777rvvZtSoUQwaNIjNmzd36/cdFOAiOSOopSnWr1/P1VdfzaBBg+jXrx9VVUeXTdq6dSsPPvggO3fuBKC6upqlS5fy9ttv89577/Hiiy9y6NAhZs6cyVNPPcX27dtpamrikUceafcZe/fu5d577+XFF19ky5YtvP322wnXN3XqVCKRCKtWrWLbtm2ceuqp3HrrraxZs4aqqipmzZrFPffc07r/4cOHqays5K677mLs2LG88sorVFdXM336dB544AEGDhzI3LlzueOOO9i2bRulpaXtPu/rX/863/ve93jjjTcYPHhwu59ImpqaeO2111i6dGncn1SSpSEUkRwxoG8+dXHCurtLU6xevZrbb78dgOnTp7N69WpGjhwJwKhRo9rNaR41ahRFRUUADBs2jNraWnr37k1xcTGDBg0C4Oabb+bhhx9m3rx5rce9+uqrjBs3jpYF8W644YbW/xSSVVNTw5tvvslVV10FQHNzM/3792/9+g033ND6fvfu3dxwww3s3buXw4cPdzk/+8CBA+zfv5/LL7+89XuZNm1a69evu+46AEaOHEltbe0J1d+WAlwkR8yfUNJuDBy6vzTFX//6VzZu3Mj27dsxM5qbmzEzysvLATjttNPa7d+zZ8/W93l5eTQ1NZFu7s6FF17Iyy+/HPfrbWu+9dZbufPOO5k8eTIVFRUsXry4W5/d8v2n6nvXEIpIjpgyvJDvXjeYwr75GFDYN5/vXje4Wxcw16xZw4wZM3j//fepra3lT3/6E8XFxUmN75aUlFBbW8u7774LwBNPPNHag20xevRofve737Fv3z4aGxt5+umnk6qzd+/e/O1vf2v9vPr6+tYAb2xs5K233op73IEDBygsjP7+rFy5Mu752jr99NM544wzWr//eN9LKqkHLpJDUr00xerVq7n77rvbtV1//fWsXr263VDE8fTq1Yvly5czbdo0mpqauPjiizvN6ujfvz+LFy/m0ksvpW/fvgwbNiypOmfOnMncuXPJz8/n5ZdfZs2aNdx2220cOHCApqYm5s2bx4UXXtjpuMWLFzNt2jTOOOMMxo8fz65duwC49tprmTp1KuvXr+eHP/xhu2NWrlzJ3Llz+eSTTzj33HNZvnx5UrUmw1quAqdDJBLxeHM6ReTE7Nixg/PPPz/TZUgKxfszNbMqd4903FdDKCIiIaUAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkFKAi+SKLUth16b2bbs2RdtDrrKykttuuw3ovNRrd9XW1vLkk0/G/axMU4CL5IrCEfD0zKMhvmtTdLtwRCarOiEdb0OPRCI89NBDwIkF+PFua+8Y4G0/K9MU4CK5orgMpq2IhvbG+6Kv01ZE209QbW1t69Kw559/PlOnTuWTTz4B4i8R+/vf/751Qaf169eTn5/P4cOHOXToEOeeey4Af/jDH7j66qsZOXIkpaWlrQ+JaLmbcvTo0SxYsKBdHRUVFUyaNCnuUq/19fVcf/31XHzxxVx88cWtS9EuXryYGTNmMGbMGGbMmEFtbS2lpaWMGDGCESNGtP4nsHDhQjZv3sywYcNYsmRJ62dBdC2YKVOmMGTIEC655BLeeOON1nPPmjWLcePGce655wYX+O6etl8jR450EUmdt99+O/mDnv9f7vf2ib52065duxzwLVu2uLv7N77xDS8vL/eGhgYvKirympoad3efMWOGL1myxBsbG724uNjd3e+66y6PRCK+ZcsWr6io8OnTp7u7+/jx433nzp3u7v7KK6/4FVdc4e7uN998s0+cONGbmpo61fHCCy/4xIkT3d393nvv9fLy8tavfe1rX/PNmze7u/v777/v5513Xut+I0aM8E8++cTd3f/+9797Q0ODu7vv3LnTW/Kq7bk7bt9yyy2+ePFid3d//vnnfejQoa3nvvTSS/3QoUNeX1/vn/nMZ/zw4cMJ/Z7G+zMFKj1OpmotFJFcsmsTVC6DsgXR1+LSbvXAAc4++2zGjBkDwE033cRDDz3EVVdddcwlYj/3uc+xY8cOXnvtNe688042bdpEc3MzpaWlHDx4kJdeeqndEqxtH+4wbdo08vLykqrvt7/9bbv1wz/++GMOHjwIwOTJk8nPjy6n29jYyC233MK2bdvIy8tLaLnaLVu28MwzzwAwfvx49u3bx8cffwzAxIkT6dmzJz179uSzn/0sH3zwQetSuqnSZYCbWS9gE9Aztv8ad7/XzL4IlBMdhjkIzHT3d1Na3Zal0fG5tn/Bdm2Cuq0wdl5KP0rkpNcy5t0ybFJcmpJhlI5PUu/qyeplZWX88pe/pEePHlx55ZXMnDmT5uZmysvLOXLkCH379mXbtm1xj+24PG0ijhw5wiuvvEKvXr2Oe74lS5Zw1lln8frrr3PkyJG4+ycjHUvnJjIG/g9gvLsPBYYBV5vZJcAjwI3uPgx4Evj3lFd3El10Ecm4uq3tw7plTLxua7dO+8c//rF1adYnn3ySsWPHHneJ2NLSUpYuXcqll15KQUEB+/bto6amhosuuog+ffpQXFzculysu/P6668nVU/HpV6/9KUvtVsx8Fj/ORw4cID+/ftzyimn8MQTT9Dc3Bz3fG2VlpayatUqIDoOf+aZZ9KnT5+k6u2OLgM8NgRzMLbZI/bLY79aKj0dSN0jn1sEcNFFJGeNndf5305xWbd/mi0pKeHhhx/m/PPP56OPPuJb3/pWuyViBw8ezCmnnNK6ROzo0aP54IMPKCuL1jJkyBAGDx7c2nNftWoVy5YtY+jQoVx44YWsX78+qXquvfZann322daLmA899BCVlZUMGTKECy64gEcffTTucd/+9rdZuXIlQ4cO5Z133mntnQ8ZMoS8vDyGDh3KkiVL2h2zePFiqqqqGDJkCAsXLmy3Zng6JLScrJnlAVXA54GH3f1uMysF1gENwMfAJe7+cZxj5wBzAM4555yR77//fvJVbrwPNj0QHbcbf0/X+4vkiEwvJ1tbW8ukSZN48803M1bDySbly8m6e3NsqKQIGGVmFwF3ANe4exGwHPjBMY59zN0j7h5peZ5dUjpedOl4I4KISI5Kah64u+8HXgC+DAx191djX3oKuCy1pdH+osv4e44OpyjERbLCwIED1fvOoC4D3MwKzKxv7H0+cBWwAzjdzAbFdmtpS62ALrqInEwSGQaVcEj2zzKReeD9gZWxcfBTgJ+6+3Nm9i/AM2Z2BPgImJVssV2Kd3GluEwXMUVievXqxb59++jXr1+X0/cku7k7+/btS2r6YpcB7u5vAMPjtD8LPJtUhSKSUkVFRezevZv6+vpMlyIp0KtXr6Ru9tGdmCIh1qNHD4qLizNdhmSIFrMSEQkpBbiISEgpwEVEQiqhOzFT9mFm9cAJ3IoJwJnAX1JYTjqp9swIa+1hrRtUe1D+m7t3uhMyrQHeHWZWGe9W0jBQ7ZkR1trDWjeo9nTTEIqISEgpwEVEQipMAf5YpgvoBtWeGWGtPax1g2pPq9CMgYuISHth6oGLiEgbCnARkZAKRYCb2dVmVmNm75rZwkzXkwgzO9vMXjCzt83sLTO7PdM1JcvM8sys2syey3QtyTCzvma2xszeMbMdZnZppmtKlJndEfv78qaZrY49VDwrmdmPzOxDM3uzTdtnzOw3ZvZfsdczMlnjsRyj9vLY35k3zOzZlmW0s1nWB3hsGduHiT5E4gLga2Z2QWarSkgTcJe7XwBcAvxrSOpu63aCWOc9eA8Cv3L384ChhOR7MLNC4DYg4u4XAXnA9MxWdVwrgKs7tC0Ennf3LwDPx7az0Qo61/4b4CJ3HwLsBBalu6hkZX2AA6OAd939PXc/DPwE+EqGa+qSu+91962x938jGiKFma0qcWZWBEwEHs90Lckws9OBMmAZgLsfjj1JKixOBfLN7FTgUwTxsPAUcfdNwF87NH8FaHmy70pgSjprSlS82t391+7eFNt8hegjJLNaGAK8EPhTm+3dhCgIAcxsINE11V/tYtdsshRYABzJcB3JKgbqgeWx4Z/Hzey0TBeVCHevA74P/BHYCxxw919ntqqkneXue2Pv/wyclcliumEW8MtMF9GVMAR4qJnZp4FngHnu/nGm60mEmU0CPnT3qkzXcgJOBUYAj7j7cODvZO+P8e3Exou/QvQ/oQHAaWZ2U2arOnEenaMcunnKZnYP0SHQVZmupSthCPA64Ow220WxtqxnZj2Ihvcqd1+b6XqSMAaYbGa1RIesxpvZf2S2pITtBna3eeD2GqKBHgZXArvcvd7dG4G1BPGw8GB9YGb9AWKvH2a4nqSY2UxgEnCjh+AmmTAE+O+BL5hZsZn9E9GLOj/LcE1dsugDCpcBO9z9B5muJxnuvsjdi9x9INHf743uHoqeoLv/GfiTmZXEmr4IvJ3BkpLxR+ASM/tU7O/PFwnJBdg2fgbcHHt/M7A+g7UkxcyuJjpsONndP8l0PYnI+gCPXVS4BdhA9C/zT939rcxWlZAxwAyivddtsV/XZLqoHHErsMrM3gCGAf87s+UkJvZTwxpgK7Cd6L/PrL2928xWAy8DJWa228xmA/cDV5nZfxH9ieL+TNZ4LMeo/f8AvYHfxP69PprRIhOgW+lFREIq63vgIiISnwJcRCSkFOAiIiGlABcRCSkFuIhISCnARURCSgEuIhJS/x9hiOQe9But2wAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "#true largest eigenvalue is 50\n",
        "r = 0\n",
        "x = A @ x0.copy()\n",
        "x = x / la.norm(x)\n",
        "rs = []\n",
        "for i in range(k-1):\n",
        "    y = A @ x\n",
        "    r = x @ y\n",
        "    x = y / la.norm(y)\n",
        "    rs.append(r)\n",
        "print(r,max_ritz)\n",
        "pt.plot(ritz_max, \"o\", label=\"Arnoldi iteration\")\n",
        "pt.plot(rs, \"x\", label=\"power iteration\")\n",
        "pt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot convergence of Ritz values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Enable the Ritz value collection above to make this work."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzUElEQVR4nO3df1Bc15Un8O/dxLOMNGA1QpFxY6CFY0ZxG0tyF0g2RrOaiZXRqOKpSlBJrKacsvUjybBEJCpbU6mt7LqSXSkjG/0gY0tCWquiIEU4mSil8gRlpBHQFjTT1g94ioRN0w2iQRI/GoMls+PUnv2j+7V5dCNo+r1+79HnU0XBezaXkxR9fDn97jmCiMAYY8x8/pPeATDGGJsdTuCMMWZSnMAZY8ykOIEzxphJcQJnjDGT+mIif1hGRgbl5uYm8kcyxpjpffDBB4NEtGjy/YQm8NzcXLjd7kT+SMYYMz0hRHe0+1xCYYwxk+IEzhhjJsUJnDHGTIoTOGOMmRQncMYYM6kZJXAhhE8I0S6EuCqEcIfupQshfi+E+Cj02aJtqIwxI3q7wYNLnkHFvUueQbzd4FHnBzj3Ad5G5T1vY/B+nI5Jx9Da36q419rfimPSsbjXHqqpwb0Wl+LevRYXhmpq4l5bFssO/L8Q0TIicoSudwE4T0RfBnA+dM3YnOV0OuH1ehX3vF4vnE5n3Gt3dx/CcKBZcW840Izu7kNxrw0A1d134AyMKe45A2Oo7r4T99oFWQ+jvPZKOIlf8gyivPYKCrIejnttAIB1BVD3rc+TuLcxeG1dEffS9oV27GzYGU7irf2t2NmwE/aF9rjXTrE/BX9lZTiJ32txwV9ZiRT7U3GvHUZE034A8AHImHSvA0Bm6OtMAB3TrfPMM88QY1oavdhDn3YGFPc+7QzQ6MWeuNfu6uqiPXv2UFdXV9TreAwNX6KGRgcNDV+Keh2vpuFRWtrURk3Do1Gv4/V+5wAtf/0cvVF/k5a/fo7e7xxQZd2wrgaiPTai8z8Ofu5qUG1pV5+Lnj/5PB28fJCeP/k8ufpcqq39SXMLdaxcRXf376eOlavok+aWWa0DwE3RcnO0mxH/EuAFcBnABwC2he6NTPjnYuL1pO/dBsANwJ2dnT3r/yPY3OD6TR11t19T3Otuv0au39Spsv6nnQHyv34pnMQnX8dLTtrnz59XLXnL5KTd6XlT1eQtk5P2bk+fqslb9kb9Tcp57Sy9UX9T1XXDzv+Y6Edpwc8qO3j5INnfsdPBywdVX/vu/v30h/w/p7v79896jakS+ExLKMVEtALAXwP4eyFEyaRdPAGIOhmCiA4TkYOIHIsWRZwEZUnmkbwncHbfbvRIbQCAHqkNZ/ftxiN5T6iyfkreAqSXLcVw7Q18fM6H4dobSC9bipS8Baqsb7PZ4HA40NjYCIfDAZvNpsq6AJBuWQWrtQw+XzWs1jKkW1aptjYAFFtS8dKjGajqvoOXHs1AsSVVtbUveQZxwtWDijWP44SrJ6ImHjdvI+A+CpS8Gvw8uSYeh9b+VpzuOI3tBdtxuuN0RE08HvdaXAicPIWM734HgZOnImricYuW1R/0AeB/ANgJLqGwWepuv0Y/e2UTOX/5c/rZK5siduRqGKn30q3XGmmk3qvqurwDjySXT+SyyeTruMnlE7lsMvk6DnL5RC6bTL6Oh1w+kcsmk69jgdmWUADMB5A64etLAL4G4B8B7Ard3wXgp9OtxQncHD74nY9u3RxW3Lt1c5g++J1PtZ/h/OXPae+GvyHnL3+u2poyuWwyUu/VpHzCNXClty52RiTr9zsH6K2LnXGvTURETVWRybqrIXg/Tkfbj0Yka1efi462H4177cEjRyKS9SfNLTR45EjMa8WTwJcAuBb6uA7gh6H7CxF8+uQjAP8KIH26tTiBm8Otm8NU84PGcBKffB0vLXfgWtbAm5qaIpJ1V1cXNTU1xb22z/d2RLIeGr5EPt/bca9NRHTQdzsiWTcNj9JB321V1mfamiqBC0rgUGOHw0HcjdAcejsCqD8iwV5ihdTox9qtdmTlx/+ov1zzXr9jF7LtBRHX8RpruIWHslIVNe9xzwg+6x1D6urH4l6fzS3HpGOwL7SjMLMwfK+1vxXSkISX7S/HtfZQTQ1S7E9h/sqi8L17LS6MS+1YuGVLTGsJIT6gzx/hDuOTmCyqrHwL7CVWuN/zwV5iVSV5A8Btz4eKZJ1tL8D6Hbtw2/OhKuunrn4s4g3LlLwFnLxZVGZ/Dpx34CwqrXbgjMXMuS94aMc24eE3byPgvwwU74h7eTlpb8jfgNMdp7F39V7FjjwectK2bNqIwMlTsFZVKXbkM8U7cDZjcvJeu9WOoq8vwdqtdtQfkdDbEdA7NGZAmh+l1/AkJgAUZhZiQ/4GHGo7hA35G1RL3gAwf2URLJs2YvCf3oJl08ZZJe8H4QRuQlr3WLjrG1XsuLPyLVi71Y67vlFV1mdzi+ZH6W0lQOk7waR94SfBz6XvKHfkcUiq58Dj+eCnUNSh5vOlLDlo/RSK5kfpiTQ5iWn258B5B25C81cWwVpVBX9lJQYOHIC/snLWtTWWHJalzcO2675wQytnYAzbrvuwLG2eKus/m5eBzUXZOHChE5uLsvFsXoYq64ZpdBJTGpIUNe/CzELsXb0X0pAU99rjUrvidSm/bsel9rjXlvGbmCY2cOAABv/pLWR89ztYVFGhdzjM4OSk/dKjGTjeN4jDT+aqdpxeLptsLsrGCVcPqsuWq5fE5Zq3XDaZfJ0E+E3MOUbz2hqbc7TqhSIn7+qy5fj+C/moLluuqInHzX9Zmazlmrj/sjrrmxgncBOSH02yVlVhUUVFuJzCSdy8tO4HDgR34Mf7BlGZsxjH+wYj+oPPVlvvx4od97N5GaguW4623o9VWR/FOyJ32rYSVR4hNDtO4BrRctJHImprLLFS0wogSRXhJD4caIYkVSA1Lf7TqcDn5ZPDT+bitSWZOPxkrqImHo9vr86LKJc8m5eBb6/Oi3ttMzPaRB4WAy1PeC3csiXiDcv5K4tiPp7LjCPdsgp2+wFIUgU8XVWQpArY7QdUayl7dfS+ouZdbEnF4SdzcXX0virrs0iGmcij1keyPUao5aQPlnhaNrOSdXrepH89v4Q6PW+qtiabmpbdCIm0n8jDO3ANaXnCiyWe1WpFXV1deC6m1+tFXV0drFarKusPB5rh99ciN7ccfn9tRE2cqU/Lv5QBPolpalqe8GLRjTXcwrhnRHFv3DOCsYZbca9ts9lQWlqKuro6XLhwAXV1dSgtLVVlKo9c87bbDyBvSWW4nMJJHJpOpZef+97ZsBPVV6qxs2Gn6r1Q+CSmCWl5wotNTeuZmERE58+fpx/96Ed0/vx51dbUuh+4qWk4kUemxUxMPolpYlqe8GJT03omptfrhdvtRklJCdxud7icEq+cnO0Rb1imW1YhJ2e7KutrSfNmVibthZKQp8WiZXWtPpJpB86i03oqvUyLmZhajlQzM81nYspM1gtFTeAd+CQa1tXY1LSeSg8Ea973XP1IXfMY7rn6I2ris+X3+xU1b7km7vf7VVlfS9XddyKe+XYGxlDdfSfuteWDO+W1V/DmuY7wqUxV+6GYsBdKQkTL6lp9GGoHnoC6GovOrDMxzUzLocayN+pvUs5rZ+mN+puqrUlEpn2tJmKocfLuwDWuq7GpZdsL8PQL69Dyq1N4+oV1qszClH3WO6aoecs18c961Tk2blbywZ1t133Y09UfPpWpZj+UE64eVKx5HCdcPer1QQFM2wuFD/IkggZ1NfZgWu7A2YPt9vTR4gtXaLenT7U1E1YDNyE+yKMljepqZne5vjtifFpvRwCX67vjXnviFPrnNmzG+h27FDXxZMXNrPShZc8iQPuDPMm7AzdpXS0Rbt0cppofNNKtm8NRr+ORqKdQzGZo+BI1NDrCz4JPvo5XImrgZqT1Uyha78CTN4E3VUUm666G4H0WTtotZzyqJW/2YHLS7vS8qWryJtJ+pJqZadWzKBEHeXgiD5uS67ddcL/ng2NdLoq+vkTvcJKCp6sKPl81cnPLkbekUu9wkkb1lWocajuE7QXbUb68XJU1h2pqkGJ/SlE2udfiwrjUHnPnUJ7Iw2LS2xGA1OiHY10upEZ/RE2cqY+bWelDq5OYiWj7zAmcRejtCKD+iIS1W+0o+voSrN1qR/0RiZO4hriZlT7k7oN7V+9F+fLycGMrszSe4wTOItz1jWLtVjuy8i0AgKx8C9ZuteOub1TnyPTldDojep94vV44nc641x4bbVMMcJAHPIyNJvfTOQA0PTVt9pOYXANnbIbk/t/ycfrJ10wjPJV+yho4J3DGYiAnbYfDAbfbzck7UeSk7XgleGYjiZI3wG9iMqYKm80Gh8OBxsZGOBwOTt6JYisJJu/GnwY/J1HyfhBO4IzFQKt+4GameT9wgE9NT2HGCVwI8QUhxBUhxNnQtU0I4RJCdAohfimE+BO1g0vIL4YJDdXURIxmutfiwlBNjU4RJYeJNe81a9aEx6slexIvyHoY5bVXwq/VS55BlNdeQUHWw+r8gIk17zU//LwJncGTeCJep7HswL8H4MaE6z0AqojocQABAK+oFlWI5r8YJpWQLmcsgpn7gWtJ837g3I1watGOZ07+AJAF4DyANQDOAhAABgF8MfTPVwGon26d2RyllzubvVF/kzucTaBWj4W5ZvRiT0Tv7087AzR6sUefgJKIZv3ATcwo3Qj3AXgVwP8LXS8EMEJEfwxd9wKwRvtGIcQ2IYRbCOEeGBiI9b8veDYvA5uLsnHgQic2F2WrO+XDxDTvcmZSD2WlYrj2RngKz7hnBMO1N/BQljp9r81Ky4k8gMb9wE1M69fptAlcCLEewF0i+mA2P4CIDhORg4gcixYtivn7+RcjunstLgROnkLGd7+DwMlTEbW2ZKX1UGOzWpY2D9uu+8JJ3BkYw7brPixLmxf32nJps7psOb7/Qn64nGKG16rW7WQ1f51G25aTsnzyvxHcYfsA3AZwH8AvkIASCjeKj07NLmdzlRZDjc1ObiG729OnaivZty52Rrwm3+8coLcudqqyvpa0bCebiG6EMbWDBfAXAM6Gvq4DsDH09dsAvjvd98eawM38i6ElNWftzUXyHMyRei/Pw5xEi4k8ZqdVO9lEzMSM6SSmEOIvAOwkovVCiCUATgFIB3AFwGYi+r8P+n4+icm0Jte85bLJ5OtkJpdNXno0A8f7BlWdiWl2WrSTVZMqJzGJ6CIRrQ993UVEhUT0OBGVTpe8GUsEHmocnZy8Dz+Zi9eWZIYHHKs1Vs3MtGonmwhf1DsAxtSUuvqxiHspeQuSfvd9dfS+YsctT6m/Ono/qXfhE9vJFmYWovCRQsW10XEzK8ZY0jomHYN9oV2RrFv7WyENSXjZ/rKOkSlxN0LGGDMp7kbIGDMnDQc6mB0ncMaYsVlXKJtXyc2trCv0jMoQOIEzZgDd3Yci5l8OB5rR3X1Ip4gMRG5eVfct4MJPkm4az4NwAmfMAFLTChRDjOUhx6lpBTpHZhA80CEqTuCMGYA8xFiSKuDpqgpPqJeHHCc9HugQFSdwllCtZ95Fj6SctN4jtaH1zLs6RWQc6ZZVsFrL4PNVw2ot4+QtM+lAh0TgBM4S6pG8J3B23+5wEu+R2nB23248kveEzpHpbzjQDL+/Frm55fD7ayNq4knLpAMdEoGfA2cJJyftp19Yh2vn3sP6HbuQbU/uWq9c85bLJpOvWXLj58CZYWTbC/D0C+vQ8qtTePqFdUmfvAFgbLRNkazlmvjYaNs038mSGfdCYQnXI7Xh2rn3sPIbG3Ht3Ht47CsFSZ/Ec3K2R9xLt6zi3Td7IN6Bs4SSyyfrd+zCcxs2Y/2OXYqaOGNzhdGm0rMYaD2qyaxuez5U1Lyz7QVYv2MXbns+1DkyxtSViKn0nMA1Yl9ox86GneEkLrettC+06xyZvgpf/GZEuSTbXoDCF7+pU0SMaWP+yiJYq6rgr6zEwIED8FdWwlpVpepgY07gGinMLMTe1Xuxs2Enqq9Um6rHMGOxeLvBEzHA+JJnEG83eHSKyDh0n0rPZq8wsxAb8jfgUNshbMjfwMmbzUkFWQ8rptDLU+oLsh7WOTL9aT2VnhO4hsw8qomxmXo2LwPVZctRXnsFb57rQHntFVSXLcezeRl6h6YrueZtrarCooqKcDlFzSTOCVwjE0c1lS8vD5dTOImbl9PphNfrVdzzer1wOp06RWQcz+ZlYHNRNg5c6MTmouykT94AMC61K2reck18XGpX7WdwAteINCQpat5yTVwaknSOjM2W1WpFXV1dOIl7vV7U1dXBarXqHJn+LnkGccLVg4o1j+OEqyeiJp6MFm7ZElHznr+yCAu3bFHtZ/BResZiICdth8MBt9uN0tJS2Gw2vcPSlVzzlssmk6+NzOwzMXkHzlgMbDYbHA4HGhsb4XA4kj55A0Bb78eKZC3XxNt6P9Y5sumZ/XFf3oEzFgOz7sCru+9gWdo8FFtSw/ecgTFcHb2P8pzFOkamPzlpb8jfgNMdpw35uC/vwNmMXa7vRm9HQHGvtyOAy/XdOkVkDHLyLi0txZo1a1BaWqqoiRvZsrR52HbdB2dgDEAweW+77sOytHk6R6Y/Mz/uywmcRfhSbhrqj0jhJN7bEUD9EQlfyk3TOTJ9+f1+xY7bZrOhtLQUfr9f58imV2xJxeEnc7Htug97uvqx7boPh5/MVezIk5WZH/flEgqLSk7a9hIrpEY/1m61IyvfondYLE57uvpR1X0HlTmL8dqSTL3DmRnnvuAE+olzML2NwYEOxTviWnri476FmYUR10bBJRQWk6x8C+wlVrjf88FeYuXkPQc4A2M43jeIypzFON43GC6nGJ51hXKEmjxizboi7qXN/rgv78BZVLwDn1vkmrdcNpl8bXhy0na8EhxqPHHEmkEN1dQgxf6U4lnwey0ujEvtMT8LzjtwNmNy8l671Y6iry/B2q12RU3cyMYabmHcM6K4N+4ZwVjDLX0CMoiro/cVyVquiV8dva9zZDNkKwkm78afBj8bPHkDiWknyztwFuFyfTe+lJum2HH3dgRw1zeKFWtzdIxseuOeEQzX3kB62VKk5C2IuGYmZcIdOPB50rZs2ojAyVOzbic71Q6cR6qxCNGSdFa+xRQllJS8BUgvW4rh2huYX5SJe65+UyTv7u5DSE0rUIxQGw40Y2y0Leq4taQiJ285adueV14b2MR2shnf/U7i28kKIVKEEK1CiGtCiOtCiP8Zum8TQriEEJ1CiF8KIf5E1cgYm6WUvAWYX5SJsQu3ML8o0/DJGwBS0wogSRUYDjQD+HxKfWpacs8KBRB82mRisraVBK/9l/WMaka0bic7bQlFCCEAzCeiT4QQDwFwAvgegO8D+DURnRJCvA3gGhG99aC1uITCEkEum5hpBw58nrSt1jL4/bWKKfXMfCa2k52/sijiOhazfhOTgj4JXT4U+iAAawC8G7p/HMDfxhQRYxqYWPN++IXccDll8hubRpRuWQWrtQw+XzWs1jLTJG+eyBOdYdrJCiG+IIS4CuAugN8D8AAYIaI/hv6VXgBRe2oKIbYJIdxCCPfAwIAKITM2tc96xxQ7brkm/lmv8Z95Hg40w++vRW5uOfz+2nA5xeh4Ik90hmsnK4RYAOCfAfx3AO8Q0eOh+48B+BciemALLy6hMBadXD6RyyaTr41OTtqbi7JxwtVjilayZqLKc+BENALg3wCsArBACCE/xZIFwPgNIRgzqLHRNkWyTresgt1+AGOjbTpHNjNmnchzTDoW0fuktb8Vx6RjOkUUm5k8hbIotPOGEOJPAXwVwA0EE/k3Q//aSwDOaBQjY3NeTs72iJ12umWVaR4hNOtEHrP3A5/Jc+CZAI4LIb6AYMI/TURnhRB/AHBKCPFjAFcAHNUwTsaYQU2ewLMyb6FpJvLIvU+M3g98KtMmcCJqA7A8yv0uAOb4X8kY08yDJvIYPYEDyn7g2wu2myZ5A9wLhSVY65l30SMp67o9Uhtaz7w7xXcwo/v26ryIRP1sXga+vTpPp4hiY+Z+4JzATWiopibiRNe9FheGamp0imjmHsl7Amf37Q4n8R6pDWf37cYjeU/oHBlLRhP7f5cvLw+XU8ySxDmBm1AiupxpJdtegPU7duHsvt14//QJnN23G+t37EK2nY+Ms8TjfuAx4OfA1aNWlzO9vH/6BFp+dQorv7ERz23YrHc4M+J0OmG1WhVDjL1eL/x+P4qLi3WMjM113A98jpnY5cyyaaOpkneP1IZr597Dym9sxLVz70XUxI3KarUqhhjLQ46t1qiHkBnTHCdwk9K6y5lW5Jr3+h278NyGzeFyihmSuDzEuK6uDhcuXAhPqJ+4Izeq6u47ESPUnIExVHff0SmiGDj3fT5OTeZtDN5PcpzATWhiV7NFFRWwVlUpauJGdtvzoaLmLdfEb3s+1DmymbHZbHA4HGhsbITD4TBF8gaAZWnzsO26L5zE5ZFqy9Lm6RzZDGg4E1NLiXjYgBO4CSWiy5lWCl/8ZsQbltn2AhS++M0pvsNYvF4v3G43SkpK4Ha7w+UUo5NHqG277sOern5zzcOU+3/XfQu48BPTDHNIyMMGRJSwj2eeeYYMo6mKqKtBea+rIXifsSi6urpoz5491NXVFfXaDHZ7+mjxhSu029OndyixO/9joh+lBT+bxCfNLdSxchXd3b+fOlauok+aW2a1DgA3RcmpybsDN+mfZUw/fr9fUfOWa+J+vzn6uDkDYzjeN4jKnMU43jcYURM3NG9jcBZmyavBz5Nr4gal+cMG0bK6Vh+G2oETBXfce2zB/6LvsUXuyBmbI5qGR2lpUxs1DY9GvTY0+XUqvz4nXxsY78C1ZCsJTrlu/Gnws8FraozN1tXR+4qat1wTvzp6X+fIZsCkMzET8bBBch/kkcsmjleCf5aZ4I0Rxpg5DNXUIMX+lKJscq/FhXGpPeapPFMd5EneBC4nbzlpT75mjDGD4JOYk5n0zzLGGJPNZKDD3FS8I/KerYR334wx00jeHThjjJkcJ3DGGDMpTuCMMWZSnMAZY8ykOIEzxpLWMelYxPi01v5WHJOO6RRRbDiBM8aSln2hXTEDU56RaV9o1zmymeEEzhiLy9sNHlzyDCruXfIM4u0Gj04RzZw8A3Nnw05UX6kODziWZ2QaHSdwxlhcCrIeRnntlXASv+QZRHntFRRkPaxzZDNTmFmIDfkbcKjtEDbkbzBN8gY4gTNmCN3dhzAcaFbcGw40o7v7kE4RzdyzeRmoLluO8torePNcB8prr6C6bDmezcvQO7QZae1vxemO09hesB2nO05H1MSNjBM4YwaQmlYASaoIJ/HhQDMkqQKpaQXTfKcxPJuXgc1F2ThwoRObi7JNlbzlskn58vJwOcUsSZwTOGMGkG5ZBbv9ACSpAp6uKkhSBez2A0i3rNI7tBm55BnECVcPKtY8jhOunoiauFFJQ5Ki5i3XxKUhSefIZiZ5uxEyZkCerir4fNXIzS1H3pJKvcOZEbnmLZdNJl+z+HE3QsYMbjjQDL+/Frm55fD7ayNq4kbV1vuxIlnLNfG23o91jmzu4x04YwYg17zlssnka5bceAfOksJYwy2Me0YU98Y9IxhruKVPQDM0NtqmSNZyTXxstE3nyJiRTZvAhRCPCSH+TQjxByHEdSHE90L304UQvxdCfBT6bNE+XMYe7KGsVAzX3ggn8XHPCIZrb+ChrFR9A5tGTs72iJ12umUVcnK26xQRM4OZ7MD/COAHRPQVACsB/L0Q4isAdgE4T0RfBnA+dM2YrlLyFiC9bCmGa2/g43M+DNfeQHrZUqTkLdA7NJZkhmpqIgYY32txYaimRrWfMW0CJ6J+Iroc+noMwA0AVgAvAjge+teOA/hb1aJiLA4peQswvygTYxduYX5RJidvposU+1OKKfTylPoU+1Oq/YyYRqoJIXIBLAfgArCYiPpD/+g2gMVTfM82ANsAIDs7e9aBMjZT454R3HP1I3XNY7jn6sd/zlvASZwl3PyVRbBWVcFfWQnLpo0InDwFa1WVYkp9vGb8JqYQ4s8A/ArADiIanfjPKPgoS9THWYjoMBE5iMixaNGiuIJliXG5vhu9HQHFvd6OAC7Xd+sU0czJNe/0sqV4+IXccDll8hubzESc+wBvo/KetzF43+DmryyCZdNGDP7TW7Bs2qhq8gZmmMCFEA8hmLx/QUS/Dt2+I4TIDP3zTAB3VY2M6eZLuWmoPyKFk3hvRwD1RyR8KTdN58im91nvmKLmLdfEP+sd0zcwNnvWFUDdtz5P4t7G4LV1hZ5Rzci9FhcCJ08h47vfQeDkqYiaeLymfQ5cCCEQrHEPE9GOCff/EcAQEe0WQuwCkE5Erz5oLX4O3DzkpG0vsUJq9GPtVjuy8vlBI6YTOWk7XgHcR4HSdwBbid5RPZBc85bLJpOvYxHPc+DPAfg7AGuEEFdDH+sA7AbwVSHERwD+KnTN5oisfAvsJVa43/PBXmLl5G1y1d134Awo/wpxBsZQ3X1Hp4hiZCsJJu/GnwY/Gzx5A8C41K5I1nJNfFxqV+1nTPsmJhE5AYgp/vFfqhYJM5TejgCkRj8c63IhNfphzbdwEjexZWnzsO26D4efzEWxJRXOwFj42hS8jcGdd8mrwc+25w2fxBdu2RJxb/7KIlXr4DE9hcKSg1w+kcsm1nyL4pqZT7ElFYefzMW26z689GgGjvcNhpO54cnlE7lsYnteeZ3E+Ci9Rsw8LPWub1SRrLPyLVi71Y67vtFpvnNuczqd8Hq9interxdOp1OniGJTbEnFS49moKr7Dl56NMMcyRsA/JeVydpWErz2X457aTO/TgFO4Jox87DUFWtzInbaWfkWrFibo1NExmC1WlFXVxdO4l6vF3V1dbBarTpHNjPOwBiO9w2iMmcxjvcNRtTEDat4R+RO21YSvB8nM79OAe5GqCn5l2FD/gac7jhtqmGpLDo5aTscDrjdbpSWlsJms+kd1rQm1rwn18BNsxPXiBlep9yNUAdmHpbKorPZbHA4HGhsbITD4TBF8gaAq6P3FclarolfHb0f99pmnkoPmPt1yglcQ2YelqqV1jPvokdStkjtkdrQeuZdnSKKjdfrhdvtRklJCdxud0RN3KjKcxZH7LSLLakoz4naASMmZp9Kb+rXKREl7OOZZ56hZOHqc9HzJ58nV58r6nWy6m6/Rj97ZRN1t1+Lem1kXV1dtGfPHurq6op6ncze7xyg5a+fozfqb9Ly18/R+50Deoc0I1q+TgePHKFPmlsU9z5pbqHBI0diXguAm6LkVE7gGjnafjTil8DV56Kj7Ud1isg45KTt/OXPTZO8iYiampoiknVXVxc1NTXpFJGxvFF/k3JeO0tv1N/UO5QZ0/J1+klzC3WsXBVO4pOvYzFVAuc3MZku3j99Ai2/OoWV39iI5zZs1jsc3XV3H0JqWoFiqMNwoBljo22mGOogl002F2XjhKuHBxqHyMfn4+1GyG9iMsPokdpw7dx7WPmNjbh27r2ImngySk0rgCRVhAcZyzMxU9MKdI5sehOn0H//hXxUly1X1MSTmSG6ETKmlh6pDWf37cb6Hbvw3IbNWL9jF87u2530SVyegSlJFfB0VZlqoDFPpZ+a1t0IuQbOEsr1m7qImnd3+zVy/aZOp4iMpdPzJv3r+SXU6XlT71CMo6mKqKtBea+rIXjfwBJRA+cdOEuowhe/iWy7siyQbS9A4Yvf1Cki4xgONMPvr0Vubjn8/tpwOSXpmbQfeCK6EfKbmIwZgFzzlssmk6+Tngn7gauJ38RkzMDGRtsUyVquiY+NJvd7A2Em7AeeCNxOljEDiPaoYLplFe++ZSbsB54IvAM3oaGamoh3s++1uDBUU6NTRMYx1nArYoDxuGcEYw239AmIxW9iP/A1Pwx+nlgTT2KcwE0oxf4U/JWV4SQuHxZIsT+lc2T6eygrVTGFXp5S/1BWcnfcMzUN+4GbHb+JaVJqnfCai+SkPb8oE/dc/Yop9cmquvsOlqXNUzS0cgbGcHX0vioNrZi2+E3MOUbrE15mlpK3APOLMjF24RbmF2UmffIGPp+JKQ9xkPuBL0ubp3Nk+uKJPEwXmp/wMrFxzwjuufqRuuYx3HP1R9TEk9HEmZh7uvp5mEOI2Sfy8ElME1LzhNdc82lngPyvX6JPOwNRr5Pdbk8fLb5whXZ7+vQOxTDkFrIHLx80bMtn8EnMuSMRJ7zM6rPeMUXNOyVvAdLLluKzXpPMf9SQaWdiaszME3n4TUzGkgDPxJwaz8RkjBmaljMxzUxO3ntX70X58nLsXb1XURM3Ot6BM8aS1jHpGOwL7Yodd2t/K6QhCS/bX9YxMiXegTPGNGHmqfQv21+OKJcUZhaqkrwTcWLa0AnczL8YjCULs0+l10oiTkwbOoHzLwYzEqfTCa/Xq7jn9XrhdDp1isgY5Ak85bVX8Oa5jvB4tWSfiSk/HeavrMTAgQPwV1aqfmLa0AmcfzGYkVitVtTV1YWTuNfrRV1dHaxWq86R6e/ZvAxsLsrGgQud2FyUza/RkKSficm/GMwobDYbSktLUVdXhwsXLqCurg6lpaWw2Wx6h6a7S55BnHD1oGLN4zjh6uGBxiFan5g2fALnXwxmJDabDQ6HA42NjXA4HJy8wVPppyLXvK1VVVhUUREup6iZxKdN4EKIY0KIu0IIacK9dCHE74UQH4U+W1SLaAL+xdDH5fpu9HYEFPd6OwK4XN+tU0TG4fV64Xa7UVJSArfbHVETT0aaT6V37ovs/e1tDN43sIScmI52vn7iB4ASACsASBPu/RTArtDXuwDsmW4dmkUvlLcudtL7nQOKe+93DtBbFztj7CTAYnHr5jDV/KCRbt0cjnqdrLq6umjPnj3U1dUV9ZpppKuBaI/t88n0k6+TAKbohTKjgzxCiFwAZ4nIHrruAPAXRNQvhMgEcJGI8qdbhw/ymEdvRwD1RyTYS6yQGv1Yu9WOrHxN/tAyDafTCavVqiibeL1e+P1+FBcX6xhZEuChxqoe5FlMRP2hr28DmLIjvBBimxDCLYRwDwwMzPLHsUTLyrfAXmKF+z0f7CXWpE/eAFBcXBxR87bZbKok7+7uQxgONCvuDQea0d19KO615wQeahxV3G9ihrb3U27jiegwETmIyLFo0aJ4fxxLkN6OAKRGPxzrciE1+iNq4kxdqWkFkKSKcBIfDjRDkiqQmlagc2QGMXmoMc/DBDD7qfR3hBCZE0ood9UMiulLLp/IZRNrvkVxzdSXblkFu/0AJKkCVmsZ/P5a2O0HeCo9oBxqbCsJTqSfeJ3EZrsD/y2Al0JfvwTgjDrhMCO46xtVJOusfAvWbrXjrm9U58jmtnTLKlitZfD5qmG1lnHylmk41HjOj1QTQpwE0AwgXwjRK4R4BcBuAF8VQnwE4K9C12yOWLE2J2KnnZVvwYq1OTpFlByGA83w+2uRm1sOv782oiaetIp3RO60bSXB+3Ey+0g1bifLmAHINW+5bDL5mmmHBzowNkOtZ95Fj9SmuNcjtaH1zLs6RWQMY6NtimQt18THRtum+c6Zqe6+EzFCzRkYQ3X3HVXWNzMzj1TjBM4S6pG8J3B23+5wEu+R2nB23248kveEzpHpKydne8ROO92yCjk521VZf1naPGy77gsncXmk2rK0eaqsb2at/a043XEa2wu243THadWm8SSiHzhPpWcJ191+jX72yiZy/vLn9LNXNlF3+zW9Q0oKTcOjtLSpjXZ7+mhpUxs1DY/qHZLu5In08iT6ydfx+KS5hTpWrqJPmluiXscCPJWeGUW2vQBPv7AOLb86hadfWIdsOz/rnAjFllS89GgGqrrv4KVHM5J+mDEASEOSouZdmFmIvav3QhqSpvnO6SV9P3BNmbRBzlzQI7Xh2rn3sPIbG3Ht3HsRNXGmDWdgDMf7BlGZsxjH+wYjauKzZebJWVqOVAO4H7h2rCuChwHkJC4fFrCu0DOqOU+uea/fsQvPbdiM9Tt2KWri8RpruIVxz4ji3rhnBGMNt1RZ36zkmvfhJ3Px2pJMHH4yV1ETjwdPzpqa1v3Ak7sGLnc1O//jpOtuphfXb+oiat7d7dfI9Zs6Vdb/tDNA/tcv0aedgajXyeqg73ZEzbtpeJQO+m6rsv77nQO0/PVz9Eb9TVr++rmILqLJKBE1cH4O/MJPgg1ySl4F1vxQ72iYCsY9IxiuvYH5RZm45+pHetlSpOQt0DusOe/Ncx04cKETFWsex/dfmLY56cw59wX/Mp54mMfbGDyJqcJhHq0M1dQgxf6Uomxyr8WFcakdC7dsiWktfg48Gm6QMyel5C3A/KJMjF24hflFmZy8E0DTyVkmLXcu3LIlouY9f2VRzMn7QZI3gU9skLPmh8HPE39JmGmNe0Zwz9WP1DWP4Z6rP6ImPls8lT46zSdnyb1P6r4V/IuZG1mFJW8C17BBDtOPXD5JL1uKh1/IRXrZUgzX3lAlifNU+ug0H6kGaNYP3OzNrLgGrpFj0jHYF9oVjyi19rdCGpJUe0SJRRpruIWHslIVZZNxzwg+6x1D6urH4l5fTtoOhwNut5un0ieKRhN55D4o8rPgk6+NgmvgCWb2Lmdmlbr6sYiad0reAlWSN6DdVHqeyPMAGpY75YM7Oxt2ovpKtarJOxFH6TmBa0TLXwymH62m0vNEngfQuNypVTOrFPtT8FdWhpP4vRYX/JWVSLE/pcr6AJL8OfAEOHj5INnfsdPBywdVW3PwyJGIZ0k/aW6hwSNHVPsZLJLWU+mHhi9RQ6ODOj1vUkOjg4aGL6myLnswuf/JwcsHVeuDIpOf/b67f/+snwEn4l4outCqy1lC/svOIvj9fkXN22azobS0FH6/X5X1tZzIw+1ko5tY8y5fXh7+q1mt16rWR+l5B64RLbucEan3X/ZoPvidj27dHFbcu3VzmD74nU+1n8EiabkDlzsRyqcxJ1/H462LnREnL9/vHKC3LnbGvbbWjrYfjXhNuvpcdLT9qCrra70D5wSuEa1/MYiI7u7fT3/I/3O6u3+/amsSBZN1zQ8aw0l88nU8tD5Kb1Zy8paT9uRrNWjVTlY+Ri8n8cnXcWuqimxz0dUQvG9giThKzwncpLTcgRN9nrRbznhUS95En/cCl5P45Ot4jV7sieh78mlngEYv9qiyvlZ8vrcjkvXQ8CXy+d5W9efs9vTR4gtXaLenT9V1Ne2FIvcskpP45Os4aLnRUvO9Kk7gc4ia/2V/kJYzHqrefp5aznhUXVfLgQ5aNrNqamqKeMOyq6uLmpqa4l47EbQe6PBG/U3Kee0svVF/U9V1iUizxnNaljoTkcD5TUwTGpfaFY3h5cbx41K7aj+jtyMAqdEPx7pcSI1+9HYEVFtby4EOKXkLwqcvPz7nC5/KVKMfipYnMbV+DlzLdrKAxr1QAM1OYmr5uC8/Rsh0oWUNnCgxI9VG6r1067VGGqn3qrqu/Ojg+fPnNXmEUKsauJbtZDWvgRNp3vpZi8d9ifgxQqaDu75RrN1qR1a+BQCQlW/B2q123PWNxr221gMdAO2aWQHancSUp9BLUgU8XVWQpArFlHoj07wXisaN57R63BfgiTxMByvW5oSTtywr34IVa3PiXvu250Os37ErXDbJthdg/Y5duO35MO61AW2bWQHancQEtH0O3NRT6TU8ian1c+A8kYexGGj5FIrZT2LyY4SRtHwKhafSMxYjLZtZaXkSU+59YrcfQN6SynA5ZfIbm/HQaiq9XDIpr72CN891hHuDyyWVuJl0oEMiHjbgBM7YDBUXF0fUvG02G4qLi+Nee2y0TVHzlmviY6PqvTeg1VR6IJjENxdl48CFTmwuylYveQOaDnQwfdfQaNtyrT64hMKYPrQ8Sk+UoKHG539M9KO04GcVadXMiksojDFVXB29j8NP5obLJsWWVBx+MhdXR+/HvbbmI9UATefXatVOVi6Z+CsrMXDgAPyVlYqSiho4gTOWBMpzFkfUvIstqSjPWRz32vwY4dS4GyFjLLlp+BSKWbqGYooSCs/EZIwlLS1n18pH5+WyyeTrWEw1EzOuBC6E+BqA/QC+AKCGiHY/6N/nBM4YSxZDNTVIsT+lSNb3WlwYl9qxcMuWmNZSPYELIb4A4EMAXwXQC+DfAWwioj9M9T2cwBljLHZaTKUvBNBJRF1E9B8ATgF4MY71GGOMxSCeBG4FcGvCdW/onoIQYpsQwi2EcA8MDMTx4xhjjE2k+WOERHSYiBxE5Fi0aJHWP44xxpJGPAncD2Big4ms0D3GGGMJEE8C/3cAXxZC2IQQfwJgI4DfqhMWY4yx6cT7GOE6APsQfIzwGBH9ZJp/fwBA9yx/XAYAlec0JQzHrg+zxm7WuAGOXSs5RBRRg07oQZ54CCHc0R6jMQOOXR9mjd2scQMce6JxLxTGGDMpTuCMMWZSZkrgh/UOIA4cuz7MGrtZ4wY49oQyTQ2cMcaYkpl24IwxxibgBM4YYyZligQuhPiaEKJDCNEphNildzwzIYR4TAjxb0KIPwghrgshvqd3TLESQnxBCHFFCHFW71hiIYRYIIR4VwhxUwhxQwixSu+YZkoIURn6fZGEECeFECl6xzQVIcQxIcRdIYQ04V66EOL3QoiPQp8tesY4lSli/8fQ70ybEOKfhRALdAxxRgyfwENta38G4K8BfAXAJiHEV/SNakb+COAHRPQVACsB/L1J4p7oewBu6B3ELOwH8Dsi+nMAT8Mk/xuEEFYAFQAcRGRH8IDcRn2jeqB3AHxt0r1dAM4T0ZcBnA9dG9E7iIz99wDsRFSAYKvsf0h0ULEyfAKHSdvWElE/EV0OfT2GYBKJ6NZoVEKILAB/A6BG71hiIYR4GEAJgKMAQET/QUQjugYVmy8C+FMhxBcBzAPQp3M8UyKiRgDDk26/COB46OvjAP42kTHNVLTYiegcEf0xdNmCYH8nQzNDAp9R21ojE0LkAlgOwKVzKLHYB+BVAP9P5zhiZQMwAOD/hMo/NUKI+XoHNRNE5AewF0APgH4AHxPROX2jitliIuoPfX0bQPxTk/XxMoB/0TuI6ZghgZuaEOLPAPwKwA4iGtU7npkQQqwHcJeIPtA7lln4IoAVAN4iouUA7sG4f8YrhOrFLyL4H6FHAcwXQmzWN6rZCw3jNd1zykKIHyJYAv2F3rFMxwwJ3LRta4UQDyGYvH9BRL/WO54YPAfg60IIH4IlqzVCiBP6hjRjvQB6iUj+a+ddBBO6GfwVAC8RDRDRZwB+DeBZnWOK1R0hRCYAhD7f1TmemAghvgVgPYD/SiY4JGOGBG7KtrVCCIFgHfYGEb2pdzyxIKJ/IKIsIspF8P/vC0Rkip0gEd0GcEsIkR+69ZcAppzTajA9AFYKIeaFfn/+EiZ5A3aC3wJ4KfT1SwDO6BhLTEJD2l8F8HUiuq93PDNh+AQeelOhHEA9gr/Mp4nour5RzchzAP4Owd3r1dDHOr2DShL/DcAvhBBtAJYB+F/6hjMzob8a3gVwGUA7gq9Pwx7vFkKcBNAMIF8I0SuEeAXAbgBfFUJ8hOBfFLv1jHEqU8ReDSAVwO9Dr9e3dQ1yBvgoPWOMmZThd+CMMcai4wTOGGMmxQmcMcZMihM4Y4yZFCdwxhgzKU7gjDFmUpzAGWPMpP4/KHnDFwC1pVUAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "for i, rv in enumerate(ritz_values):\n",
        "    pt.plot([i] * len(rv), rv, \"x\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}