ACVI_Higher_dimension_bg.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# High dimensional bilinear game on the simplex\n",
        "\n",
        "Results on the following constrained bilinear game, where $\\mathbf{x}_1, \\mathbf{x}_2 \\in \\mathbb{R}^{500}$, and $\\eta \\in (0,1)$:\n",
        "\n",
        "$ \n",
        "min_{\\mathbf{x}_1} max_{\\mathbf{x}_2} \\quad \\eta \\cdot \\mathbf{x}_1^\\intercal + (1-\\eta) \\cdot \\mathbf{x}_1^\\intercal \\mathbf{x}_2 - \\eta \\cdot \\mathbf{x}_2^\\intercal \\mathbf{x}_2\n",
        "$\n",
        "\n",
        "$\n",
        "\\hspace{2em} \\text{s.t.} \\qquad \\mathbf{x} \\geq \\mathbf{0}\n",
        "$\n",
        "\n",
        "$\n",
        "\\hspace{5em} \\mathbf{e}^\\intercal \\mathbf{x}_1 = 1\n",
        "$\n",
        "\n",
        "$\n",
        "\\hspace{5em} \\mathbf{e}^\\intercal \\mathbf{x}_2 = 1\n",
        "$"
      ],
      "metadata": {
        "id": "_GkYPL_Y3C1l"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title install latex for plots description\n",
        "! sudo apt install texlive-latex-extra\n",
        "! sudo apt install dvipng\n",
        "! sudo apt-get install texlive-full"
      ],
      "metadata": {
        "cellView": "form",
        "id": "bjpSpMxL3CLo"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "M3dAiwtRbbO0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "cellView": "form",
        "outputId": "6210775a-4cc4-4287-847c-e594cd12bc54"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Init [[0.00221029]\n",
            " [0.00288036]]\n"
          ]
        }
      ],
      "source": [
        "#@title problem setup\n",
        "import numpy as np\n",
        "import scipy.optimize as opt\n",
        "from scipy import linalg\n",
        "import matplotlib.pyplot as plt\n",
        "from cvxopt import matrix,solvers\n",
        "solvers.options['show_progress'] = False\n",
        "\n",
        "eps = .02\n",
        "dim = 500  # dim(x1) == dim(x2) == dim\n",
        "x_opt = np.ones((2*dim,1))/dim\n",
        "_seed = 0\n",
        "\n",
        "\n",
        "temp1 = np.concatenate((np.ones((1,dim)), np.zeros((1,dim))),axis=1)\n",
        "temp2 = np.concatenate((np.zeros((1,dim)), np.ones((1,dim))),axis=1)\n",
        "C = np.concatenate((temp1,temp2), axis=0)\n",
        "d = np.ones((2,1))\n",
        "temp = np.linalg.inv(np.dot(C, C.T))\n",
        "temp = np.dot(C.T, temp)\n",
        "dc = np.dot(temp, d)\n",
        "Pc = np.identity(2*dim) - np.dot(temp,C)\n",
        "\n",
        "G = matrix(-np.identity(2*dim), tc='d')\n",
        "h = matrix(np.zeros(((2*dim,1))), tc='d')\n",
        "AA = matrix(C, tc='d')\n",
        "b = matrix(d, tc='d')\n",
        "\n",
        "def f_y(y,c,mu): \n",
        "  \"\"\" Update y in Alg. ACVI. \"\"\"\n",
        "  return -mu * np.sum(np.log(y)) + beta/2 * np.sum((y-c)**2)\n",
        "\n",
        "def proj(z):\n",
        "  P = matrix(np.identity(2*dim), tc='d')\n",
        "  q = matrix(-z, tc='d')\n",
        "  res = solvers.qp(P, q, G, h, AA, b)\n",
        "  return res['x']\n",
        "\n",
        "np.random.seed(_seed)\n",
        "init = np.random.rand(2*dim, 1)\n",
        "init[:dim] = init[:dim] / np.sum(init[:dim])\n",
        "init[dim:] = init[dim:] / np.sum(init[dim:])\n",
        "\n",
        "print('Init {}'.format(init[:2]))"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## For fixed $\\epsilon$ and varying rotation component $\\eta$ ($x$-axis), compute number of iterations ($y$-axis) needed to be at most $\\epsilon$ far away from the solution"
      ],
      "metadata": {
        "id": "BBhOhgFJ16Oy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from numpy.core.arrayprint import format_float_scientific\n",
        "\n",
        "eps = .02\n",
        "gamma = .3  # step-size for EG, OGDA, LA\n",
        "max_iters = 500\n",
        "T, K = 10, 1\n",
        "last_K = max_iters - (T-1) * K\n",
        "n = (T-1) * K + last_K\n",
        "num_acvi, num_gda, num_eg, num_la, num_ogda = [], [], [], [], []  \n",
        "rot = np.linspace(.01, .99, 5)  #  rot = np.linspace(.01, .99, 30) \n",
        "\n",
        "for a in rot:\n",
        "  A1 = np.concatenate((a*np.identity(dim), (1-a)*np.identity(dim)), axis=1)\n",
        "  A2 = np.concatenate((-(1-a)*np.identity(dim), a*np.identity(dim)), axis=1)\n",
        "  A = np.concatenate((A1, A2), axis=0)\n",
        "\n",
        "  # ACVI ------------\n",
        "  beta, mu, delta = .5, 1e-6, .5\n",
        "  y, lmd = init, np.zeros(((2*dim,1)))\n",
        "  finished, cnt = False, 0\n",
        "  for t in range(0, T):\n",
        "    mu *= delta\n",
        "    flag = 1\n",
        "    for k in range(0, K if t+1 < T else last_K):\n",
        "      cnt += 1\n",
        "      x = linalg.solve(np.identity(2*dim) + np.dot(Pc,A)/beta,\n",
        "                       np.dot(Pc, y-lmd/beta) + dc)\n",
        "      if np.linalg.norm(x-x_opt)/np.linalg.norm(x_opt) <= eps or \\\n",
        "         cnt == max_iters:\n",
        "        finished = True\n",
        "        num_acvi.append(cnt)\n",
        "        break\n",
        "      res = opt.minimize(lambda z: f_y(z,x+lmd/beta,mu), y, method='SLSQP')\n",
        "      if not res.success:\n",
        "        print('Failed to update y. t={}. k={}'.format(t, k))\n",
        "        flag = 0\n",
        "        break\n",
        "      y = res.x.reshape(-1,1)\n",
        "      lmd += beta*(x-y)\n",
        "    if flag == 0 or finished:\n",
        "      break\n",
        "\n",
        "  # GDA ------------\n",
        "  z = np.copy(init)\n",
        "  for t in range(n):\n",
        "    z = proj(z - gamma * np.dot(A,z))\n",
        "    if np.linalg.norm(z-x_opt)/np.linalg.norm(x_opt) <= eps or t+1 >= max_iters:\n",
        "      num_gda.append(t+1)\n",
        "      break\n",
        "\n",
        "  # EG ------------\n",
        "  z = np.copy(init)\n",
        "  for t in range(n):\n",
        "    zz = proj(z-gamma*np.dot(A,z))\n",
        "    z = proj(z-gamma*np.dot(A,zz))\n",
        "    if np.linalg.norm(z-x_opt)/np.linalg.norm(x_opt) <= eps or \\\n",
        "       (t+1)*2 >= max_iters:\n",
        "      num_eg.append(2*(t+1))\n",
        "      break\n",
        "\n",
        "  # OGDA ------------\n",
        "  z0 = np.copy(init)\n",
        "  z1 = np.copy(init)\n",
        "  for t in range(n):\n",
        "    z = proj(z1 - 2*gamma*np.dot(A,z1) + gamma*np.dot(A,z0))\n",
        "    z0 = np.copy(z1)\n",
        "    z1 = np.copy(z)\n",
        "    if np.linalg.norm(z-x_opt) / np.linalg.norm(x_opt) <= eps or \\\n",
        "       t+1 >= max_iters:\n",
        "      num_ogda.append(t+1)\n",
        "      break\n",
        "\n",
        "  # LA ------------\n",
        "  z = np.copy(init)\n",
        "  la_k, la_alpha = 5, 0.5\n",
        "  for t in range(n):\n",
        "    z_old = np.copy(z)\n",
        "    for _ in range(la_k):\n",
        "      #z = z - gamma * np.dot(A,z)\n",
        "      z = proj( z - gamma * np.dot(A, z))\n",
        "    z = proj(la_alpha * z_old + (1-la_alpha) * z)\n",
        "    if np.linalg.norm(z - x_opt) / np.linalg.norm(x_opt) <= eps or \\\n",
        "       (t+1) * la_k >= max_iters:\n",
        "      num_la.append(min(la_k*(t+1), max_iters))\n",
        "      break"
      ],
      "metadata": {
        "id": "6OZBWWsDEU9n"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ms = 3\n",
        "print('GDA:    {}'.format(num_gda))\n",
        "print('EG:     {}'.format(num_eg))\n",
        "print('OGDA:   {}'.format(num_ogda))\n",
        "print('LA-GDA: {}'.format(num_la))\n",
        "print('ACVI:   {}'.format(num_acvi))\n",
        "\n",
        "plt.plot(rot, num_gda, 'ys-', label='GDA', markersize=ms)\n",
        "plt.plot(rot, num_eg, 'ch-', label='EG', markersize=ms)\n",
        "plt.plot(rot, num_ogda, 'mD-',label='OGDA', markersize=ms)\n",
        "plt.plot(rot, num_la, 'bo-', label='LA{}-GDA'.format(la_k), markersize=ms)\n",
        "plt.plot(rot, num_acvi, 'gv-', label='ACVI', markersize=ms)\n",
        "plt.xlabel(r'$\\eta$')\n",
        "plt.ylabel('number of iterations')\n",
        "plt.legend(loc=\"upper right\")\n",
        "plt.savefig('acvi_hd_bilinear_eta.pdf', bbox_inches='tight')\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 370
        },
        "id": "_amS4pZuEne9",
        "outputId": "301a5c22-16e3-4b5a-a4f1-ebee0c6e2355"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "GDA:    [500, 66, 23, 14, 10]\n",
            "EG:     [150, 64, 44, 36, 30]\n",
            "OGDA:   [63, 32, 23, 17, 14]\n",
            "LA-GDA: [95, 60, 45, 40, 35]\n",
            "ACVI:   [18, 9, 6, 6, 5]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU5b348c93lswkZGFJCCQTBGSRRQ2IAlapa69aK1RaRdtKra3XCgrU3tb23p/X1tpqq1Wpu7VKN6BqXapWa6FYvQIKssmOCGQSEkJYQjJJZnt+f8zJkECWSTKTyfJ9v17nNWc/35PAfPM85znPI8YYlFJKKQBbsgNQSinVdWhSUEopFaVJQSmlVJQmBaWUUlGaFJRSSkU5kh1AR2RnZ5uhQ4cmOwyllOpW1q5de9AYk9PUtm6dFIYOHcqaNWuSHYZSSnUrIrK3uW1afaSUUipKk4JSSqkoTQpKKaWiuvUzBaVUzxYIBPB6vdTW1iY7lG7J7Xbj8XhwOp0xH6NJQSnVZXm9XjIyMhg6dCgikuxwuhVjDBUVFXi9XoYNGxbzcQmtPhKRPSKySUTWi8gaa11/EXlHRHZan/2s9SIiC0Vkl4hsFJGJiYxNKdX11dbWMmDAAE0I7SAiDBgwoM2lrM54pnChMabQGDPJWr4TWGaMGQkss5YBLgdGWtPNwBOdEJtSqovThNB+7fnZJaP6aDpwgTW/CFgB/NBa/3sT6ct7lYj0FZHBxpj98Q7g/fcHEAweAsDpzOVznyuN9yWUUqpbSnRJwQD/EJG1InKztS63wRd9KZBrzecDRQ2O9VrrGhGRm0VkjYisKS8vb1dQ9QkBIBAoa9c5lFK9Q1lZGddffz3Dhw/nrLPOYurUqbz88susWLGCrKwsJkyYwOjRo5k2bRqvv/56o2ODwSA5OTnceeedzZy960l0UjjPGDORSNXQHBGZ1nCjVSpo0yg/xpinjTGTjDGTcnKafEtbKaXiwhjDjBkzmDZtGrt372bt2rUsWbIEr9cLwPnnn8+6devYvn07CxcuZO7cuSxbtix6/DvvvMOoUaN44YUX6C4DmiU0KRhjiq3PA8DLwDlAmYgMBrA+D1i7FwMFDQ73WOvizunMbTCfnYhLKKV6gOXLl5OSksItt9wSXXfKKadw2223nbRvYWEhd911F48++mh03eLFi5k3bx5Dhgxh5cqVnRJzRyXsmYKI9AFsxphj1vwXgJ8CrwGzgfusz1etQ14D5orIEmAycDQRzxMAPve5Umpr97Fq1XAGDfpmIi6hlIqznTvnU1W1vtX9Kis/wJgAACJOMjPPbXbf9PRCRo58uNntmzdvZuLE2BtCTpw4kV/96ldApOXUP//5T5566imOHDnC4sWLOffc5mPpKhJZUsgF3heRDcCHwBvGmLeIJINLRWQncIm1DPAmsBvYBTwD3JrA2HC7h5CTM5OSkmcIBqsSeSmlVCeqTwgnzsfDnDlzOPPMMzn77LObufbxKqLXX3+dCy+8kNTUVGbOnMkrr7xCKBSKazyJkLCSgjFmN3BmE+srgIubWG+AOYmKpykez3zKy/9CaenzeDxzO/PSSqk2aukv+oZWrGjcDHPChBXtvua4ceN46aWXosuPPfYYBw8eZNKkSU3uv27dOsaMGQNEqo7ef/996rv3r6ioYPny5Vx66aXtjqcz9Oq+j7KyppKRMZni4kcwJpzscJRScdD4mWFuC3u27qKLLqK2tpYnnjj+2pTP52ty340bN3LPPfcwZ84cKisree+999i3bx979uxhz549PPbYYyxevLhD8XSGXt/Nhcczn61br6Oi4k2ys69MdjhKqQ6K53tHIsIrr7zCggUL+OUvf0lOTg59+vTh/vvvB+C9995jwoQJ+Hw+Bg4cyMKFC7n44otZtGgRF110ES6XK3qu6dOn84Mf/IC6urpG67sa6S7NpJoyadIk09FBdsLhAKtXDyc1dRSFhctaP0Ap1Wm2bt0arY5R7dPUz1BE1jboZaKRXl19BGCzOcnPn8uRI8upqtqY7HCUUiqpen1SABg8+DvYbGl4vbE9yFJKqZ5KkwLgdPZn0KDZlJX9Cb//QOsHKKVUD6VJwZKffzvG+CkpeTLZoSilVNJoUrD06XMa/ftfTnHx44TDdckORymlkkKTQgMezwICgTIOHFiS7FCUUiopNCk00K/fJaSljaOo6KFu06OhUiqx7HY7hYWF0em++yI98wSDQX784x8zcuTI6LZ77703ydF2XK9/ea0hEcHjmc+OHd/hyJF36dfvgmSHpJRKstTUVNavP7kjvv/5n/+htLSUTZs24Xa7OXbsGA8++GASIowvLSmcIDf3azgcA7R5qlLdVEUgwNMlJVQE4tsZXkM+n49nnnmG3/zmN7jdbgAyMjK4++67E3bNzqIlhRPY7ank5d3Cvn0/x+fbRVraiGSHpJQC5u/cyfqqlns09oVCfFxVhQFu3bGDienppNntze5fmJ7OwyNHtnjOmpoaCgsLo8s/+tGPGDNmDEOGDCEjI6NN99AdaEmhCfn5tyLioLj4N8kORSnVBr5wpGPL8AnLHVFffVQ/XXvttSft89xzz1FYWEhBQQFFRUVNnKX70JJCE1yuPAYOvJbS0t8xbNhPcTiykh2SUr1ea3/RA2z3+Th77VqCxuAQ4aXx4xmdlhb3WEaMGMG+ffs4duwYGRkZ3Hjjjdx4442MHz++W4yZ0BItKTTD45lPKFTF/v3PJjsUpVSMRqel8dmUKTw8YgSfTZmSkIQAkJaWxk033cTcuXOpra0FIBQK4ff7E3K9zqRJoRkZGWeRlXU+Xu9CwuFgssNRSsVogNPJzXl5DHA643K++mcK9dOdd94JwL333svgwYMZP348EyZM4Pzzz2f27Nnk5eXF5brJotVHLfB45rN580wqKl4lJ2dmssNRSiVBc9VBTqeT++67L/reQk+hJYUWZGdPx+0eqs1TlVK9hiaFFojYyc+/naNH36eysmOD+SilVHegSaEVgwffhN2eoaUFpVSvoEmhFQ5HJoMGfYvy8qXU1ZUkOxyllEooTQox8Hhuw5gQxcWPJTsUpZRKKE0KMUhNPZXs7OmUlDxFKORLdjhKKZUwmhRi5PHMJxisoKzsT8kORSnVibxeL9OnT2fkyJGceuqpzJs3L/qS2ocffsgFF1zAyJEjmThxIl/84hfZtGkTAHfffTf5+fkUFhYycuRIrr76arZs2dLo3OvXr0dEeOuttzr9vpqjSSFGWVnTSE+fgNf7sI61oFQvYYzh6quvZsaMGezcuZMdO3ZQVVXFf//3f1NWVsY111zDz3/+c3bu3MnHH3/Mj370Iz799NPo8QsWLGD9+vXs3LmTa6+9losuuojy8vLo9sWLF3PeeeexePHiZNxek/TltRjVj7WwbdtsDh9+h/79v5DskJRSTSj7cxnbbtzGac+dRu71uR061/Lly3G73dx4441AZMCdhx56iGHDhgEwe/Zszj333Oj+5513XrPnuvbaa3njjTf485//zLx58zDG8MILL/DOO+9w/vnnU1tbG+2GO5k0KbTBwIHX8umnP8DrfUiTglKdbOf8nVStb7nrbH+Zn5odNRCGrd/Yyp579pCSm9Ls/umF6Yx8uPmO9jZv3sxZZ53VaF1mZiZDhgzh008/Zfbs2W26h4kTJ7Jt2zYAPvjgA4YNG8app57KBRdcwBtvvMHMmcnvOUGrj9rAZnORnz+HQ4feorp6a7LDUUqdoGZ7zfF+s8PWcieZPHkyY8aMYd68ec3u07DqefHixcyaNQuAWbNmdZkqJC0ptFFe3n+yd++9FBcvZNSoJ5IdjlK9Rkt/0dcr+3MZ27+9nXBNGFuqjdG/Hd2hKqSxY8fy4osvNlpXWVnJvn37uOSSS/j444+ZPn06AKtXr+bFF1/k9ddfb/Z869atY9KkSYRCIV566SVeffVV7r33XowxVFRURLviTiYtKbRRSspAcnO/TmnpIgKBimSHo5RqIPf6XEb/djSSIh1OCAAXX3wxPp+P3//+90Ckc7w77riDb37zm9xxxx08//zzfPDBB9H9fb7mm6y/9NJL/OMf/+C6665j2bJlnHHGGRQVFbFnzx727t3LzJkzefnllzsUbzwkPCmIiF1E1onI69byMBFZLSK7RGSpiKRY613W8i5r+9BEx9ZeHs88wuEaSkqeSXYoSqkT5F6fy+frPt/hhACRBiYvv/wyL7zwAiNHjmTUqFG43W5+/vOfM2jQIJYuXcqPfvQjRowYwbnnnsuLL77I3Llzo8c/9NBD0Sapf/zjH1m+fDk5OTksXryYL3/5y42uNXPmzC5RhSSJbl4pIt8DJgGZxpgrReQvwF+NMUtE5ElggzHmCRG5FTjDGHOLiMwCvmyMOXncuwYmTZpk1qxJTkd169dfgs+3jSlTPsNmi0+/7UqpxrZu3cqYMWOSHUa31tTPUETWGmMmNbV/QksKIuIBvgj81loW4CKgvpJuETDDmp9uLWNtv9jav0sqKFiA319MefmLre+slFLdRKKrjx4GfsDx9gADgCPGmPqhzLxAvjWfDxQBWNuPWvs3IiI3i8gaEVnT8CWQzta//+Wkpo7C631IX2ZTSvUYCUsKInIlcMAYszae5zXGPG2MmWSMmZSTkxPPU7eJiA2PZx7Hjn1EZeXKpMWhlFLxlMiSwueAq0RkD7CESLXRI0BfEalvCusBiq35YqAAwNqeBXTp5j25uTfgcPTVsRaUUj1GwpKCMeZHxhiPMWYoMAtYboz5GvAv4CvWbrOBV63516xlrO3LTRevl3E40hk8+DuUl79Ebe3eZIejlFIdloz3FH4IfE9EdhF5ZvCstf5ZYIC1/nvAnUmIrc3y8+cCQnHxo8kORSmlOqxTkoIxZoUx5kprfrcx5hxjzAhjzFeNMXXW+lpreYS1fXdnxNZRbvcQcnJmUlLyDMFgy/2yKKW6n/T09Ga3zZ8/n/z8fMLhcHTdihUryMrKorCwkMLCQn760582e3xZWRnXX389w4cP56yzzmLq1KnRF9jqzzNhwgRGjx7NtGnTTnpbOhgMkpOTw513xu9vaH2jOQ48ngWEQkcpLX0+2aEopTpJOBzm5ZdfpqCggHfffbfRtvPPP5/169ezfv167rrrriaPN8YwY8YMpk2bxu7du1m7di1LlizB6/U2Os+6devYvn07CxcuZO7cuSxbtiy6/Z133mHUqFG88MILcWsFqUkhDrKyppCRMZni4kcwJtz6AUqphFm5En7xi8hnIq1YsYJx48bx3e9+t11vIi9fvpyUlBRuueWW6LpTTjmF2267rcn9CwsLueuuu3j00eNV1YsXL2bevHkMGTKElXG6Ye0QL048nvls3XodFRVvkp19ZbLDUarHmT8f1q9veZ+jR2HjRgiHwWaDM86ArKzm9y8shIfb2Xhw8eLFXHfddUyfPp0f//jHBAIBnM5I7wYrV67kzDPPJC8vjwceeIBx48addPzmzZuZOHFim645ceJEfvWrXwFQW1vLP//5T5566imOHDnC4sWLG43t0F5aUoiTnJyZuFwevN6Hkh2KUr3W0aORhACRz6NHE3Mdv9/Pm2++yYwZM8jMzGTy5Mm8/fbbQOSLe+/evWzYsIHbbruNGTNmtHK2iDlz5nDmmWdy9tlnN7tPwyqi119/nQsvvJDU1FRmzpzJK6+8QigU6tiNEUNJQUS+CrxljDkmIv8DTAR+Zoz5uMNX70FsNif5+XPZvftOqqo2kp5+RrJDUqpHieUv+pUr4eKLwe+HlBT4059g6tT4x/L2229z5MgRTj/9dCDSO2pqaipXXnklmZmZ0f2uuOIKbr31Vg4ePMjSpUt55plIJ5pvvvkm48aN46WXXoru+9hjj3Hw4EEmTWqySyIg0vV2fT9Gixcv5v3332fo0KEAVFRUsHz5ci699NIO3VssJYX/ZyWE84BLiDQd1YEEmjB48Hew2dL0ZTalkmTqVFi2DO65J/KZiIQAkS/k3/72t+zZs4c9e/bw2Wef8c477+Dz+SgtLY3+Rf/hhx8SDocZMGAAc+bMiT58zsvL46KLLqK2tpYnnjj+ddpS19sbN27knnvuYc6cOVRWVvLee++xb9++aAyPPfZYXHpZjeWZQn155IvA08aYN0TkZx2+cg/kdPZn0KDZ7N//LMOH30dKysBkh6RUrzN1anyTgc/nw+PxRJdvvfVW3nrrLZ588snouj59+nDeeefxt7/9jfLycp544gkcDgepqaksWbKEpvr2FBFeeeUVFixYwC9/+UtycnLo06cP999/f3Sf9957jwkTJuDz+Rg4cCALFy7k4osvZtGiRVx00UW4XK7ovtOnT+cHP/gBdXV1jda3VatdZ1vjIBQDlxKpOqoBPjTGnNnuq8ZJMrvObo7Pt50PPzyNoUN/wtChTTdFU0rFRrvO7rhEdJ19DfA28B/GmCNAf+C/OhpoT5WWNpr+/a+guPhxwuG6ZIejlFJt0mpSMMb4iPRPVC0iQwAnsC3RgXVnHs98AoEyDhxYkuxQlFKqTVpNCiJyG1AGvAO8YU3Nj0yt6NfvEtLSxlFUpGMtKKW6l1iqj+YBo40x44wxp1uTtrdsgYjg8cynunoDR4682/oBSinVRcSSFIqIjIKm2iA392s4ndnaPFUp1a3E0iR1N7BCRN4Aok9OjTG/TlhUPYDdnkpe3i3s3XsvPt8u0tJGJDskpZRqVSwlhX1EniekABkNJtWKvLzvIuKguPg3yQ5FKdUBr7zyCiLCtm3H29h8+OGHTJs2jdGjRzNhwgS+/e1vs2XLFjweT6OutCHSmd3q1au5++67eeCBBzo7/DaJpfXRT4wxPwEeBB5ssKxa4XLlMXDgtZSW/o5gUGvglOquFi9ezHnnnRd9Y7isrIyvfvWr3H///Wzfvp1169Zx2WWXMWDAAIYMGcJ7770XPXbbtm0cO3aMyZMnJyv8Noml9dF4EVkHbAY2i8haETm5yz/VJI9nPqFQFfv3P9v6zkqpDjnoO8jdK+6OTgd9Bzt8zqqqKt5//32effZZliyJNDN/7LHHmD17NlMbvDr9la98hdzcXK677rrofgBLlixh1qxZHY6js8TyTOFp4HvGmH8BiMgFwDNAx/to7QUyMs4iK+t8vN6F5Offjs2mvZUr1R7z35rP+tKW+84+VHOITQc2RZf/uvWv9E/t3+z+hYMKefiylhuDvPrqq1x22WWMGjWKAQMGsHbtWj755BNmz57d5P7XXHMNhYWF/OY3v8HhcLB06VJeeOGFFq/RlcTyTKFPfUKAyNCaQJ+ERdQDeTwLqKvbS0XFq8kORakerZ+7H26HG4BURyr93P06fM7FixdH/9KfNWtWq53O5ebmMn78eJYtW8b69etxOByMHz++w3F0lphaH4nI/wP+YC1/nUiLJBWj7OyrcLuH4vU+TE7OzGSHo1S31Npf9PVe3fYqM5bOYMlXlnDV6Ks6dM1Dhw6xfPlyNm3ahIgQCoUQEWbPns3atWuZPn16k8fVVyHVVyd1J7GUFL4F5AB/taYca52KkYid/PzbOXr0fSoru1YHfkr1NFeNvoqVN63kS6O+1OFzvfjii3zjG99g79697Nmzh6KiIoYNG8Yll1zCokWLWL16dXTfv/71r5SVlQFw9dVX8+abb7J06dJu9TwBYmt9dNgYc7sxZqI1zTPGHO6M4HqSwYNvwm7P0JfZlEowEWGKZ0qT3VW31eLFi/nyl7/caN3MmTNZsmQJS5Ys4fvf/z6jR49mzJgxvP3222RkRFrr9+3bl6lTp5Kbm8vw4cM7HEdnarbrbBF52BgzX0T+Bpy0kzGmY+WyOOiKXWe3ZOfO+ZSUPMaUKXtxufKSHY5SXZ52nd1xbe06u6VnCvXPELr2mxbdiMdzO8XFCykufozhw+9NdjhKKXWSZquPjDFrrdlCY8y7DSegsHPC61lSU4eTnT2dkpKnCIWaH3ZPKaWSJZYHzU01xv1mnOPoNTye+QSDFZSV/SnZoSjVLWj38+3Xnp9ds0lBRK6znicME5HXGkz/Ag51IM5eLStrGunpE/B6H9Z/7Eq1wu12U1FRof9X2sEYQ0VFBW63u03HtfRM4QNgP5BNpN+jeseAjW2OUAHHx1rYtm02hw+/Q//+X0h2SEp1WR6PB6/XS3l5ebJD6Zbcbjcej6dNxzTb+qg76G6tj+qFw3WsWjWU9PRCzjjj78kORynVy7TU+iiWDvGmiMhHIlIlIn4RCYlIZfzD7D1sNhd5ebdy6NBbVFdvTXY4SikVFcuD5keB64CdQCrwbeCxRAbVG+Tl/SciLoqLFyY7FKWUioolKWCM2QXYjTEhY8xzwGWtHSMibhH5UEQ2iMhmEfmJtX6YiKwWkV0islREUqz1Lmt5l7V9aPtvq+tLSRlIbu7XKS1dRCBQkexwlFIKiC0p+Kwv7vUi8ksRWRDjcXXARcaYM4m813CZiEwB7gceMsaMAA4DN1n73wQcttY/ZO3Xo3k88wiHaygpeSbZoSilFBDbl/s3rP3mAtVAAdBqV58mospadFqTAS4CXrTWLwJmWPPTrWWs7RdLPDov6cLS00+nX79LKC5+lHA4kOxwlFKq5aQgInbg58aYWmNMpTUU5/es6qRWiYhdRNYDB4iM8/wpcMQYE7R28QL51nw+UARgbT8KDGjinDeLyBoRWdMTmql5PPPx+4spL3+x9Z2VUirBWkwKxpgQcEp9vX9bWc8gCgEPcA5wWnvOc8I5nzbGTDLGTMrJyeno6ZKuf//LSU0dhdf7kL6go5RKupgG2QH+T0ReI1J9BIAx5texXsQYc8R6E3oq0FdEHFZpwAMUW7sVE6ma8oqIA8gCevwTWBEbHs88du6cQ2XlSrKydJRTpVTyxPJM4VPgdWvfjAZTi0QkR0T6WvOpwKXAVuBfwFes3WYD9WNUvsbxfpa+Aiw3veRP59zcG3A4+upYC0qppGu1pGCMqW9KmmaMaUvXnoOBRdZzCRvwF2PM6yKyBVgiIj8D1gHPWvs/C/xBRHYR6Vupew1X1AEORzqDB99MUdED1Nbuxe0+JdkhKaV6qVjeaJ5qfZFvs5bPFJHHWzvOGLPRGDPBGHOGMWa8Mean1vrdxphzjDEjjDFfNcbUWetrreUR1vZeNQ50fv4cQCgufjTZoSilerFYqo8eBv4Dq37fGLMBmJbIoHojt3sIOTkzKSl5hmCwqvUDlFIqAWJ9o7nohFWhBMTS63k8CwiFjlJa+nyyQ1FK9VKxJIUiETkXMCLiFJHvE3lgrOIsK2sKGRmTKS5+BGPCyQ5HKdULxZIUbgHmEHm5rJhIlxW3JjKo3qygYAE1NbuoqHgz2aEopXqhWJLCaGPM14wxucaYgcaYrwNjEh1Yb5WdfTUulwev96Fkh6KU6oViSQq/iXGdigObzUl+/lyOHFlOVZUOcKeU6lzNvqcgIlOBc4EcEfleg02ZgD3RgfVmgwd/hz17forX+zCnnfa7ZIejlOpFWioppADpRBJHwzeZKzn+RrJKAKezP4MGzaas7E/4/QeSHY5SqhdptqRgjHkXeFdEnjfG7O3EmBSRsRZKSp6gpORJhg69K9nhKKV6iZaqjx42xswHHhWRk/ogMsZcldDIerm0tNH0738FxcWPM2TID7HZXMkOSSnVC7TU99EfrM8HOiMQdTKPZz4bN36BAweWMGjQ7NYPUEqpDmqp+mit9flu54WjGurX7xLS0sbh9T5Mbu4N9PCB6JRSXUBM3Vyo5BARPJ75VFWt58gRzc1KqcTTpNDF5eZ+DaczW8daUEp1imaTgoj8wfqc13nhqBPZ7ank5d1CRcVr+HwxDY2tlFLt1lJJ4SwRyQO+JSL9RKR/w6mzAlSQl3crIg6Ki/VFcqVUYrWUFJ4ElgGnAWtPmNYkPjRVz+UazMCBsygt/R3B4NFkh6OU6sGaTQrGmIXGmDHA74wxw40xwxpMwzsxRkWkeWooVMX+/c+2vrNSSrVTqw+ajTHftYbgnGtNZ3RGYKqxjIyJZGVNw+tdSDgcTHY4SqkeKpYxmm8H/gQMtKY/ichtiQ5MnczjmU9d3V4qKl5NdihKqR4qliap3wYmG2PuMsbcBUwBvpPYsFRTsrOvwu0eps1TlVIJE0tSEBqPyRyy1qlOJmInP/92jh59n8pKfdavlIq/WJLCc8BqEblbRO4GVgH6tDNJBg/+FnZ7hpYWlFIJEcuD5l8DNwKHrOlGY4x+IyWJw5HJ4ME3UV6+lLq6kmSHo5TqYWLq5sIY87HVRHWhMWZdooNSLcvPvw1jQhQXP57sUJRSPYz2fdQNpaYOJzt7OiUlTxIK1SQ7HKVUD6JJoZvyeBYQDFZQVvbHZIeilOpBWkwKImIXkX91VjAqdllZ55OePgGv92GMOWlgPKWUapcWk4IxJgSERSSrk+JRMYqMtbAAn28Lhw+/k+xwlFI9RCzVR1XAJhF5VkQW1k+JDky1buDAa0hJGaTNU5VScdPSGM31/mpNqoux2Vzk5d3Knj13UV29lT59xiQ7JKVUNxfLewqLgL8Aq4wxi+qn1o4TkQIR+ZeIbBGRzfWD9VjjMbwjIjutz37WerFKIbtEZKOITOzozfUGeXm3IOKiuFgLb0qpjoulQ7wvAeuBt6zlQhF5LYZzB4E7jDFjifSXNEdExgJ3AsuMMSOJjNdwp7X/5cBIa7oZeKKN99IrpaTkkJv7dUpLFxEIVCQ7HKVUNxfLM4W7gXOAIwDGmPVAq+MpGGP2G2M+tuaPAVuBfGA6UF/SWATMsOanA783EauAviIyOPZb6b08nnmEwzWUlDyT7FCUUt1cLEkhYIw5cbivcFsuIiJDgQnAaiDXGLPf2lQK5Frz+UBRg8O81roTz3WziKwRkTXl5eVtCaPHSk8/nX79LqG4+FHC4UCyw1FKdWOxJIXNInI9YBeRkSLyG+CDWC8gIunAS8B8Y0xlw20m0sC+TY3sjTFPG2MmGWMm5eTktOXQHs3jmY/fX0x5+YvJDkUp1Y3FkhRuA8YBdcBioBKYH8vJRcRJJCH8yRhT34KprL5ayPo8YK0vBgoaHO6x1qkY9O9/Oampo/B6H9KX2ZRS7RZL6yOfMea/gTIBS5MAACAASURBVIuBC40x/22MqW3tOBERIl1sb7V6Wq33GjDbmp8NvNpg/Q1WK6QpwNEG1UyqFSI2PJ55HDv2EZWVK5MdjlKqm4ql9dHZIrIJ2EjkJbYNInJWDOf+HPAN4CIRWW9NVwD3AZeKyE7gEmsZ4E1gN7ALeAa4te2307vl5t6Aw9FXX2ZTSrVbLC+vPQvcaox5D0BEziMy8M4ZLR1kjHmf5kdou7iJ/Q0wJ4Z4VDMcjnQGD76ZoqIHqK3di9t9SrJDUkp1M7E8UwjVJwSIftkHExeS6oj8/LmAUFz8aLJDUUp1Q80mBRGZaL1V/K6IPCUiF4jI50XkcWBFp0Wo2sTtLiAn5yuUlDxDMFiV7HCUUt1MS9VHD56w/L8N5rV5Sxfm8cynvHwpZWWLyM/XGjmlVOyaTQrGmAs7MxAVP1lZU8jImIzX+wh5ed9FRMdSUkrFptUHzSLSF7gBGNpwf2PM7YkLS3VUQcECtmyZRUXFm2RnX5nscJRS3UQsrY/eBFYBm2hj9xYqebKzr8bl8uD1PqRJQSkVs1iSgtsY872ER6LiymZzkp9/G7t3/5Cqqo2kp7fYglgppYDYmqT+QUS+IyKDrbEQ+otI/4RHpjps8OBvY7Ol4fU+kuxQlFLdRCxJwQ/8ClgJrLWmNYkMSsWH09mfQYNmU1b2J/z+A60foJTq9WJJCncAI4wxQ40xw6yp1fEUVNfg8czDmDpKSp5MdihKqW4glqSwC/AlOhCVGGlpo+nf/wqKix8nHK5LdjhKqS4ulqRQDay33mpeWD8lOjAVPx7PAgKBMg4cWJLsUJRSXVwsrY9esaYeoyIQ4KXycmbm5DDA6Ux2OAnXr9/FpKWNw+t9mNzcG4j0aq6UUidrNSkYYxa1tk93st3nY8KaNdSFw8zftYuPzzqL0/r0SXZYCSUieDzz2bHjOxw58i79+l2Q7JCUUl1ULOMpfCYiu0+cOiO4RNjhizweCQM14TCXb9zI70tL8Yd79nt5ublfw+nM1rEWlFItiuWZwiTgbGs6H1gI/DGRQSXSqLQ0HCKk2my4bTYcIszeto1hq1Zx/759HAn0zIHv7fZU8vJuoaLiNXy+XckORynVRcUyHGdFg6nYGPMw8MVOiC0hRqel8dmUKTw8YgTeqVPZMXkyfz/9dMakpXHn7t0UrFrFgl272FNTk+xQ4y4v71ZEHBQX/ybZoSiluihpbZB3a0yFejYiJYfvGmPOTGRgsZg0aZJZsyZ+79GtP3aMB71elhw4QNgYvpqTwx0FBZydmRm3ayTb1q03cPDgy0yd6sXhyEp2OEqpJBCRtcaYSU1uiyEp/KvBYhDYAzxgjNketwjbKd5JoZ63tpaFxcU8VVJCZSjEtKws7igo4MoBA7B185Y7x459zNq1Z3HqqQ9SUKBdWinVG3UoKXRliUoK9SqDQZ7dv5+HvV721dUxKjWV7xUUcENuLql2e8Kum2jr1n2e2tq9TJ68C5stllbJSqmepKWkEEvrI5eIXC8iPxaRu+qn+IfZ9WQ6HCwoKODTyZNZPGYMGXY7t+zYwZBVq7j7s88o9/uTHWK7eDzzqavbS0XFa8kORSnVxcTS+uhVYDqRqqPqBlOv4bDZmJWby0dnncWKwkKmZGbyk717GbJqFbds3852X/fqBSQ7+yrc7mF4vQ8lOxSlVBcTS92BxxhzWcIj6QZEhM/37cvn+/ZlW3U1v/Z6eb60lKf27+eqAQO4o6CA87OyuvwbwyJ28vNv59NPF1BZuYbMzCZLkUqpXiiWksIHInJ6wiPpZk7r04enR49m39Sp3HXKKfzf0aN8fv16Jn/8MUsPHCDYxV+GGzz4W9jtGfoym1KqkViSwnnAWhHZLiIbRWSTiGxMdGDdxcCUFH4ybBj7pk7liZEjORIMMmvLFkasXs3DRUUcCwaTHWKTHI5MBg++ifLypdTVlSQ7HKVUFxFLUrgcGAl8AfgScKX1qRpIs9u5JT+fbeecwyvjx1PgdrPg008pWLmSH376KcV1Xa/b6vz82zAmRHHx48kORSnVRWiT1ARaXVnJg0VFvFRejk2E6wcO5I6CAs5IT092aFGffHI1R478m6lTi7DbU5MdjlKqE3SoSapqv8mZmfxl3Dh2TZ7MrXl5vFRezplr1vCFDRt4+9AhukJC9njmEwxWUFbWbbuzUkrFkSaFTjAsNZVHRo6kaOpUfjFsGJ9UV3PZxo2cuWYNz+/fT10SH0pnZZ1PevoEvN6Hu0SSUkollyaFTtTP6eTOU05hz5QpPH/aaRjgxu3bGbZqFb/Yu5fDSeihNTLWwgJ8vi0cPvxOp19fKdW1aFJIghSbjdmDBrFx0iTePuMMxvfpw48/+4yClSu5fedOdndyD60DB15LSsogbZ6qlEpcUhCR34nIARH5pMG6/iLyjojstD77WevFGvt5l9XsdWLzZ+45RIQv9O/PP848k/WTJjEzJ4cnS0oYuXo112zezOrKyk6Jw2ZLIS9vDocO/Z3q6q2dck2lVNeUyJLC88CJb0LfCSwzxowEllnLcLzZ60jgZuCJBMbVJZ2Zns6iMWP4bMoU/quggH8cOsSUjz/m/HXreKW8nFCC6/vz8v4TERfFxQsTeh2lVNeWsKRgjPk3cOiE1dOB+jGfFwEzGqz/vYlYBfQVkcGJiq0ry3e5uO/UUymaOjUyEFBdHV/evJkxH37IE8XF+EKhhFw3JSWH3NyvU1q6iECgIiHXUEp1fZ39TCHXGLPfmi8Fcq35fKCowX5ea91JRORmEVkjImvKy8vbFcTKlfCLX0Q+u6oMh4N5Hg87zzmHpWPH0tfh4NadOxmyciV3ffYZZQnoodXjmU84XENJyTNxP7dSqntI2oNmE2n/2OY6EWPM08aYScaYSTk5OW2+7sqVcMEF8OMfw7Rpkc81a6Crjr7psNm4ZuBAVk+cyL8LC/lcVhY/27uXU1au5Dvbt7O1On4d1qanj6dfv0soLn6UcLhnjlWtlGpZZyeFsvpqIevzgLW+GChosJ/HWhd3K1ZAfXdEwWCkxHD22ZCeDmPHwnXXwX33wd//DiUl0FWa7osI5/fty6unn87Wc87hm4MG8ceyMsZ+9BFXbtzIisOH4/KegcezAL+/mPLyF+MQtVKqu+nsYbdeA2YD91mfrzZYP1dElgCTgaMNqpni6oILwOUCvx9SUuD558Fmgw0bItMHH8CSJcf3z86GM89sPI0ZEzk2WUanpfHk6NHcM2wYj5eU8GhxMRdu2MBZ6encUVDAV3JycNral+/797+M1NRReL0PMXDgrC7fDbhSKr4S1veRiCwGLgCygTLgf4FXgL8AQ4C9wDXGmEMS+eZ5lEhrJR9wozGm1U6N2tv30cqVkRLDBRfA1Kknbz98GDZuPJ4oNmyATz6B+j7tnM5IYigsbJwssrPbHEpc1IRC/KGsjAeLithRU8MQl4t5Hg/fHjyYTEfb835x8ePs3DmHCRM+ICuriR+QUqpb0zGa4yAYhB07GieKDRtgf4PyTF7eyaWKkSOhHd/L7RI2hjcqKniwqIh3jx4l027n5rw8bs/Pp8Dtjvk8oVA1K1d66NfvUsaN+0sCI1ZKJYMmhQQqLz85UWzZcvy5hdsN48c3ThRnnAF9+yY2ro+sHlpfLC9HRLg2J4c7CgqYkJER0/GffvpDiooeYMqU3bjdpyQ2WKVUp9Kk0Mn8fti69XiSWL8+8lnRoPn/KaecXKoYPjzyfCOe9tTU8EhxMb/dv5+qUIiL+vbl+wUFXNa/f4vPC2pri1i1ahgFBQs49dRfxTcopVRSaVLoAoyJtGY6sVSxYwfUd5Kang6nn944UZx+emR9Rx0JBHh6/34e8Xop8fsZl5bG9woK+FpuLq5mMtHmzbM4dOgtpk714nB0nTEglFIdo0mhC/P5YPPmk5NFfbdHInDqqSeXKoYMiWxrK384zNIDB3igqIiN1dUMSklhbn4+t+TlMcDpbLTv0aOrWLduKiNHPkp+/pw43K1SqivQpNDNGAN7956cKD799Pg+fftGnk00TBTjxkFqjIOnGWNYdvgwDxQV8fbhw6TZbNw4aBALCgo4tcFJPv54KoFABeecsw0R7VRXqZ5Ak0IPcewYbNrUOFFs2gT1LzXb7TB69MmlikGDWi5VbKqq4tdeL38qKyNoDF/Ozub7BQVMzcriwIGlbNkyi/Hj/0Z29pWdc6NKqYTSpNCDhcOREsSJpYp9+47vk5NzcqI47bSTX8Arqavj0eJiniwp4XAwyNTMTL7nySN150WsclzBvAk/Izsl9qatSqmuSZPCCfYv2s+Om3dw2nOnkXt9busHdEOxvIA3duzJySI7G6qCQZ4rLeUhr5fPamsRwtgJYSfEXfyC//Bcwil9J5GVNQWns19yb1Qp1WaaFBoo+3MZ2765DRMwIJAxJYMBVwygz9g+pI1JI3VEKjZnz6w7b/gCXn0z2Q0boLT0+D4NX8A7/QzD8v5FPLujDDYOgMIjMC7yBNyJn2wOMtDmIz/FQUFqX4amFzA03YPH5SLP5SIvJQW33Z6ku1VKNUeTQgPvut7F+Ju/Z3EIqSNTSRubFk0Ufcb2IXV0KnZ3z/yCO3Dg5OqnrVuPv4AX7czWBgVnrWfo6cMI9DmEL/0Q1X2qqEwPcCTLTSDLBplBsB//+Q5wOMh3uY5PKSnkW0mjfj7b6dQ+lpTqRJoUGij7cxnbv72dcE0YW6qN0b8dTfb0bHzbfFRvrca3xUf1lmp8W33U7KoB6x0CbOAe5o4kirFp9BkT+Uw7LQ1HRmf3K5h4dXWRxHDvvfDii2EiHeoaMjKOYEw/qqqaPk4kTGpGNa6sozgzfdiy/Ji+NoJ9HdRmCdWZIcgKNJqc6SHy3U0njPopLyWFVC11KBUXmhROUPbnMrbduK3VZwrhujC+HT58W61EsSUy79vui1Q/WVwFrkaJor6E4ezvbPbc3cXKlXDxxcd7lV22LNKJYG1t5A3tgwePT+Xlx+cPHPBTVlZJeXmAigoHhw9nEgi4mryGzWFI6RvE3jdIODOAP6OOUJb/pOSR2T/M4IE2CgYKQ7JSTip91Jc6bFrqUKpFmhTiLBwMU/tp7UklC99WH+GacHQ/Z66zURVU2tg00sakkZKb0q2qS1rrVTYWoVCQsrIt7N27gX37dlJSUkR5eZAjR3KorMzF5zuNqqrhHDs2mCNHsjhY4eBQBYTDTf+cxB3CZAUgs3HisPUNkDXAkJMNg3KEgoE2hufaGTHIySnpkQSS53KRpqUO1YtpUugkJmyo3Vt7vGSx1RdNGqHK42MrO/o5mixZuApc3SpZdFRdXSmVlSuprPyAo0c/4NixtRgTaR7ldg+jT59zEbkAv/9cqqtHceiQo1GJ5EC5oaQ8TGm5oeKgcPSQUFPZQiOB9OPJw9E3RHq/EH2zDdlWAvEMtDEs186IQQ7G5jkZkZOCw957fh+q99CkkGTGGPz7/Y2qoOrnAwePD3tpT7eTdlpao0SRNjaN1GGpSC/4cgqH6zh2bF00SVRWfoDfH+mb3GbrQ2bmOWRmnktW1rlkZk7B6ex/0jn8/pOrtYrKQnxWGsJ7IEzpQcPBg3CkQqg+ZKP2sANT10wisRlsWQFc/UL06Rcma4AhO8eQmx1JIENz7Ywc5MCTayc7O/I+SFoavP1+kKfe9PGfV6Rx2fk973mT6v40KXRh/nJ/kyULf4k/uo+4hLTRjaug+oztE2k+m9Izm89CJJnW1e2LJoijRz+gqmoDECl1paWNITNzqpUkziUtbXSbu+IwJtL/VOmBMDtKg+zcH+CzshBFB8LsPxBJIIcrhKrDNuoO2wkfccJRJzRXreUMYwLWNoG80/wMHCikugS3G1JThDS3kOaGNJeQ5hLcbsHlijyzcbmOTw2X27LN4Whfv1iq99Ck0A0FjwYjzywaJArfVh+1n9VG9xGHkDoi9aSSRdroNOypPbPOPBSqprLyowaliZUEg4cAcDj6kZk5JVqayMg4J+69ux4LBvHW+tlR7mfH/gC7S4MUlRv2HwhTfhD2r8gguDETEMDA4Fro74eADQJifTYxH4pfchcx0UQR+ZQOJ5p4JCyns/lkFY/nVip2mhR6kFB1CN/2xlVQ1VurI81n6x9bSIPms2MaVEedloYjs2dVZxhjqKnZ0aA0sRKfb7O11UZ6+hkNqpzOxe0emtDnNkv+Vct1lzsjX/ROw/0vH2PauUJ1OIwvFKI6FKI6HKY6FIosW/NVgRDHag1VtWGO1Rqqa8L4/FBda/DVGmrqoLbOEPTLyUnFLxBsOtmI34YjZMcRtGEP2rEFbNiCNiQQmQjYMAEbJiCYgBD2CyFrCgaEUCC+P6umEkY4HOmWxZjIeCITJ0aq4uJRaoplW0pK7ytZaVLoBcJ1YXw7fSeVLHzbfY1e1nN5XI2qoOofeDsHdP/ms/UCgcNUVq5u8AB7NaFQ5MUKpzM3miCyss4lPX0idnt8+3N6673EPVMIhMP4mkgq1aFQdH2Ty60kpfopdOIFw1gJp+kk5AjYcIcduIJ2XKHIpzNkxxm0ppAde6CZhBSMnHPPJw4+22KnvnSVP8SQMxACdULAD/468PuFujqiUyBw0o+mQ5zOzktCsZbEUlKaH3SroyUrTQq9WDgYpnZ37UklC99WH2Ffg+azA52NSxZWy6iUQSkcWHwgpvc6uipjQlRXf9KoNFFbG+mHXCSFjIyJDUoTU3G58pIccfL4m0syMSSVlpJS/Xy4qYtuzoQ7zoyWrnhwQ7Q7lXpum63R5BIbrpAdV9BOSshOipWEHPWJKHg8GdmDtkhCspKRREtVkWRnAjaMP1JKCgePl5RCfiHojySmujrB76dRYmpq2e9v6gbbz+k8OWGEQuD1Rko3Ltfxd4faQpOCOokJG2r31R4vWTR45yJ09Pjfi5IqmDoT+YvRAf2v6E/WlCzs6fbIlGE/Pn/ilGZHbF2zXO73l3H0aMPmsGuizWFdrlMalSb69DkDm61nVbslgzGGunD4pATzj0OH+Mkbh6hbl0XKhKN899J0TktLoyYcprYd04nHBTr4HWfj5KTU7CQ2UkJ2nEEbDitJRSYbNitRNSwx0SBJGX8kQRG0nVSN1zDxrFsHmzYZQLDbDffcI/zoR227J00KKmbGGPyl/miC2DV/F03/eRcbWx8b9nQ7jgxH88mj4dRSkkm3Y+9jT0jz3HDYT1XVukYtnfz+ksg92NKaaA47IO4x9FbbfT7OXruWoDE4RPjorLMYnZYWt/OHrGTUnoQScyIKhZo9Xwf++wDgEGmUeAKb0im7fWy0ZLX47wFmXdi2KlBNCqrdTuwratQzo8iZmUOoKtT6dCyGfayppU4KT2RLtbWePNqSaNLt2ByNK28jzWGLokmisnIlx46to/5pfmrq6AaliamkpY3Rkek6oCIQ4KXycmbm5Jw0LGx3F4xzItpVU8P/fWAIre9L2sRKlnzVw5eys9sUkyYF1SGx9hXVEWF/OOYEEmuyCdfG/jeauKTV0oykhQmllOJ37KHOtoMathBOKYfUGuzpDtJzRpEx8AzCH57F/jsdjP7daAZ9bVBCfl6q94pHyUqTguqVwsE2JJoYk03Dh/OtM+AIQJofcQUhJYS4QojLYHMZxA02l2BzCza3PTKl2rG7ndhTndhTU7CnurCnuXGkuXGkpWJPdWFLtWFzxzA5tOTSU3W0ZNVSUuh1T88O+g7yyKpHALDb7Mw9Zy7ZaW0reqnuweawYetrw9k3ftURJmQI+ZpPIFu+tgWi41AIhO2kXVVEuDZMuNZg6iBcC6FaG8ZnwxyxQ50d/CkNJjv4TaSZCT5rage7QVxhbC6QVKzkY8PmtmN3O7C5HdhTndhSHdjd9taTTKzJyJra0sigM0qjPckAp5Ob8xLTSq7XJYU1JWu49717MdbAMb9e+Wv6p/YnPSW9ySkjJaPZbSft68ogxZ7SSgSqOxN7pJqpuTE0TNCw9VvroM4Nrlr4/q8452fLWjynMWHC4RpCoSqCwWOEQlWR+boqgr4qgj5fZKqpIeirJeSrI1RTPwUI1wYJ1QQJ14YaJB8DfifGn0KoPtnUuY4nnpoUONogEQVc4K/f7oQ6J5iOlTTEKTElF3+Zn2NrjkEItt6wlf3P7afP2D6IQ5qenM2sT9Bkc9rARpfqrDKRSbTXJYX/OPU/yM/Mx1vppZ+7H9eNv47qQDVV/qrodNB3MDp/zH8MXyD2v9ScNmdsCaQNycbtcHepf5CqebnX57J9x9OEf34LfP9XOC/f3OoxIjbs9j7Y7X1ISYnPf3BjDOFwbTTBRKZjJyxXEQodaWKdlZz8PisB+QnVBAj5/OB3nFCqaX4ygT6YYB/C/j6YQCrhQCr43RBwIbUpmMpIAgpsySDy4hoQgiP/PEzlR4ci7xCEwASBUPL//XdmImop8VVtrKL8L+WYgGH7t7cDxDUx9MpnCq9ue5UZS2fw6qxXuWr0Va3uHwqH8AV8jRJHw+mY/1iz21raL1Z2scdeWokx2aQ502JKNAd9B3n0w0ejy1rd1nsZYzDG33QCaSKxND2dvK955zx44L8ala64ZPkJFwfCNgjZ2zA5IeTCZlwQdiGhFGwmBUIuJJyChF0QdiIhJxJOgZATCTshlIKEnNY2B4QcELY+Qw4k7IicP2iH8AnXDNowIVukL6uQPfLWdsiGCQkEJfIZEkyw/pNGE9F5c3wKtPwdLSnC5+s+36bfpT5TOMFVo69i5U0rmZw/Oab97TY7Ga4MMlwZcYshbMLUBGo6lGT2V+0/aV3YxPYgVJCYkkdZdRkvbH4Bg0EQDtUconBQIQ6bI6GTTZt3djkigogLm80V1/c0wuf7+TeXwy9/GE0I55yzjXA4gDF+jAkQDkc+m5tvbbsxfut8DeerG8zXn6e5+ZPPGR27PEFEHIg4sYkTwYWEXUjYTfifnyP4wHWR6j5XLXLn00DbkkKL1+2NJYWeyhhDbbA2biWZKn8VlbWVdPz1m7aziS2xiUfif06n3Um1v5q/7fgbgiAiXD3mavq6+3b6z687OVJ7hCdWfJVwuBYBLh6UxXmT3o4kIQSb2No1bxNb9PfQ3vmWrgHhDialQLuSUkVNJc+8/DH833nwuff58vl7mH5p277Hu01JQUQuAx4B7MBvjTH3JTmkbkVESHWmkupMJadPTlzOaYzhpS0v8dUXv8qiGYu4ZPglBMPBuEyBUKB9x5rW9/GH/PgCvnbHFi/PrX8ubufqLf5edhQ2TEl2GDHrSGJpT4LzBXx47V6YtgcBxh6L7/10maQgInbgMeBSwAt8JCKvGWO2JDey3k1EmDl2ZrS6rTc88DbGEDbh2BJbuHFiC4QCXPfSdRQfKyY/I5/npj/XK35mHWGM4ZuvfpOSYyXkZeTx2y/9FohUsRpM9PfRnvmwCUeehXRgPmlxNHOOsAnzt+0vUBuCvFQ4d+DAuP4+ukxSAM4BdhljdgOIyBJgOqBJIclEhCme7vOXW0eJCHaxY7fZceFq8/GPXfEYM5bO4PEvPs6lp16agAh7nseveJwZS2fwxBef4PKRlyc7nC6vvrHM49Nf5bwYGsu0RZd5piAiXwEuM8Z821r+BjDZGDP3hP1uBm4GGDJkyFl79+7t9FiVaokxhtXFq3tNySoe9GfWNh39ebX0TKHbNfEwxjxtjJlkjJmUkxOfenOl4qm+ZKVfbrHTn1nbJPLn1ZWSQjFQ0GDZY61TSinVSbpSUvgIGCkiw0QkBZgFvJbkmJRSqlfpMg+ajTFBEZkLvE2kServjDGt9xGglFIqbrpMUgAwxrwJvJnsOJRSqrfqStVHSimlkkyTglJKqagu855Ce4hIOdCWFxWygYMJCqcr0/vufXrrvet9x+YUY0yTbfq7dVJoKxFZ09wLGz2Z3nfv01vvXe+747T6SCmlVJQmBaWUUlG9LSk8newAkkTvu/fprfeu991BveqZglJKqZb1tpKCUkqpFmhSUEopFdUjk4KIXCYi20Vkl4jc2cR2l4gstbavFpGhnR9l/MVw398TkS0islFElonIKcmIM95au+8G+80UESMiPaLJYiz3LSLXWL/zzSLy586OMRFi+Hc+RET+JSLrrH/rVyQjzngTkd+JyAER+aSZ7SIiC62fy0YRmdiuCxljetREpDO9T4HhQAqwARh7wj63Ak9a87OApcmOu5Pu+0IgzZr/bm+5b2u/DODfwCpgUrLj7qTf90hgHdDPWh6Y7Lg76b6fBr5rzY8F9iQ77jjd+zRgIvBJM9uvAP4OCDAFWN2e6/TEkkJ0WE9jjB+oH9azoenAImv+ReBi6f6je7R638aYfxljfNbiKiJjVnR3sfy+Ae4B7gdqOzO4BIrlvr8DPGaMOQxgjDnQyTEmQiz3bYBMaz4LKOnE+BLGGPNv4FALu0wHfm8iVgF9RWRwW6/TE5NCPlDUYNlrrWtyH2NMEDgKDOiU6BInlvtu6CYif1V0d63et1WMLjDGvNGZgSVYLL/vUcAoEfk/EVklIpd1WnSJE8t93w18XUS8RHpdvq1zQku6tn4HNKlLdZ2tOoeIfB2YBHw+2bEkmojYgF8D30xyKMngIFKFdAGRUuG/ReR0Y8yRpEaVeNcBzxtjHhSRqcAfRGS8MSac7MC6g55YUohlWM/oPiLiIFLErOiU6BInpuFMReQS4L+Bq4wxdZ0UWyK1dt8ZwHhghYjsIVLX+loPeNgcy+/bC7xmjAkYYz4DdhBJEt1ZLPd9E/AXAGPMSsBNpMO4ni4uQxr3xKQQy7CerwGzrfmvAMuN9aSmG2v1vkVkAvAUkYTQE+qXoZX7NsYcNcZkG2OGGmOGEnmWcpUxZk1ywo2bWP6dv0KklICIZBOpTtrdmUEmQCz3vQ+4GEBExhBJCuWdGmVyvAbcYLVCmgIcNcbsb+tJelz1kWlmWE8R+SmwxhjzGvAskSLlLiIPbmYlL+L4iPG+NM3IdgAAAVxJREFUfwWkAy9Yz9X3GWOuSlrQcRDjffc4Md7328AXRGQLEAL+yxjTrUvEMd73HcAzIrKAyEPnb/aAP/oQkcVEkny29bzkfwEngDHmSSLPT64AdgE+4MZ2XacH/KyUUkrFSU+sPlJKKdVOmhSUUkpFaVJQSikVpUlBKaVUlCYFpZRSUZoUlFJKRWlSUEopFaVJQak4EpEfisgfRORtEdkrIt9PdkxKtYUmBaXi63QgTOTN0kuAG5IbjlJt0+O6uVAqyU4HrjbGhEQkRMv93yvV5WhJQak4EREnkGOM+dRadQawKYkhKdVmmhSUip/RwNYGy4VEhotUqtvQpKBU/JwObGywrElBdTvaS6pSCWJ1zX66MaYm2bEoFSstKSiVACLSF/BrQlDdjZYUlFJKRWlJQSmlVJQmBaWUUlGaFJRSSkVpUlBKKRWlSUEppVSUJgWllFJRmhSUUkpF/X/r5bOgdD+F0AAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## For fixed rotation $\\eta=.05$ and varying CPU time ($x$-axis), compute relative error ($y$-axis)"
      ],
      "metadata": {
        "id": "R4VL4ldh3dlu"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import time\n",
        "\n",
        "aa = .05\n",
        "max_iters = 5000\n",
        "T, K = 10, 1\n",
        "last_K = max_iters - (T-1) * K\n",
        "n = (T-1) * K + last_K\n",
        "la_k, la_alpha = 5, .5  # hyperparameters for LA-GDA \n",
        "gamma = .1  # step-size for GDA, EG, OGDA, LA-GDA\n",
        "A1 = np.concatenate((aa*np.identity(dim), (1-aa)*np.identity(dim)), axis=1)\n",
        "A2 = np.concatenate((-(1-aa)*np.identity(dim), aa*np.identity(dim)), axis=1)\n",
        "A = np.concatenate((A1,A2), axis=0)\n",
        "\n",
        "# tm = [60, 90, 120, 150, 180, 210, 240, 270]\n",
        "tm = [60, 180, 270]\n",
        "eps_acvi, eps_gda, eps_eg, eps_ogda, eps_lak = [], [], [], [], []\n",
        "\n",
        "for dt in tm:\n",
        "  # ACVI ------------\n",
        "  t0 = time.time()\n",
        "  timesup = False\n",
        "  beta, mu, delta = .5, 1e-7, .5\n",
        "  y, lmd = init, np.zeros(((2*dim,1)))\n",
        "  for t in range(T):\n",
        "    mu *= delta\n",
        "    flag = 1\n",
        "    for k in range(K if t+1<T else last_K):\n",
        "      x = linalg.solve(np.identity(2*dim)+np.dot(Pc,A)/beta,np.dot(Pc,y-lmd/beta)+dc)\n",
        "      res = opt.minimize(lambda z: f_y(z,x+lmd/beta,mu), y, method='SLSQP', jac=None)\n",
        "      if not res.success:\n",
        "        print('Failed to update y. t={}. k={}'.format(t,k))\n",
        "        flag = 0\n",
        "        break\n",
        "      y = res.x.reshape(-1,1)\n",
        "      lmd += beta*(x-y)\n",
        "      if time.time()-t0 >= dt:\n",
        "        eps_acvi.append(np.linalg.norm(x-x_opt)/np.linalg.norm(x_opt))\n",
        "        timesup = True\n",
        "        break\n",
        "    if flag == 0 or timesup:\n",
        "      break\n",
        "  \n",
        "\n",
        "  # GDA ------------\n",
        "  z = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    z = proj(z - gamma*np.dot(A, z))\n",
        "    if time.time()-t0 >= dt:\n",
        "      # Relative error\n",
        "      eps_gda.append(np.linalg.norm(z-x_opt)/np.linalg.norm(x_opt))\n",
        "      break\n",
        "\n",
        "  # EG ------------\n",
        "  z = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    zz = proj(z - gamma*np.dot(A, z))\n",
        "    z = proj(z - gamma*np.dot(A, zz))\n",
        "    if time.time()-t0 >= dt:\n",
        "      eps_eg.append(np.linalg.norm(z-x_opt) / np.linalg.norm(x_opt))\n",
        "      break\n",
        "\n",
        "  # OGDA ------------\n",
        "  z0 = np.copy(init)\n",
        "  z1 = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    z = proj(z1 - 2*gamma*np.dot(A,z1) + gamma*np.dot(A,z0))\n",
        "    z0 = np.copy(z1) \n",
        "    z1 = np.copy(z0) \n",
        "    if time.time()- t0 >= dt:\n",
        "      eps_ogda.append(np.linalg.norm(z-x_opt) / np.linalg.norm(x_opt))\n",
        "      break\n",
        "\n",
        "  # LA5-GDA ------------\n",
        "  z = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    z_old = z\n",
        "    for _ in range(la_k):\n",
        "      z = proj(z - gamma*np.dot(A,z))\n",
        "    z = proj(la_alpha * z_old + (1-la_alpha) * z)\n",
        "    if time.time() - t0 >= dt:\n",
        "      eps_lak.append(np.linalg.norm(z-x_opt)/np.linalg.norm(x_opt))\n",
        "      break"
      ],
      "metadata": {
        "id": "bzcBqhO8xQry"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ms = 3\n",
        "print('CPU-time: {}'.format(tm))\n",
        "print('GDA:      {}'.format(eps_gda))\n",
        "print('EG:       {}'.format(eps_eg))\n",
        "print('OGDA:     {}'.format(eps_ogda))\n",
        "print('LA{}-GDA:  {}'.format(la_k, eps_lak))\n",
        "print('ACVI:     {}'.format(eps_acvi))\n",
        "\n",
        "plt.plot(tm,eps_gda,'ys-', label='GDA', markersize=ms)\n",
        "plt.plot(tm,eps_eg,'ch-', label='EG', markersize=ms)\n",
        "plt.plot(tm,eps_ogda,'mD-', label='OGDA', markersize=ms)\n",
        "plt.plot(tm,eps_lak,'bo-', label='LA{}-GDA'.format(la_k), markersize=ms)\n",
        "plt.plot(tm,eps_acvi,'gv-', label='ACVI', markersize=ms) \n",
        "plt.legend(loc='upper right')\n",
        "plt.xlabel('cpu time (s)')\n",
        "plt.ylabel('relative error')\n",
        "plt.yscale('log')\n",
        "plt.savefig('acvi_hd_bilinear_relErr_vs_time.pdf', bbox_inches='tight')\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "16AOAFRIh5ne",
        "outputId": "f2b70acc-d81f-4204-8807-9075cf57122c"
      },
      "execution_count": 24,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU-time: [60, 180, 270]\n",
            "GDA:      [0.5145986619628568, 0.43888856687104627, 0.3961038745965842]\n",
            "EG:       [0.3012852130692706, 0.08117955504120694, 0.03117351039238496]\n",
            "OGDA:     [0.5848684785035604, 0.5848684785035604, 0.5848684785035604]\n",
            "LA5-GDA:  [0.28825187323310963, 0.07110590692614277, 0.026617190845794735]\n",
            "ACVI:     [0.02904948188322766, 0.00016789817370326595, 4.061120650482063e-06]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAELCAYAAADHksFtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeVxUZfvH8c/FjoJLgoQii7vmlmFqLpktmmVmWqn1pD6ZlZlaWWmLWWZaVj9NrcwWs0fRbHNJM5fc0swlcsUNQdHcMFdknfv3xwyICIgyw8BwvV+v82LmcObMxXHky33u+9xHjDEopZRSeXFzdgFKKaWKNw0KpZRS+dKgUEoplS8NCqWUUvnSoFBKKZUvDQqllFL50qBQSimVLw0KpZRS+fJwdgFXIiJlgY+BVGCFMWaGk0tSSqlSxSktChH5UkSOici2HOs7isguEdkrIsNsqx8AvjPGPAHcV+TFKqVUKeesFsU0YBIwPXOFiLgDk4E7gQRgg4jMA0KArbbNMgqy84CAABMeHm7HcpVSyvVt2rTphDEmMOd6pwSFMWaViITnWH0zsNcYEwsgIrOALlhDIwSIJp8WkIj0B/oDhIaGsnHjRvsXrpRSLkxE4nNbX5w6s6sCB7M9T7Ct+wHoJiKfAPPzerEx5jNjTKQxJjIw8LJAVEopdY2KfWe2MeY80Lcg24pIZ6BzzZo1HVuUUkqVIsWpRXEIqJbteYhtXYEZY+YbY/qXL1/eroUppVRpVpxaFBuAWiISgTUgegC9rmYH2qJQyvWlpaWRkJBAcnKys0spsXx8fAgJCcHT07NA2zslKEQkCmgHBIhIAvCGMeYLERkILAbcgS+NMduvZr/GmPnA/MjIyCfsXbNSqnhISEjA39+f8PBwRMTZ5ZQ4xhgSExNJSEggIiKiQK9x1qinnnmsXwgsvNb9aotCKdeXnJysIVEIIkKlSpU4fvx4gV9TnPooCq2wfRRHZx5lpfdKjs48aufKlHKc0vi51ZC4XFpiGmc3nSUtMe2K217t8StOfRROdXTmUXb23gnpsPM/O4kbFYdXkJezy1IqX6lHU7mw+wJYSs/n1u0NN5LckpxdRrFiSbNgkg0AyfHWvhvPSgXrfygIl2pRiEhnEfns9OnTV/3amL4xkG57YoELuy7YtzilHODCLmtIAPq5LWJHTxylzwt9uOGOG2j1QCtue/g25i2Zx6r1qwiODKZl15Y06diEux69i0W/Lbrktenp6YS1DOP1D163Sy2ZIQGABZLj7NvR71ItisJ0Ztf9qi47/xsNKd7gnYy8NJmqA3sj4lGAxf0qt3OpfFZOdHTmUXb124XlggU3XzfqfF6HoF5Bzi7LoXbu3EmZOmWcWoMxhl59etG7d2++XfAtAPHx8cybN4/rQ6+nza1tWLBgAQDR0dHcf//9VKhVgdtvvx2ARYsWUbtebX5a9hPvT3m/0KfS0hLTrC0JC+AGPmE+hdpfTi4VFIUR1CuInTveg/dehqHjMO2Xs2PHAge9m9gleOBqA8oe29mrNnc9z2wHmaEQ0zemVIREcbF8+XK8vLx46qmnstaFhYXx7LPPsmLFiku2bdKkCSNGjGDSpElZQREVFcXgwYP55JNPWLduHbfcckuh6sk8zZQcl4xPmI9dTzuBiwVFoUc93bHcuthERm7BmHSMybB9zX+Bgm1X0P3lt53FknTV+8teX/HgbucQK3khC26FDsygXkGlNiD27BnCuXPRV9zuzJm1GGPt5BXxpFy5vH8x+/k1oVat8fnub/v27TRt2rTAdTZt2pRx48YB1lFbS5cuZcqUKZw6dYqoqKhCBwVYw8LeAZHJpYKisNdReHoGkZZ2NOuxn19De5ZXbBhjAIvDQuxatrvWkLVYkq/5fS+e3HeuwobYmTN/YEyqbV9eVKrUiUtD+OLj/Ndfus46obMj1l9NTZnrCxemmSGR87G9PPPMM6xZswYvL6+sQLj0/S/2ISxYsIDbbrsNX19funXrxqhRoxg/fjzu7u52r8teXCooCqtVqyPOLqFIWP/Tudv+A3s7uxynMcbYOcicE7KZIWH9mVK5cCE26/UX687Itr+MbD9TRrZti0dw5k6yQsPf/yfOnk0GhODgPoDYFrIeW3Pl4voNG264ZG916nyW7/YXLuzPdX3m41q1KjNnzkxSUg4DwocfvsaJEydp1aojaWknsVhSSE09nvWajRtXU6dODdLSTjJjxtesXbue8PBQABITE1myZD533HF7tp/10ve72vX2Pq2rQaFKLRGx/eVasv8brFhx6S+FZs3+vqb9ZAZnzgDJK1gKtv7ykCrs+vPn/fH0rJRZNWCw/sFusi3Y1lsfe3gEkJ5+Iuuxdb+Xv/7iX/657ydzfatWNblw4Rwff/wR/fp1B+D06SMYk0Fa2nEslmRSUqwzdm/btofRo99l4sRXOXZsC2vW/M7OnQvw9rYOY/7f/+YxY8aXtGoVck3/brkR8cTPr7Hd9ley/4fkoFdmq9Io5ynTa1VSgnPnzp34+IRe1Wtaty74VcgFNX/+Ep577jk++uhBAgMDKVu2DO+++z4+PkGsW7eFtm37kZSUROXKgUyY8BGdOt3D9Onf0L797VSs2ITM0OnWrTIjRjTFzS0Mb29vcg8s63Pr+svDzBhDampCVm32Pr0m2c+duYrIyEijNy5SyjXt3LmTevXqObuMYufs2Ut/5/n7R+a7fW7HUUQ2GWMue6EO6M8mMS2Nzw4fJjHN/p1dSinlSCKeuT62h+LdxixCu5KSaDRtF+nR5Xm2yTbmPFSN2ypUwN9DD5FSqvizZ59ETi71W7AwfRTzVqaSOqgxZAip7oYuKw5B08OUr5ZGRDhUL+9NuI8P4T4+RNi+hvn4aJAopVyeS/2WK8x1FP/8WQYstiFnGcCcajCnGqeBaDFsC0zFEnwBS/AFqHIWgo9BlWQqVksl4noPInx9LgaJr681SLy98dMgUUqVcPpbzObBu7z45F1DaqrBywvmfAvXXQexsbBvnxAb601srBd7/i7P0V8uDkf8FzhTJoPtVVJIC06yBklwIlRJhioXqBSSQYT/xdZIzhZJ2WJ8kY1SSoEGRZaWLWH5MmHFCmjXzvoc4NIr660BkZQEcXGZIQKxse7ExpZh3z5f9m+A5OSLQXLSzXA+KI0dVS6QHGQLkiqnocoFCL5AQCW5pDWSPUzCfHwoo0GilHIyDYpsWra8GBD5KVMG6te3LpcSLBY4ciR7iAixsV7s2+dF7KZyHD166cVRp/0yiAlJYWdwMklB57EEn4cqJyD4AgSlUNnX47IQyR4kvhokShU5d3d3Gja8OMVPjx49GDZsGOnp6YwYMYI5c+ZQtmxZAB588EFeffVVZ5VqFxoUdubmBlWqWJfWrXN+Vzh/HvbvzwwRa2tk374yxMaWYf/vFUlNvRgkbu6GC8Fp7K2aQsz1FzgfdI6M4FPWEKmSDP7pBHl65h4kvr6EentrkCjlAL6+vkRHXz4Z4WuvvcaRI0fYunUrPj4+nD17lg8++MAJFdqXBkURK1sWGjSwLjlZLMLhw9lDRKwtkVgvYtf6c+Z45Uu29y2fQUZIGvurJLMrOImzlc+Scf1xa4hUTgF3w/VeXnm2SEK9vfHRIFGlQGJaGt8fP063wEAqeTpmhtWkpCSmTp1KXFwcPj7W+0H4+/szcuRIh7xfUXKpoCjpU3i4uUFIiHW59dbLv3/2bGaAZG+NuBMb60PcigpkZLtO0N3DUL5qOm5VUzlY5QK7g85zuvI5MoJPWvtHylrnugnOL0h8fPB202syVfE1ZM8eos+dy3ebpIwMNp87hwEG7N5NUz+/fPv+mvj5Mb5WrXz3eeHCBZo0aZL1fPjw4dSrV4/Q0FD8/f2v6mcoCVwqKAo7zXhx5+8PjRtbl5wyMiAhIXuICPv2eRIb60nsb2U5mRhwyfZ+FS2Ur5aGV0gKh4OT2FP5PCcrn8USfBQCUsD2/6hKPkFSTYNElQBJFuusuBasH+ski6XQg0RyO/W0ZcuWS55/9dVXTJgwgcTERNauXUu1atUK9Z7O5FJBUZq5u0NYmHW57bbLv3/69KWtkX373KxDfnd5c/DXcqRnu5eRh6ehUrUM/Kul4lElhaPXJ7G38jlOBCZiCU6GMtbWiJBPkPj6Us3bGy8NEuVAV/rLH6yzLjTbtIl0Y/AQ4fsGDahTxv63Uq1ZsyYHDhzg7Nmz+Pv707dvX/r27UuDBg3IyMi48g6KMQ2KUqJ8ebjxRuuSU3o6HDyYPUSE2FgPYmM92LeoDKdOVbxk+woBFiqFpVMmJBUJTubE9Unsq3yWYwFHsVRKyZpBTICq3t55t0i8vfHUIFEOVqdMGfa3aOHwPooyZcrw+OOPM3DgQKZMmYKPjw8ZGRmkpqZe+cXFnAaFwsMDIiKsy+23X/79f//NHiIQG+tmHfL7txcH5vthyXa/Gy9vw/VhFipUS8O7agomOJl/rz/HvoCzHKl0HON7cWM38g+SEA0SZSeVPD3pX6WK3faXs4+iY8eOjB07ltGjR/P666/ToEED/P398fX1pXfv3lSx43s7g04zrgolLQ0OHMgeIhcf79tn7YDPLjDIUDksg3LVUvGsmkJG8AXOVT7H8YAz/ON/DpPtMhM3ICSPIAm3tUg8NEhKHZ1m3D6uZppxbVGoQvH0hBo1rMudd176PWPg5MmcIWI9rbVvgwcHfyiDMRdPa/n6GkLCDZVC0/GrloZ7lWRSr0/idOVzLKt0hsMcJfufNe7kHyQhGiRK2YUGhXIYEahUybo0a3b591NTIT4+Z4hYrx3Z9rsX586VBSplbV+liqFKuIXrwtLxrZoKtmlR/q18lqVlTnE4NeWyIKmWR4iE+/hQ1ctLg0SpAij2QSEi1YFXgfLGmO7OrkfZj5cX1KplXXIyBk6cuPR0ljVE3Nmx2p1Dh7wx5uJ49TJloH51Q3C4bdhv1WQs118g6frznKh0liVJJzmcmnpJkHiIUC2fFklVb2/c7XyTeqVKIocGhYh8CdwLHDPGNMi2viMwAesffZ8bY8bmtQ9jTCzwuIh858haVfEiAoGB1qVFi8u/n5x8cWLGS0drubN2uTtJST5Ahax9Va0KrWsYKodm4B+SinvVFNKDkzhb+Rz/ZJxn8UlrkGTnIUJoPkFSRYNElRKOblFMAyYB0zNXiIg7MBm4E0gANojIPKyhMSbH6/9rjDnm4BpVCeTjA3XrWpecjIFjxy7vXI+NFdYt8+DwYQ+gDGDtH/H3h+rVITLCEBBqHfbrVjWZlKDz/Fv+HActySw6eZJ/cgSJZwGCxE2DRLkAhwaFMWaViITnWH0zsNfWUkBEZgFdjDFjsLY+lCoUEQgKsi6XThNvlX2a+Owd7Xt3C7/+4klysidg7R9xc4Nq1aBudegQYeG60HS8Q1IxwRe4EHSeIz5JxKck8/PJkxzJJUjC8ukjCfby0iBRJYIz+iiqAgezPU8Amue1sYhUAkYDN4rIcFug5LZdf6A/QGhoqP2qVS4n72niuWSa+JxDfn9Z6MaRI16AF+AHBFK+vLU10qo6hEVYKBeajleVZNKvT+JMpfMczEgmLjmZ+SdOcDQt7ZL38rpCkFyvQVKsJSQk8Mwzz7Bjxw4sFgv33nsv48aNw8vLiz///JOXXnqJQ4cO4e/vT3BwMGPHjqVhw4aMHDmSqVOnEhgYyPnz52nYsCFvv/029bN9IKOjo7nxxhtZtGgRHTt2dOJPaVXsO7ONMYnAUwXY7jMR+Qfo7OXldZPjK1OuKP9p4smaJj5niGzfDgsWuJGSkhkk5XB3h9BQa5B0qQEh4RZr/0iVZFKDkzjmfYG4ZGuQzD1xgmM5gsS7AEEiIkUyM6q6lDGGBx54gKeffpq5c+eSkZFB//79efXVVxk6dCgPPfQQM2fO5BZbk3bNmjXs27cv6x4Wzz33HEOHDgVg9uzZtG/fnq1btxIYGAhAVFQUrVu3JioqqtQGxSEg++xYIbZ1hebqkwIq58t/mng4fPjyEImNhR9/hOPH3QAf21KBihWtIVKjBrSvDlXDLZSpmoJb1WTOXpfEwfTkrCD568QJjucIEh83N4I9PTmQkgLAoD17eDsigpblyxPm7U2wdrZnOTrzKDF9Y6j7VV2CegUVen/Lly/Hx8eHvn37AtYbGf3f//0fERERAPTu3TsrJABa5/ZXh83DDz/Mzz//zMyZMxk8eDDGGObMmcOSJUto06YNycnJWdOWO4szgmIDUEtEIrAGRA+glz12XNKnGVclW/Zp4tu2vfz7Z8/mvGmVdfnrL2uQpKW5Ab6ALx4eFQkLs4bITdXhwepQNSID76rW/pGjntbWyNrTp4nb7I2JrkBGk1O8aO36A6x9JNW8vQmz3Q0x3MeHsGzPXWGurT1D9nAuOv9pxlOPpnJh9wWwwM7/7CRuVBxeQV55bu/XxI9a4/OfbHD79u3cdNOlJy7KlStHaGgo+/bto3fv3gX/IYCmTZsSExMDwNq1a4mIiKBGjRq0a9eOn3/+mW7dul3V/uzN0cNjo4B2QICIJABvGGO+EJGBwGKsI52+NMZst8f7aYtCFWf+/tCokXXJKSMDDh3KfSqUOXMgMRGs/12sQVKpkjVE3MumY1a5gUXAwzD8k3O06JTGYUsy8SkpxNtaJL/aRm1lv47EDaji7Z0VHuG2AMl87Cp3SLyw6wJZP7jF+jy/oLC35s2bc+bMGe666y4mTJiQ6zbZp1KKioqiR48egPUWq9OnT3ftoDDG9Mxj/UJgoSPfW6mSJLM/IzQU2rW7/Ps5p4nPDJFNmzwgwwACacKYfuWyrhsJD7cut4dbJ3ysEmrBKziFtMBkDmUkZ4VIfHIya8+cYfaxY+ScDLuyp+elIZIjVMp5OLeb80p/+YP1tNOufruwXLDg5utGnc/rFPr0U/369fnuu0sv7Tpz5gwHDhzgjjvuYPPmzXTp0gWA9evX891337FgwYI89/fXX38RGRlJRkYG33//PXPnzmX06NEYY0hMTMyautxZin1n9tXQU0/KVeU1Tfy6dXD77UJqqnUW4BdftH7dv986BHj1apg5E9sMv9ZTW25uvpcESVvb12qhFnyrppEWcMEaJCkpWUHy97lzzDtxgpQck4hW8PAgzHYtSViOFkmYtzeVPD0RJ/eTZIZCTN8Yu4QEwO23386wYcOYPn06jz32GBkZGbzwwgv06dOHF154gebNm9OhQ4esfoqkpKQ89/X999/z66+/8sEHH7Bs2TIaNWrE4sWLs77fu3dvfvzxRx577LFC132tdPZYpUq4detgxQprS6Rly8u/n5ZmPa0VF3f5sn+/9c6I2aeKz+xryQySzCU0zOBfJZW0wJTLWiSZoXIuxw16yri5XdYiyR4q1zIEuLjMHnvw4EEGDBhATEwMFouFTp068f777+Pt7c0ff/zByy+/zKFDh6hcuTIBAQGMGDGCyMjIy4bHNmjQgNGjR1O/fn369u1L8+bNeeqpiwM9582bxyeffMKiRYvsWv/VzB7rUkGRrUXxxJ49e5xdjlIlQlqaNSxyC5K4uLyDJCLi0iAJCzNUDEknLSCZQ7YRW5n9JJmhcjL7rRSxXksSmqOTPXuo5DYDcHEJipKu1E4zrp3ZSl09T8+LN67KTVqa9Q6IuYXI8uXWILH+vSmAJ+7unoSE+GcFSKNw6GILlcBq6aQFpFwMkmwtkoW5XN2eeU+SsGzhcV9GBmfS0/ESwcvNTS9KLAIuFRRKKfvz9LRe71G9eu7fT029tEWS2T8SFwfLlllPe108ceGBu7sH1aqVzQqSeuFwt+3x9aEZWLIHSbYWyapTp0hISaF1+fLsvnDhYn0iWaHhbfua/bFeS1J4LhUU2pmtVNHz8rpykOTWItm/H5YssV6keDFI3PHwKEO1amWygqRmONxhexwSbuH8+RjCfX1JMYZUi4VUi4UUY0jKyOCUMeQ8me4Bl4VH9mBxF3F6h3tx51JBoaeelCp+vLwu3gUxNykpeZ/a+vXXnEHixqJFAhYPvLzA29u6f3/bVy8vg3gYUo2FVGNIsVi/plosJFssnLFYsOR4fzfAO0d4ZA8WDw0S1woKpVTJ4+0NNWtal9ykpFjvy54ZHuXLW6dSSU21Xl9y6cwm1l/qXl5ulwSJn+2rt3dmkNhaI5lhYnt8Li3tsmtJBGunu7ctPLIeZ2uduHqQaFAopYo1b+9L74S4c+elp7ksFmtopKRc/jXvIAEvL/esICmbFSTgltkisVzeKjmdnk5aLiNFcwuP7M9Leoe7SwWF9lEoVfq4uVlvZJXXvHlXFyTYgsQdb2/3rPAok6114u5hKF/On0OnTl0WJmczMhjzwgss++knFuzciZttaG/0mjU816MH1cLCEBHuvf9+Xnn99awgyd7hfvToUZ577jn++OMPKlasiJeXFy+99BJdu3ZlxYoVdOnSherVq5OUlERQUBAvvfQS99578VY+6enpBAcH8/jjjzN2bJ43D70qLhUU2kehlMrJEUFiscDR2Iv9JOVtXz08LKz9+WciQkM5vmEDzdu2JdUY9rq5cdMttzBxzhxSbR3ue7KN3PKwjdzyFKHHfffx8H/+w8dff42Xmxv/HDjAz/PnZ23bpk2brOlAoqOjuf/++/H19eX2228HYMmSJdSuXZs5c+YwZswYu5wWc6mgUEqp3OR39fqVgiQjwxoc2UNExLo+Z5Bs2LCCkJAbuPvuh5k+9VtuqHkH3t7gm+yNj3GnjpcfHp6GdNuprJQcfSUrly8HT0/a9+7N3uRk604DAmjTty/bz58nITmZCxYLR1JS8HJzo1bDhrz62mtMmjSJW2+7jX/T05lhm678k08+Yd26dZdMd36tNChsTiSdYNKfk7KeD7x5IAFlApxYkVLqSoYMgejo/Lc5fRq2bLG2HNzcrLP3li+f9/ZNmsD48Refu7uDr691ySQCmRc1Zw+S8eOjeOihntx6axcmTnyF48fTEPHkyBH44491NGjQmMDAKrz44vvUr3/Dxc52L6jkDWd2x9ImMpLGZcte1j+SarGQZnuckO3CxDL16hE9bhzR58+TkpzM4qVLmfjJJ5w6dYqoqCgNCnvaeHgjb618C4PBTdxoEdKCjjWdf2cppVThnD59cQoSi8X6PL+guFqZQeLunspvvy3kk08+xN/fn1atmnP06GLuvvteqlVryr33xuPp6ccvvyxkyJD7WbRoD//+C9lnNTl8GI4dE3btsI7aevvtQWzatAYvLy9Wr95ANU9fynt40MQWJKkWC6c8PXETQYA1v/xCZJs2iLc33bp1Y9SoUYwfPx73Qk4X71JBUZjO7A41OuDv7c+ZlDMIwpYjW2gd2ho/Lz/7F6qUsovsf/nnxTrDrvUvfi8vmDEj98kTC2vx4sWcOnUq63anSUlJ+Pr6cu+99xIUVC5ru8ce68SIEQMIDDzB7NmzmTp1KhYLfPvtQlq0uIH33/8eX1/raa7nn5/MiRMneOyxSGJirFPLnz4NMTsyh/+6s2XNNurWrAenPVk84wf+3vQ7jWvVQoDExESWL1/OnXfeWaifrWTf3ioHY8x8Y0z/8tfw54KIMP3+6QA0DGrIy8teJmJCBGPXjOVsyll7l6qUKiItW1qnEhk1yvrVESEB1hsOff7558TFxREXF8f+/ftZsmQJSUlJHDlyJOvmRH/++ScWi4VKlSrxzDPPEB0dzZYt0dStW4Vu3dpjTDK//voJ9etbT4PVqJGEp6f1OpPKla1h5+trbYmsX7+FDz4YxT2dBnJ2byrRf/7OgnkH2L7NWsPkyZOJiooq/A9njHG55aabbjLXwmKxmHUH12V97TSjk2Ek5rp3rzOjVo4ypy6cuqb9KqXsZ8eOHc4uwYiIqVq1atYyevRoU7FiRXP69OlLtuvatauZNWuWmThxoqlfv75p1KiRad68ufn999/z3Pfhw4fNww8/bMLDw02zZs1Mu3btzKxZs4wxxvz222+mXLlypkmTJqZ27dqmdevW5scf55n4eGPeeGOaufPOh82GDcYcPmzdV2JiogkICDDJycmXvU9uxxHYaHL5nepS04xnsuf9KDL7Lubvnk957/IMaTGEwc0HU9G3ol32r5S6OjrN+OXOnYPduy922NeuDX5XOGt+NdOMu9SpJ0eIrBLJvJ7z2Nx/M+0j2vPmyjcJnxDO68tfJzEp0dnlKaUUfn7WcKhatWAhcbU0KAroxuAb+eHhH4h+Mpq7atzF26vfJnxCOMOXDuf4+ePOLk8pVcr5+UFwsP1DAlwsKESks4h8dvr0aYe9R+PrGzPnwTlsfXor99a+l3d/f5eICRG8tOQljp476rD3VUopZ3GpoDCFGPV0tRpUbkBUtyh2PLODrvW68sG6D4iYEMHzi5/nn7P/OPz9lVKqqLhUUDhD3YC6fNP1G3Y+s5OHbniIj9Z/RPWPqjN40WAOnTnk7PKUUqrQNCjspHal2ky7fxq7Bu7ikYaP8PHGj6n+UXWe+fkZDpw+4OzylFLqmmlQ2FmN62rw+X2fs3vgbvo07sPUzVOp+VFNnpz/JHGn4pxdnlLKjn766SdEhJiYmKx1f/75J23btqVOnTrceOON9OvXjx07dhASEoLFcun99Zo0acL69esZOXIk77//flGXX2AaFA4SUTGCKZ2nsHfQXvo17ce0v6dRa2It+s3rR+y/sc4uTyllB1FRUbRu3Trr6uejR4/y4IMP8u6777Jr1y7++usvOnbsSKVKlQgNDWX16tVZr42JieHs2bM0b97cWeUXmAaFg4WWD+Xjez5m36B9PB35NP/b8j9qT6xN37l92ZO4x9nlKeXyTiSdYOSKkVnLiaQTdtnvuXPnWLNmDV988QWzZs0CYPLkyfTu3ZuW2eYJ6d69O0FBQfTs2TNrO4BZs2bRo0cPu9TiaC41KWBxFlIuhI/u/ojhrYczbu04Pt34KdP/nk6vhr14tc2r1A2o6+wSlSpxhvwyhOgj+c8zfvLCSbYe25r1/IedP3Cd73V5bt/k+iaM73jl2Qbnzp1Lx44dqV27NpUqVWLTpk1s27aN3r1757r9Qw89RJMmTZg4cSIeHh7Mnj2bOXPmXPF9igNtURSxYP9gPuzwIfsH79T6dBQAACAASURBVOf5Fs/zw84fqD+5Pr2+78WO4zucXZ5SLqeiT0V8PKx3JfL18KWij32m34mKispqEfTo0eOKk+8FBQXRoEEDli1bRnR0NB4eHjRo0MAutThaiWhRiMj9wD1AOeALY8yvTi6p0IL8ghh31zheavUSH677kEkbJjFr2yy61+/O621fp2FQQ2eXqFSxV5C//AHmxszl/tn3M6v7LO6rc1+h3/fkyZMsX76crVu3IiJkZGQgIvTu3ZtNmzbRpUuXXF+Xefop81RUSeHwFoWIfCkix0RkW471HUVkl4jsFZFh+e3DGPOTMeYJ4CngYUfWW9QCywYy5o4xxA2O45U2r/DL3l9o9Gkjun3b7YpNaqVUwdxX5z7WPb6OzrU722V/3333Hf/5z3+Ij48nLi6OgwcPEhERwR133MHXX3/N+vXrs7b94YcfOHrUOmvDAw88wMKFC5k9e3aJ6Z+Aojn1NA245FZxIuIOTAbuBuoDPUWkvog0FJEFOZbK2V76mu11LqdSmUq83f5t4ofE88atb7Asdhk3TrmRLrO6sPGwfWbCVaq0EhFahLRAROyyv6ioKLp27XrJum7dujFr1ixmzZrF0KFDqVOnDvXq1WPx4sX4+/sDUKFCBVq2bElQUBDVq1e3Sy1FId9pxm2/0LcbYwrV0yoi4cACY0wD2/OWwEhjTAfb8+EAxpgxebxegLHAEmPM0jy26Q/0BwgNDb0pPj6+MCU73ankU0xcP5H/++P/+Df5XzrV6sSItiNoHlL8h9Ip5Ug6zbh92G2acWNMBrBLRELtWyJVgYPZnifY1uXlWeAOoLuIPJXbBsaYz4wxkcaYyMDAQPtV6iQVfCrw+q2vEzckjnfav8P6hPW0+KIFHf/XkbUH1zq7PKVUKVKQU08Vge0iskxE5mUuji4sO2PMR8aYm4wxTxljPs1ru6KYPbaolfMux/A2w9k/eD/v3vEum//ZTKsvW3HH9DtYHb/6yjtQSqlCKkhQvA7cC7wFfJBtKYxDQLVsz0Ns6wqlKGePLWr+3v681Ool9g/ezwd3fcC2Y9toO60tt319G7/t/w1XvFOhUnnRz3vhXO3xu2JQGGNWAjGAv23ZaVtXGBuAWiISISJeQA+g0K0UV2xR5FTWqyzPt3ye/YP3M77DeHad2EX76e25ddqtLI1dqv+BlMvz8fEhMTFRP+vXyBhDYmIiPj4+BX7NFe+ZLSIPAeOAFYAAbYAXjTHfFegNRKKAdkAAcBR4wxjzhYh0AsYD7sCXxpjRBa76Cux5z+ziLjk9mS82f8HY38eScCaBliEtGXHrCDrU6GC3ER5KFSdpaWkkJCSQnJzs7FJKLB8fH0JCQvD09LxkfV6d2QUJir+BO40xx2zPA4GlxpjG9ivbPkSkM9C5Zs2aT+zZU7rmUUpJT+Gr6K8Ys2YMB04foFmVZoy4dQT31LpHA0MpVSDXNOopc5vMkLBJLODripwr91FcibeHN09FPsWeZ/cwtfNUjicdp3NUZyKnRjI3Zq4205VS16wgv/B/EZHFItJHRPoAPwMLHVvWtSkNfRRX4uXuRb+m/dg9cDdf3vclp5NPc//s+7lxyo18v+N7LMZy5Z0opVQ2V7rgTrCOSGoGtLatXm2M+bEIartmpamP4krSLelEbY3i7dVvsztxNw0qN+D1tq/TrV433N3cnV2eUqoYKUwfxVZjTImaoU6D4nIZlgy+3f4to1aNYueJndQLqMfrbV/noRse0sBQSgGF66PYLCLNHFCT3empp7y5u7nTs2FPtj69ldndZ+Pu5k6vH3pR/+P6fPP3N6Rb0p1dolKqmCpIiyIGqAnEA+exDpE1xphGji/v2miL4sosxsKPO3/krVVvseXoFmpUrMGrbV7l0UaP4unueeUdKKVczjWderL1UbTBGhKXMMYU21n3NCgKzmIszN81n7dWvcXmfzYTUSGCV9q8wmONH8PL3cvZ5SmlitC1TgpogMnGmPici8MqVUXKTdzoUrcLG5/YyPye8wkoE8AT85+g1sRafLrxU1LSU5xdolLKybSPQgHW+frvrX0v6/utZ9Eji6jiX4Wnf36amhNrMvnPySSn61WwSpVWBe2jqAXEoX0UpYYxhmX7l/HmyjdZc2ANwX7BvNzqZfrf1B9fT19nl6eUcoDCDI8Ny219cT79pEFhP8YYVsSt4K1Vb7EibgVBZYN4qdVLPHnTk5T1Kuvs8pRSdnTNw2NtgVANaG97nFSQ1ynXICLcFnEbv/X+jZV9VtKgcgNe+PUFIiZE8N7v73Eu9ZyzS1RKOdgVf+GLyBvAy8Bw2ypP4H+OLEoVT23D2rL0saWs6buGpsFNeXnpy4SPD2fM6jGcSTnj7PKUUg5SkJZBV+A+rP0TGGMOY70vRbGjndlFo1VoK3559Bf+ePwPmoc055XlrxA+PpxRK0dxKvmUs8tTStlZQYIi1TZM1gCISLE9MV2aZ491huYhzfm5189seGIDbcLaMGLFCMLHhzNyxUj+vfCvs8tTStlJQYLiWxGZAlQQkSeApcBUx5alSpLIKpHM7TGXzf030z6iPW+ufJOw8WG8tvw1EpMSnV2eUqqQrjjqCUBE7gTuwjo0drExZomjCysMHfXkXFuObuHtVW/z3Y7vKOtVloHNBvJ8y+cJLBvo7NKUUvm45uGxJZEGRfGw/dh2Rq8ezaxts/D19GVA5ACG3jKUIL8gZ5emlMpFYWaPVeqa3FD5BmZ2m8mOZ3bwQL0H+PCPD4mYEMHzi5/nn7P/OLs8pVQBaVAoh6sbUJdvun5DzDMxPNzgYT5a/xEREyIYtGgQCWcSnF2eUuoKChQUIuIrInUcXUxh6fDY4q1WpVp81eUrdg3cxaONHuWTjZ9Q46MaDPh5AAdOH3B2eUqpPBTkgrvOQDTwi+15ExGZ5+jCroUOjy0ZalxXg8/v+5w9z+6hT+M+fL75c2p+VJMn5z9J3Kk4Z5enlMqhIC2KkcDNwCkAY0w0EOHAmlQpEV4hnCmdp7B30F6eaPoE0/6eRq2JtXh87uPsO7nP2eUppWwKEhRpxpic53Jcb6iUcprQ8qFMvmcysYNiGRA5gJnbZlJnUh36/NSHPYl7nF2eUqVeQYJiu4j0AtxFpJaITATWOrguVQpVLVeVCXdPIHZQLIOaD+Lb7d9Sd3Jd/vPjf4g5EePs8pQqtQoSFM8CNwApwEzgNDDEkUWp0i3YP5gPO3zI/sH7eaHlC/yw8wfqT65Pz+97sv3YdmeXp1SpU5D7UTQ1xmwuonrsQi+4cy3Hzx/nw3UfMmnDJM6nnqd7/e681vY1GgUV23tnKVUiFeaCuw9EZKeIjBKRBg6oTal8BZYNZMwdY4gbHMerbV5l8b7FNP60MQ/MfoC//vnL2eUp5fIKcuOi24DbgOPAFBHZKiKvObwyGxGpJyKfish3IvJ0Ub2vKn4qlanEqPajiBscxxu3vsHy/ctp+llT7ou6j42HtQWplKMU6II7Y8wRY8xHwFNYr6kYUZDXiciXInJMRLblWN9RRHaJyF4RGXaF995pjHkKeAhoVZD3Va6tom9FRrYbSfyQeEbdNoo1B9bQbGoz7pl5D+sT1ju7PKVcTkEuuKsnIiNFZCuQOeIppID7nwZ0zLE/d2AycDdQH+gpIvVFpKGILMixVLa95j7gZ2BhQX8w5frK+5TntbavETckjnfav8P6hPW0+KIFHf7XgbUHdWCeUvZSkM7sdcBs4Fvb3e2u7g1EwoEFxpgGtuctgZHGmA6258MBjDFjCrCvn40x9+Txvf5Af4DQ0NCb4uPjr7ZUVcKdSz3HJxs+YdzacRxPOs7tEbcz4tYRtA1r6+zSlCoRrrkz2xjT0hgz/lpCIg9VgYPZnifY1uVKRNqJyEe2myfl2aIwxnxmjIk0xkQGBup9D0ojPy8/Xmz1IvsH7+eDuz5g27Ft3DrtVtpNa8dv+3/DFafUV6oo5BkUIvKt7etWEdmSbdkqIluKqkBjzApjzCBjzJPGmMn5bauTAiqAsl5leb7l8+wfvJ8JHSew5+Qe2k9vT9tpbVmyb4kGhlJXKb8WxWDb13uBztmWzOfX6hBQLdvzENu6QtNJAVV2vp6+DGo+iH2D9jHp7knEnYrjrv/dxS1f3sKiPYs0MJQqoDyDwhiTeWeZAcaY+OwLMKAQ77kBqCUiESLiBfQA7DIbrbYoVG58PHx45uZn2PvsXj6951MOnz1Mp5mdaP55cxbsXqCBodQVFGR47J25rLu7IDsXkShgHVBHRBJE5HFjTDowEFgM7MTaSW6XeRm0RaHy4+3hzZORT7Ln2T1M7TyVE0kn6BzVmZs+u4mfYn7SwFAqD3mOerJd3DYAqA5kn/PZH/jdGPOo48u7OrZ7Z3SuWbPmE3v26KyjKn9pGWnM2DqD0atHs/fkXhoHNeb1tq/TtV5X3ERv/qhKn7xGPeUXFOWBisAYIPtFcWeNMScdUqWd6FxP6mqkW9KZtW0Wb696m12Ju2hQuQGvt32dbvW64e7m7uzylCoyVz081hhz2hgTZ4zpaeuXuID1PhR+IhLqwFqVKlIebh482uhRtg/YzswHZpJhyeDh7x6m4ScNmbnV+lyp0qxAt0IVkT3AfmAlEAcscnBd10Q7s1VhuLu507NhT7YN2Ma33b/F3c2dR354hPof12f639NJt6Q7u0SlnKIgJ2LfBloAu40xEcDtwB8OreoaaWe2sgc3cePBGx7k76f+5vuHvsfXw5feP/Wm7qS6fPXXV6RlpDm7RKWKVEFvhZoIuImImzHmN+Cyc1hKuRo3ceOBeg/w15N/8dPDP1Hepzz/nfdfak+qzdRNU0nNSHV2iUoViYIExSkR8QNWATNEZAJw3rFlXRs99aQcQUToUrcLG5/YyIKeCwgsE0j/Bf2pNbEWn278lJT0FGeXqJRDFWRSwLJAMiDAI0B5YIatlVEs6agn5UjGGH7d9ytvrnyTdQnrqOpflWGth9GvaT98PHycXZ5S1+yqh8eWZBoUqigYY1i2fxlvrnyTNQfWEOwXzMutXqb/Tf3x9fR1dnlKXbWrHh4rImdF5Ey25Wz2r44tV6niT0S4o/odrOqzit96/0adgDoMWTyEiAkRfLD2A86nFssztEpdNZdqUeiV2crZVsWvYtSqUSyNXUpgmUCG3jKUAc0G4Ofl5+zSlLqiQp16EpHWQC1jzFciEgD4G2P2O6BOu9BTT8rZ1h5cy1sr32LxvsVU8q3E8y2fZ+DNAynnXc7ZpSmVp2u+cZGIvAG8DAy3rfIC/mff8pRyLbdUu4VfHv2FPx7/gxYhLXh1+auEjw9n1MpRnEo+5ezylLoqBRke2xW4D9uQWNud7vwdWZRSrqJ5SHMW9FrAxic20iasDSNWjCB8fDhv/PYGJy8U6ynTlMpSkKBINdbzUwayhssWS3odhSqubqpyE3N7zGVz/83cXv123lr1FuHjw3lt+WskJhXbkeZKAQULim9t96uuICJPAEuBqY4t69roFB6quLsx+Ea+f+h7tjy1hbtr3c07q98hfEI4w5YO4/j5484uT6lc5duZLSKC9ValdYG7sF50t9gYs6Royrs22pmtSortx7YzevVoZm2bha+nLwMiBzD0lqEE+QU5uzRVCl3zqCcR2WqMaeiwyhxAg0KVNDEnYnhn9TvM2DoDb3dvnrzpSV5s9SJV/Ks4uzRVilzzqCdgs4g0c0BNSimbugF1md51OjHPxPBwg4eZ+OdEqk+ozrMLnyXhTIKzy1OlXEFaFDFATSAe68gnAYwxppHjy7s22qJQJV3sv7G8s/odvv77a9zEjcdvfJxhrYcRWl7vGaYcpzCnnsJyW2+7612xpEGhXEXcqTjGrhnLl399CUCfJn14pc0rhFcId25hyiWVikkBdQoP5aoOnj7Iu7+/y9TNU7EYC481eoxX2rxCjetqOLs05UJKRVBk0haFclWHzhzivd/f47PNn5GWkUb3+t0JLBNIpTKVABh480ACygQ4uUpVUmlQKOVC/jn7D++vfZ9Jf04i1WK9054gLHxkIR1rdnRydaqkKsyoJ6VUMRPsH8wHHT7gwHMHqOBTAQCDYfKfk/kjoVje0l6VYBoUSpVgQX5BTOsyDYCeDXqyNmEtLb9oyW1f38aSfUtwxTMGquhpUChVwt1X5z7WPb6OGQ/MIH5IPB/e9SG7E3dz1//uotnUZvyw8wcsxuLsMlUJpkGhVAknIrQIaYGI4Oflx3MtnyN2UCxTO0/lVPIpun3bjQYfN+Dr6K9Jy0hzdrmqBNKgUMoFeXt4069pP2IGxhDVLQpPd0/6zO1DzYk1mfTnJC6kXXB2iaoEKRFBISJlRWSjiNzr7FqUKkk83Dzo0aAH0U9Gs6DnAqqVq8azi54lfEI4Y9eM5XSyTsmvrsyhQSEiX4rIMRHZlmN9RxHZJSJ7RWRYAXb1MvCtY6pUyvWJCPfUvoc1/13Dqj6raBrclOHLhhM6PpRXl73KsfPHnF2iKsYceh2FiLQFzgHTjTENbOvcgd3AnUACsAHoCbgDY3Ls4r9AY6AS4AOcMMYsuNL76nUUSl3Z5n82M3bNWL7b8R0+Hj70a9qPobcM1fmkSjGnXXAnIuHAgmxB0RIYaYzpYHs+HMAYkzMkMl8/GigL1AcuAF2NyX8IhwaFUgW368Qu3vv9PaZvmQ7Ao40e5eVWL1M3oK6TK1NFrThdcFcVOJjteYJtXa6MMa8aY4YAM4GpeYWEiPS39WNsPH5c7xSmVEHVCajDF12+YN+gfQyIHMDsbbOpP7k+3b/tzqbDm5xdnioGSkRnNoAxZlp+p52MMZ8BbwKbvby8iq4wpVxEaPlQJtw9gfgh8bzS5hWWxi4lcmokHf7XgZVxK/XivVLMGUFxCKiW7XmIbV2h6T2zlSq8wLKBvN3+bQ48d4Cxt48l+kg07b5uR6svW7Fg9wINjFLIGUGxAaglIhEi4gX0AObZY8ci0llEPjt9Wof8KVVY5bzL8XLrl4kbHMfkTpM5fPYwnaM60/jTxkRtjSLdku7sElURcfTw2ChgHVBHRBJE5HFjTDowEFgM7AS+NcZst8f7aYtCKfvz9fRlQLMB7Hl2D9Pvn066JZ1eP/Si7qS6fLbpM1LSU5xdonIwnWZcKXVVLMbCvF3zeGf1O2w4vIEq/lV4vsXzPBn5JH5efs4uTxVCcRr15DB66kkpx3MTN+6vez/r+61n6X+WUjegLkOXDCVsfBgjV4wkMSnR2SUqO9MWhVKq0NYnrGfMmjHM3TWXsp5leSryKZ5v+TxV/Ks4uzR1FbRFoZRymOYhzfmpx09sfXorXet1Zfwf44mYEMGT859k38l9zi5PFZK2KJRSdrf/3/2MWzuOL//6kjRLGg/f8DDDWg+jUVAjZ5em8lEqWhRKqeIhomIEH9/zMfsH72doy6HM3z2fxp82pnNUZ9YdXOfs8tRVcqmg0FNPShUvwf7BvHvnuxwYcoC32r3FuoPruOXLW2g3rR2/7vtVL94rIfTUk1KqyJxPPc/UzVN5f+37HDp7iJuCb2J46+F0rdcVN3Gpv1tLJD31pJRyurJeZRnSYgj7Bu3j886fczrlNN3ndKf+5PpMi56mt2otpjQolFJFztvDm8ebPk7MMzHM7j4bX09f+s7tS42PajBx/USS0pKcXaLKxqWCQvsolCpZ3N3ceeiGh9jcfzMLey0kvEI4g34ZRPj4cN5Z/Q6nkk85u0SF9lEopYqZ1fGrGbNmDIv2LqKcdzkGRA5gSIshBPkFObs0l6d9FEqpEqFNWBsWPrKQzf0307FmR979/V3CJ4Tz7MJniT8V7+zySiUNCqVUsXRj8I3M7j6bmIExPNLwEaZsmkLNiTXp81Mfdh7f6ezyShUNCqVUsVa7Um0+v+9zYgfHMrDZQObsmMMNH99At2+7sfGwnmIuCi4VFNqZrZTrCikXwv91/D/ih8TzWtvXWL5/Oc2mNuOub+5iRdwKvXjPgbQzWylVIp1JOcOUjVP4YN0HHD1/lBYhLXil9SvcU/sevXjvGmlntlLKpZTzLseLrV4kbkgcH3f6mCPnjnDfrPto/GljZm6dqbdqtSMNCqVUiebj4cPTzZ5mz7N7+KbrNxhjeOSHR6gzqQ5TNk4hOT3Z2SWWeBoUSimX4OHmwaONHmXL01v46eGfCCgTwFM/P0XEhAjeX/s+Z1POOrvEEkuDQinlUtzEjS51u/DH43+w7LFlNKjcgBeXvEjY+DDe+O0NvVXrNXCpoNBRT0qpTCJC+4j2LPnPEv7s9yftwtvx1qq3CBsfxvOLn+fQmUPOLrHE0FFPSqlSY/ux7bz7+7vM3DoTN3Gjd+PevNz6ZWpeV9PZpRULOupJKVXq3VD5BqZ3nc7eQXt5oukTfLPlG+pMqkPP73vy95G/nV1esaVBoZQqdcIrhDP5nsnEDYnjxVte5OfdP9NkShPunXkvvx/43dnlFTsaFEqpUut6v+sZe8dYDjx3gLdve5v1h9bT+qvW3DrtVhbvXaxXe9toUCilSr0KPhV4te2rxA2OY0LHCcT+G0vHGR256bOb+G7Hd2RYMpxdolNpUCillE1Zr7IMaj6IfYP28eV9X3I+7TwPznmQ+h/X56u/viI1I9XZJTqFBoVSSuXg5e5F3xv7smPADuY8OIeynmX577z/UvOjmny0/qNSd6tWDQqllMqDu5s73et3Z1P/TfzyyC9EVIxg8C+DCRsfxuhVo0vNrVqLfVCISDsRWS0in4pIO2fXo5QqfUSEDjU7sLLPStb0XUPzqs157bfXCP2/UIYtHcbRc0edXaJDOTQoRORLETkmIttyrO8oIrtEZK+IDLvCbgxwDvABEhxVq1JKFUSr0FYs6LWA6Cejuaf2PYxbO47wCeE88/MzxJ2Kc3Z5DuHQK7NFpC3WX/LTjTENbOvcgd3AnVh/8W8AegLuwJgcu/gvcMIYYxGRIOBDY8wjV3pfvTJbKVVU9iTu4b3f3+Prv7/GYiz0atiLYa2HUT+wvrNLu2p5XZnt8Ck8RCQcWJAtKFoCI40xHWzPhwMYY3KGRM79eAEzjTHd8/h+f6A/QGho6E3x8XoTdqVU0Uk4k8CH6z5kyqYpJKUlcX/d+3ml9Ss0q9rM2aUVWHGawqMqcDDb8wTbulyJyAMiMgX4BpiU13bGmM+MMZHGmMjAwEC7FauUUgURUi6EDzt8SPyQeEa0HcHKuJXc/PnN3PnNnSzfv7xEX7xX7DuzjTE/GGOeNMY8bIxZkd+2OnusUsrZAsoE8OZtbxI/JJ5xd45j27Ft3D79dlp+0ZK5MXOxGIuzS7xqzgiKQ0C1bM9DbOsKzRgz3xjTv3z58vbYnVJKXTN/b3+G3jKU/YP38+k9n3Ls/DHun30/jT5pxIwtM0rUrVqdERQbgFoiEmHrd+gBzLPHjrVFoZQqbnw8fHgy8kl2P7ubGQ/MQER49MdHqT2xNp9u/LRE3KrV0cNjo4B1QB0RSRCRx40x6cBAYDGwE/jWGLPdHu+nLQqlVHHl4eZBr4a9+Pupv5nXYx5BfkE8/fPTREyIYNzv44r1rVpd6sZFItIZ6FyzZs0n9uzZ4+xylFIqT8YYVsav5J3V77AkdgkVfCrw7M3PMqj5IALKBDilJqcNj3UGvY5CKVWSbDy8kTFrxvDDzh8o41mG/k3788ItLxBSLqRI6yhOw2MdRvsolFIlUWSVSL5/6Ht2DNjBg/UfZOKfE6k+oTr95vVjT6Lzz45oi0IppYqZ+FPxvL/2fT7/63NSM1LpXr87w1sPp8n1TRz6vqWiRaGUUq4grEIYEztNJG5wHC/d8hK/7P2FG6fcSKcZnVhzYE2R1+NSQaGnnpRSriTIL4gxd4whfkg8o9uPZuPhjbT5qg1tv2rLoj2Liuxqbz31pJRSJURSWhJfbP6CcWvHcfDMQZpc34ThrYfTrV433N3cC71/HfWklFIuIjUjlZlbZzJ2zVh2Je6ieoXq3FD5BhoHNcbdzZ2BNw+8piG2GhRKKeViMiwZ/BTzE8OWDmPvv3sBEISFjyykY82OV72/UtGZrX0USqnSxN3NnW71u7Fr4C6q+FUBILxCOB1qdLDr+7hUUOgUHkqp0sjNzY2P7/kYgPEdxyMidt2/h133ppRSyinuq3Mf6x5fR/Oqze2+bw0KpZRyASJCi5AWDtm3S516UkopZX8uFRTama2UUvbnUkGhndlKKWV/LhUUSiml7E+DQimlVL40KJRSSuXLJafwEJHjwHnghLNrKcYC0ONzJXqM8qfH58pK2jEKM8YE5lzpkkEBICIbc5uzRFnp8bkyPUb50+NzZa5yjPTUk1JKqXxpUCillMqXKwfFZ84uoJjT43Nleozyp8fnylziGLlsH4VSSin7cOUWhVJKKTvQoFBKKZUvlwgKEYkTka0iEi0iG23rrhORJSKyx/a1orPrLEoi8qWIHBORbdnW5XpMxOojEdkrIltEpKnzKi8aeRyfkSJyyPY5ihaRTtm+N9x2fHaJiH1vH1YMiUg1EflNRHaIyHYRGWxbr58hm3yOket9jowxJX4B4oCAHOveA4bZHg8D3nV2nUV8TNoCTYFtVzomQCdgESBAC2C9s+t30vEZCQzNZdv6wN+ANxAB7APcnf0zOPj4BANNbY/9gd2246CfoSsfI5f7HLlEiyIPXYCvbY+/Bu53Yi1FzhizCjiZY3Vex6QLMN1Y/QFUEJHgoqnUOfI4PnnpAswyxqQYY/YDe4GbHVZcMWCM+ccYs9n2+CywE6iKfoay5HOM8lJiP0euEhQG+FVENolIf9u6IGPMP7bHR4Ag55RWrOR1TKoCB7Ntl0D+H3hXNtB26uTLbKcrS/XxEZFw4EZgPfoZylWOYwQu9jlylaBobYxpCtwNPCMibbN/01jbfToOOBs9Jrn6BKgBNAH+AT5wbjnOJyJ+wPfAEGPMmezf08+QVS7HyOU+Ry4RFMaYQ7avlT3kPQAABDpJREFUx4AfsTbnjmY2fW1fjzmvwmIjr2NyCKiWbbsQ27pSxRhz1BiTYYyxAFO5eFqgVB4fEfHE+gtwhjHmB9tq/Qxlk9sxcsXPUYkPChEpKyL+mY+Bu4BtwDygt22z3sBc51RYrOR1TOYBj9lGrrQATmc7vVBq5Din3hXr5wisx6eHiHiLSARQC/izqOsrSiIiwBfATmPMh9m+pZ8hm7yOkUt+jpzdm17YBaiOdSTB38B24FXb+krAMmAPsBS4ztm1FvFxicLa7E3Dei708byOCdaRKpOxjsLYCkQ6u34nHZ9vbD//Fqz/qYOzbf+q7fjsAu52dv1FcHxaYz2ttAWIti2d9DNUoGPkcp8jncJDKaVUvkr8qSellFKOpUGhlFIqXxoUSiml8qVBoZRSKl8aFEoppfKlQaGUHYnIEBEpk+35QhGp4ID3CRaRBfl830tEVomIh73fW5U+GhRK2dcQICsojDGdjDGnHPA+z2O96jdXxphUrNc7POyA91aljAaFKrVE5DHbxG1/i8g3tnXTRORTEdkoIrtF5F7b+j4iMinbaxeISLsc+xsEVAF+E5HfbOviRCRARMJFJMa2/90iMkNE7hCR3233drjZtn1Z20Ryf4rIXyLSJY/yuwG/2F5zg237aNvPU8u2zU/AI3Y7YKrU0mapKpVE5AbgNeAWY8wJEbku27fDsc7PUwPrL/2aBdmnMeYjEXkeuO3/27t71qiCKIzj/6Oo6VNZGPGlsBCJiB/AfABBZC0sU6WxNkQ78QOY0koELVTUwlfQoIUICqImqxgxAVsDFhI28WUfi5lNxnVzkWtQNM8PFvbODnPPFruHuXPvGUlzPbrsBBrAMPAUOEp6uvcgMEYq2X0CmJA0nC9ZPYmIe5Lmi9i3AR8lLeamEeCMpAsRsRFYn9ungP2/ErtZFc8obK0aAi53/tAllXtTXJLUlvQWmAF2rdI5ZyVNKhWLawL3lUojTJKSE6RaZaMR8Rx4APQBA13jbAY+FMePgbGIOA5sldTK3+kb8LlTC82sLicKs59117UR8JUffy99NcZdLN63i+M2y7P7AA5LGsyvAUmvu8ZpleeXdJE0K2kBtyJiqOi7CVioEavZEicKW6smgEZE9EPaC7r4rBER6yJiB6no5BvSdruDuX0LK+9M9om0LWZdd4FjuTIpEbG3R59plmcgRMR2YEbSOKma657c3g/MSfryG/GYOVHY2iSpCZwGHkbEC6Aspf2eVP75NjAiaQF4BMwCr4Bx4NkKQ58F7nQWs2s4BWwAXkZEMx93xz4PvCvWTo4AU/ly1W7gfG4/ANysGYfZElePNStExDnghqQrfzuWKhFxCNgn6WRFn6vAqKTpPxeZ/Y9815PZP0jStc5ls17y3U/XnSRsNXhGYWZmlbxGYWZmlZwozMyskhOFmZlVcqIwM7NKThRmZlbpO8QbTwalPAxSAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## For fixed rotation $\\eta=.05$ and varying relative error ($x$-axis), compute CPU time ($y$-axis)"
      ],
      "metadata": {
        "id": "L2FDRl8G3c9R"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "id": "ALte7iV5x1Bh"
      },
      "outputs": [],
      "source": [
        "import time\n",
        "\n",
        "aa = .05\n",
        "max_cpu_time = 1500  # in sec\n",
        "max_iters = 3000\n",
        "T, K = 10, 1\n",
        "last_K = max_iters - (T-1) * K\n",
        "n = (T-1) * K + last_K\n",
        "la_k, la_alpha = 5, .5  # hyperparameters for LA-GDA \n",
        "gamma = .1  # Step-size for GDA, EG, OGDA, LA-GDA\n",
        "A1 = np.concatenate((aa*np.identity(dim),(1-aa)*np.identity(dim)),axis=1)\n",
        "A2 = np.concatenate((-(1-aa)*np.identity(dim),aa*np.identity(dim)),axis=1)\n",
        "A = np.concatenate((A1,A2),axis=0)\n",
        "t_acvi, t_gda, t_eg, t_ogda, t_lak = [], [], [], [], []\n",
        "eps = [.5, .1, .05, .01, .001, .0005]\n",
        "\n",
        "for de in eps:\n",
        "  # ACVI ------------\n",
        "  beta, mu, delta = .5, 1e-6, .5\n",
        "  y, lmd = init,np.zeros(((2*dim,1)))\n",
        "  t0 = time.time()\n",
        "  for t in range(T):\n",
        "    mu *= delta\n",
        "    flag = 1\n",
        "    for k in range(K if t+1 < T else last_K):\n",
        "      x = linalg.solve(np.identity(2*dim)+np.dot(Pc,A)/beta,np.dot(Pc,y-lmd/beta)+dc)\n",
        "      if np.linalg.norm(x-x_opt)/np.linalg.norm(x_opt) <= de or \\\n",
        "         time.time()-t0 >= max_cpu_time:\n",
        "        flag = 0\n",
        "        t_acvi.append(time.time()-t0)\n",
        "        break\n",
        "      res = opt.minimize(lambda z: f_y(z,x+lmd/beta,mu), y, method='SLSQP', jac=None)\n",
        "      if not res.success:\n",
        "        print('Failed to update y. t={}. k={}'.format(t,k))\n",
        "        flag = 0\n",
        "        break\n",
        "      y = res.x.reshape(-1,1)\n",
        "      lmd += beta*(x-y)\n",
        "    if flag == 0:\n",
        "      break\n",
        "\n",
        "  # GDA ------------\n",
        "  z = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    z = proj(z-gamma*np.dot(A,z))\n",
        "    if np.linalg.norm(z-x_opt) / np.linalg.norm(x_opt) <= de or \\\n",
        "        time.time()-t0 >= max_cpu_time:\n",
        "      t_gda.append(time.time()-t0)\n",
        "      break\n",
        "\n",
        "  # EG ------------\n",
        "  z = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    zz = proj(z-gamma*np.dot(A,z))\n",
        "    z = proj(z-gamma*np.dot(A,zz))\n",
        "    if np.linalg.norm(z-x_opt)/np.linalg.norm(x_opt) <= de or \\\n",
        "        time.time()-t0 >= max_cpu_time:\n",
        "      t_eg.append(time.time()-t0)\n",
        "      break\n",
        "  # OGDA ------------\n",
        "  z0 = np.copy(init)\n",
        "  z1 = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    z = proj(z1 - 2*gamma*np.dot(A,z1) + gamma*np.dot(A,z0))\n",
        "    z0 = np.copy(z1)\n",
        "    z1 = np.copy(z)\n",
        "    if np.linalg.norm(z-x_opt)/np.linalg.norm(x_opt) <= de or \\\n",
        "        time.time()-t0 >= max_cpu_time:\n",
        "      t_ogda.append(time.time()-t0)\n",
        "      break\n",
        "\n",
        "  # LAK-GDA ------------\n",
        "  la_k, la_alpha = 5, .5 \n",
        "  z = np.copy(init)\n",
        "  t0 = time.time()\n",
        "  for t in range(n):\n",
        "    z_old = z\n",
        "    for i in range(la_k):\n",
        "      z = proj(z-gamma*np.dot(A, z))\n",
        "    z = proj(la_alpha*z_old + (1-la_alpha)*z)\n",
        "    if np.linalg.norm(z-x_opt) / np.linalg.norm(x_opt) <= de or \\\n",
        "        time.time()-t0 >= max_cpu_time:\n",
        "      t_lak.append(time.time()-t0)\n",
        "      break"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "ms = 3\n",
        "print('GDA:      {}'.format(t_gda))\n",
        "print('EG:       {}'.format(t_eg))\n",
        "print('OGDA:     {}'.format(t_ogda))\n",
        "print('LA{}-GDA: {}'.format(la_k, t_lak))\n",
        "print('ACVI:     {}'.format(t_acvi))\n",
        "\n",
        "# plt.plot(eps, t_gda, 'ys-', label='GDA', markersize=ms) <- GDA too slow\n",
        "plt.plot(eps, t_eg, 'ch-', label='EG', markersize=ms)\n",
        "plt.plot(eps, t_ogda, 'mD-', label='OGDA', markersize=ms)\n",
        "plt.plot(eps, t_lak, 'bo-', label='LA{}-GDA'.format(la_k), markersize=ms)\n",
        "plt.plot(eps, t_acvi, 'gv-', label='ACVI', markersize=ms)\n",
        "plt.xscale('log')\n",
        "plt.xlabel('relative error')\n",
        "plt.ylabel('cpu time (s)')\n",
        "plt.legend(loc='upper right')\n",
        "plt.savefig('acvi_hd_bilinear_time_vs_relErr.pdf', bbox_inches='tight')\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "3QGVmddj5Lq8",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 372
        },
        "outputId": "208fd1b3-2315-4277-9b80-31637c4e1a24"
      },
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "GDA:      [78.64775919914246]\n",
            "EG:       [15.134716033935547, 161.77575206756592, 226.96552443504333, 376.7856979370117, 587.2638783454895, 656.7127356529236]\n",
            "OGDA:     [7.774897813796997, 84.55064845085144, 116.01626467704773, 194.55628848075867, 293.76169443130493, 329.068234205246]\n",
            "LA5-GDA: [14.562126159667969, 157.92764544487, 218.11382746696472, 357.96415185928345, 564.7706151008606, 627.1676275730133]\n",
            "ACVI:     [0.13395404815673828, 33.48216962814331, 47.89238142967224, 78.29657125473022, 125.12867450714111, 140.62643718719482]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3zURfrA8c+k916oCemAQAIkFA85TlE5QQFRKaci9oZiOQUscCp2zwLqKRbU3xFUUEGw0rzzqEE6pFASCCWkkYT07D6/P76bkEASSsqmzPv12ld2vzu7O+HOPDvzzDyjRARN0zRNA7Cxdgc0TdO0lkMHBU3TNK2KDgqapmlaFR0UNE3TtCo6KGiapmlVdFDQNE3TqthZuwMN4efnJ926dbN2NzRN01qVLVu2ZImIf23Pteqg0K1bNxISEqzdDU3TtFZFKZVW13N6+kjTNE2rooOCpmmaVkUHBU3TNK1Kq84paJrWtpWXl5Oenk5JSYm1u9IqOTk50aVLF+zt7c/7NTooaJrWYqWnp+Pu7k63bt1QSlm7O62KiJCdnU16ejohISHn/bp2OX2UXV7Oh0ePkl1ebu2uaJpWj5KSEnx9fXVAuAhKKXx9fS94lNXugkJSURHdNmxgakoKIRs2kFRUZO0uaZpWDx0QLt7F/Nu1u6CQXFREudlMmQglZjN/FBRYu0uaprVgtra2xMTEVN1efvllACoqKpg5cyYRERFVz82ZM8fKvW24dpdTiHRxwcHGBpPZTLkI0/btw9fenqt8fKzdNU3TWiBnZ2e2bdt21vWnn36a48ePs3PnTpycnCgoKOCNN96wQg8bV7sbKUS5uHBw0CDejYjg1z598LW35+odO7g/OZlCk8na3dM0rYGaI2dYVFTE/PnzmTt3Lk5OTgC4u7sze/bsJvvM5tLuRgoAvvb23N2pEwBb+vfn6YMHeTM9nV9ycvisRw/+5Olp5R5qmnamaSkpbDt1qt42RSYTf5w6hQD3JyfTz80NF1vbOtvHuLnxVkREve9ZXFxMTExM1eMZM2bQo0cPgoKCcHd3v6DfoTVodyOFMznb2vJGeDhrYmIwAUO3bmX6/v2Ums3W7pqmaReoyPLfrfmMxw1ROX1UeRs/fvxZbT799FNiYmLo2rUrhw8fbvBnWlO7HCnU5s9eXuyIjeXR/ft55fBhfsjJ4fPu3Ylpg98ENK01Otc3ejBWF8Zt2UKFCHZKsaRXL6JcXBq9L+Hh4Rw6dIiCggLc3d2ZMmUKU6ZMoVevXpha+TR0ux8pVOduZ8f8qCiW9+5NZnk5A/74gzlpaVToUYOmtQqVOcO3wsM5OGhQkwQEABcXF+644w4efPDBqn0AJpOJsrKyJvm85tQug8L69fDSS8bP2oz09WVXXBzX+/nx9MGDDNm6lWS9n0HTWoXKnKHvBZR2qE9lTqHyNn36dADmzJlDx44d6dWrF3379uWyyy5j8uTJdLLkK1srJSLW7sNFi42NlQs9T2H9erj8cigtBScnWLUKBg+uu/2ijAzuT0mhxGzmldBQHujcGRu9mUbTmsXevXvp0aOHtbvRqtX2b6iU2iIisbW1b3cjhbVrjYAgAsXF8O9/199+QmAgu+Li+LOXFw/t28eV27dzSBfn0jStjWp3QWHYMGOEYGP5zd97Dx56CPLy6n5NJ0dHfujdmw8iI9mYn0/vzZv57PhxWvMoS9M0rTbtLigMHmxMGb3wAvz8M9x/P8ybBz16wJdfGiOI2iiluLtTJ3bExRHt5sZtiYmM3bWLE20gsaRpmlap3QUFMALDjBlw1VVGQNi4ETp1ggkTYMQI2Lev7teGOjuzJiaG10JD+TEnh0s2b+abzMzm67ymaVoTapdB4UxxcUZgmDsXNmyAXr3gH/+AulIHtkrxeFAQf8TGEuToyLjdu7l1715O6lLcmqa1ck0aFJRSXkqpxUqpRKXUXqXUYKWUj1LqV6VUiuWnt6WtUkq9o5Tap5TaoZTq15R9O5OtLTz4ICQmwtixMHs29OkDK1fW/ZpLXF3Z0K8fzwYHszAjg94JCfyak9NsfdY0TWtsTT1SeBv4SUS6A9HAXmA6sEpEIoBVlscAfwUiLLe7gfebuG+16tgR4uPhl1+M/MKVV8KkSXD8eO3t7W1s+EdICOv79cPN1parduzgAV1cT9PajPT0dEaPHk1ERARhYWE8/PDDVZvUNm3axLBhw4iIiKBfv36MHDmSnTt3AjB79mw6d+5MTEwMERERXH/99ezZs6fGe2/btg2lFD/99FOz/151abKgoJTyBIYCHwOISJmInARGA59Zmn0GjLHcHw18LoYNgJdSqmNT9e9crrwSdu6EWbNgyRKIioJ334W6/tbHeXjwR//+TOvShfeOHiUmIYF19S1p0jStxRMRrr/+esaMGUNKSgrJycmcOnWKp556ioyMDG666SZefPFFUlJS+OOPP5gxYwb79++vev0jjzzCtm3bSElJYfz48Vx++eVkVstBxsfHM2TIEOLj463x69WqKUcKIUAm8KlSaqtS6iOllCsQKCLHLG2OA4GW+52B6pWk0i3XrMbJyZhG2rULBg40ppcGDYK69ss529ryZng4a6KjKTebuWzrVmYcOKCL62laM8pYmMFvjr+RsTCjwe+1evVqnJycmDJlCmAcuPPmm2/yySef8PrrrzN58mQuvfTSqvZDhgxhzJgxtb7X+PHjueqqq1i4cCFgBJyvv/6aBQsW8Ouvv17wsZlNpSkL4tkB/YCpIrJRKfU2p6eKABARUUpd0GJ/pdTdGNNLBAUFNVZf6xURYSxf/fJLeOQRGDAAHnjAWNZaW5XtYd7e7IiL45F9+3j50CF+yM7m8x49iHZza5b+alpblDIthVPb6i+dXZZRRnFyMZhh7y17SX0+FYdAhzrbu8W4EfFW3YX2du/eTf/+/Wtc8/DwICgoiP379zN58uQL+h369etHYmIiAOvWrSMkJISwsDCGDRvGihUrGDdu3AW9X1NoypFCOpAuIhstjxdjBImMymkhy88TluePAF2rvb6L5VoNIvKhiMSKSKy/v3+Tdf5MShlLVhMTjYDw7rvQvTssWlT73gYPOzs+7t6dZb16kVFWRtyWLbyUlkZGaWmTHwCiae1VcVLx6brZZsvjZjJw4EB69OjBww8/XGeb6hte4+PjmTBhAgATJkxoMVNITTZSEJHjSqnDSqkoEUkCrgD2WG6TgZctP5daXrIMeFAptQgYCORVm2ZqMTw9jaWrkyfDvffCxInw8cfGzujaKvte6+fHLg8P7ktJYebBgzx98CAOSvH4/v1s7t+/yao4alpbU983+koZCzNIujMJc7EZG2cboj6KInBS4DlfV5eePXuyePHiGtfy8/M5dOgQw4cP548//mD06NEAbNy4kcWLF7N8+fI632/r1q3ExsZiMplYsmQJS5cuZc6cOYgI2dnZVaW4rampVx9NBf6tlNoBxAAvYgSDK5VSKcBwy2OAH4ADwD5gPnB/E/etQWJjjb0N8+bBpk3Qu7eRf6htWtDPwYGvevbksS5dMAMlIpSZzewuLGzubmtamxY4KZCoj6JQDqrBAQHgiiuuoKioiM8//xwwymM/9thj3HbbbTz22GMsWLCAdevWVbUvqqea8pIlS/jll1+YOHEiq1atok+fPhw+fJjU1FTS0tIYN24c3377bYP62yhEpNXe+vfvLy3B0aMiEyeKgEh4uMjPP9feLrGwUNx++01s16wR1qyRkHXrZE1OTvN2VtNakT179li7C3Lo0CEZNWqUhIeHS2hoqDz44INSUlIiIiLr16+XoUOHSlhYmAwePFiuvfZa2bx5s4iIzJo1Szp16iTR0dESHh4uY8aMkd27d4uIyG233Sbvv/9+jc9ZunSpjBgxotH7X9u/IZAgdfxdbXels5vSypVGLaWUFCP/8M9/GvseqssuL2dJZiZetrZMP3iQgyUl3NahA6+FhuLnUHdCTNPaI106u+F06WwrGj4cduwwSmR8+62RiJ43r+behsoDQG6ylOSeHhTE/2Vk0H3TJl15VdM0q9NBoZE5OcGzzxp7GwYNgqlTjT0OtQ1oXGxteSk0lD/69yfSxYXbEhO5fPt2kvQpb5qmWYkOCk0kPBx++slYsnr0qLG34cEHjfIZZx4F2tvNjd/79uVfkZFsLSigz+bN/CM1VW960zSt2emg0ISUgvHjYe9eIyC89x5cfTU8/TRccUXNwGCjFPd06kTigAGM8/dndmoqfTZvZm1urvV+AU3T2h0dFJqBpye8846RhAYwm42jQL/++uy2HRwdWdizJz/16UO5CH/Zvp3b9u4lSx/mo2laM9BBoRn97W/g7Hz6KNC5c438Q3Etmy6v9vFhV1wcM4KC+PeJE3TftIkFx47pRLSmaU1KB4VmVP0o0O+/N6aWnn/eONSntsq5Lra2vBgaylbLzucpSUn8Zds2EvWmN01rNm711CybNm0anTt3xlwt/7d27Vo8PT2JiYkhJiaG5557rs7XZ2RkMGnSJEJDQ+nfvz+DBw+u2sBW+T59+/YlKiqKoUOHnrVbuqKiAn9/f6ZPn17b21+cujYwtIZbS9m81hCrVolERRkb3264QSQ9vfZ2JrNZPjhyRLz++19xWLtWZh04IMUVFc3bWU1rZi1h85qrq2ut100mkwQFBcnAgQNl9erVVdfXrFkjI0eOPOf7ms1mGTRoUI1NbKmpqfLOO+/U+j5bt26V4OBgWblyZdW1H374QS699FIJDQ0Vs9lc6+dc6OY1PVKwsssvh+3bjdHD8uXG3oY334SKiprtbJTi7k6d2BsXxzh/f/6RlkZ0QgJrdCJa02pYv/7sFX5NYe3atVxyySXcd999F1XMbvXq1Tg4OHDvvfdWXQsODmbq1Km1to+JieHZZ59l3rx5Vdfi4+N5+OGHCQoKYn0j/cJNWTpbO0+OjvDUU0ZxvalT4dFH4bPP4F//MvY6VFeZiL6tQwfuS07m8u3bmRwYyOthYXpHtNamTZsG27bV3yYvz9hAajYbubs+fWovb18pJgbeeuvi+hMfH8/EiRMZPXo0M2fOpLy8HHt7ewDWr19PdHQ0nTp14vXXX+eSSy456/W7d++mX78LO3W4X79+vPbaawCUlJSwcuVKPvjgA06ePEl8fHyNsx0ulh4ptCChocZoYckSyM6GSy+Fe+6B2o59vsqSiJ5pSURHbdrEpzoRrbVzeXlGQADjZ1MdflhWVsYPP/zAmDFj8PDwYODAgfz888+A8Yc7LS2N7du3M3Xq1DoP3TnTAw88QHR0NHFxcXW2qf7f9/Lly/nLX/6Cs7Mz48aN47vvvsPUGMcA1zWv1BpubSGnUJf8fJFHHxWxtRXx9xdZsECkjilD2XXqlPxpyxZhzRr58x9/yN5Tp5q3s5rWRC40p7BunYizs/HfjbOz8bihasspLFu2TFxcXCQ4OFiCg4PF399fJk2aVOvrg4ODJTMzU+bNmyfR0dESHR0tR44ckZUrV8rQoUNrtM3MzJTg4GARqT038fHHH8vYsWNFROT666+XgICAqj44OzvLL7/8ctbnX2hOwep/2Btya8tBodL27SKDBxv/Sw0dKrJrV+3tTGazfGhJRNuvXSvP6kS01gZcTKJ53TqRF19snIAgUntQmDhxoixcuLDq8alTp8Tf318KCwvl2LFjVUnfjRs3SteuXWtNApvNZhkwYIC89957VdfS0tLqDArbt2+Xbt26ycqVKyUvL0/8/f2rqrWKiHzyyScyZcqUsz5HB4U2yGQSmT9fxMdHxM5O5MknReoaDBwvLZVJu3cLa9ZIxIYNskqX5tZasZaw+kgpJZ07d666zZkzR7y9vSUvL69Gu7Fjx8qiRYtk7ty50rNnT+nTp48MHDhQ/ve//9X53kePHpXx48dLt27dJC4uToYNGyaLFi0SESMoeHh4SExMjERGRsqQIUNk2bJlIiKyYMECGT9+fI33ys7OFj8/vxqBQkSXzm7TMjPhySfh008hONjYJX3ddbW3/SUnh/uSkzlQUsKtlkS0v05Ea62MLp3dcLp0dhvm7w+ffAL/+Q+4ucHo0TBmDBw6dHbb6onohZYd0ToRrWnaueig0Apddhls3Qqvvgq//go9ehj3y8trtnO2tWVOaCjbYmPp4eLC7UlJDNu2jb16R7SmaXXQQaGVsreHv//dqMB61VXGtFLfvvDf/57d9hJXV/7Tty/zIyPZUVhIdEICzx48SEljLF/TNK1N0UGhlQsKMk55W7YMTp2CoUNhyhQj/1CdjVLcaSnNfZO/P8+npdEnIYFVeke0pmnV6KDQRlx7LezeDdOnw//9n1EuY/780xt5KgU6OPB/PXvyS58+mEUYvn07t+7dS6Yuza1pGk0cFJRSqUqpnUqpbUqpBMs1H6XUr0qpFMtPb8t1pZR6Rym1Tym1Qyl1Yfu/NVxdjZov27cblVfvvhuGDDEen+lKHx92xsXxVFAQiyyJ6E90IlrT2r3mGCn8RURiqi1/mg6sEpEIYJXlMcBfgQjL7W7g/WboW5vUsyesXWvUT0pJgf79jXpKBQU12znb2vJCaChbY2Pp6erKHToRrWm1+u6771BKkZiYWHVt06ZNDB06lKioKPr27cudd97Jnj176NKlS41S2mAUs9u4cSOzZ8/m9ddfb+7uXxBrTB+NBj6z3P8MGFPt+ueWvRUbAC+lVEcr9K9NUApuvRWSkuDOO42iXz16GHWVzhwMXOLqym8xMXwUFcVOnYjWtLPEx8czZMiQqmqoGRkZ3HjjjbzyyiskJSWxdetWRowYga+vL0FBQfy32oqPxMRECgoKGDhwoLW6f0GaOigI8ItSaotS6m7LtUAROWa5fxwItNzvDByu9tp0yzWtAXx8jGqr69YZ+xxuuAFGjoQDB2q2s1GKOzp2JHHAAMYHBPB8Whq9ExJYWVs1Pk1robKKspi9dnbVLasoq8HveerUKX7//Xc+/vhjFi1aBMC7777L5MmTGTx4cFW7G264gcDAQCZOnFjVDmDRokVMmDChwf1oLk1dOnuIiBxRSgUAvyqlEqs/WbmF/ELe0BJc7gYICgpqvJ62cYMGwebNMG8ePPMMXHKJUa777383SndXCnBw4IsePZgcGMh9KSlcuWMHNwcG8kZYGAF6R7RmRdN+msa24/XXzs4pzmHniZ1Vj7/Z+w0+zj51to/pEMNbI+qvnb106VJGjBhBZGQkvr6+bNmyhV27djF58uRa2990003ExMQwd+5c7Ozs+PLLL/m6tgPZW6gmHSmIyBHLzxPAt8AAIKNyWsjy84Sl+RGga7WXd7FcO/M9PxSRWBGJ9ff3b8rutzl2dkZN+sREY7XSM89AdDSsXn122+E+PuyIjeXp4GC+tCSiPz52DLNORGstmLeTN052TgA42znj7eTd4PeMj4+v+qY/YcKEcx6oExgYSK9evVi1ahXbtm3Dzs6OXr16NbgfzaXJRgpKKVfARkQKLPevAp4DlgGTgZctP5daXrIMeFAptQgYCORVm2bSGlHnzvDVV8a50A8+CFdcAX/7G7zxBgQGnm7nbGvL8yEhTAoI4J7kZO5MSuKz48f5V2QkPV1drfcLaO3Sub7RV1qauJQxX45h0Q2LuC6qjuJg5yknJ4fVq1ezc+dOlFKYTCaUUkyePJktW7YwevToWl9XOYVUOZ3UmjTlSCEQ+F0ptR3YBKwQkZ8wgsGVSqkUYLjlMcAPwAFgHzAfuL8J+6YBI0bAzp3w7LPw9dcQFQXvvQdn5pd7uLqy1pKI3lVYSExCAs8cPEixTkRrLdB1Udex/o71XBt5bYPfa/Hixdxyyy2kpaWRmprK4cOHCQkJYfjw4Xz22Wds3Lixqu0333xDRkYGANdffz0//PADX375ZavKJ0ATBgUROSAi0ZbbJSIyx3I9W0SuEJEIERkuIjmW6yIiD4hImIj0FpH2U/7Uipyd4R//MIJDbCw88ICRf9iypWa7MxPRL1h2ROtEtNbSKKUY1GUQSqkGv1d8fDxjx46tcW3cuHEsWrSIRYsW8fjjjxMVFUWPHj34+eefcXd3B8DLy4vBgwcTGBhIaGhog/vRnHTpbK2KCCxaZOxpOHEC7r8fXnih9jNuV+bkcF9KCvuKi/lbQAD/DA/XiWit0enS2Q2nS2drF00pmDjRSETff78xldS9O8THn723YbiPDztjY3kmOJivMjPpvmkTHx09qhPRmtbK6aCgncXTE+bOhU2boEsXmDTJqMSanFyznZOtLc+FhLA9NpZerq7clZzMn7dtY4/eEa1prZYOClqd+veHDRvg3XeNPQ69e8OsWVBSUrNdZSL646go9lgS0U8fOKAT0ZrWCumgoNXL1taYSkpMhBtvhOeeM4rt/fxzzXY2SnG7JRE9MSCAOYcO0XvzZn7ViWitgVpz3tPaLubfTgcF7bx06GCU5F650ggUI0bATTfBkTO2F/o7OPBZjx6sio7GRimu2rGDm/fs4YQuza1dBCcnJ7Kzs3VguAgiQnZ2Nk5OThf0Or36SLtgpaXw2mswZ45xAtxzzxmb4OzO2ApZYjLx0qFDvHToEG62trwaGsrtHTti0whLBbX2oby8nPT0dErOnLPUzouTkxNdunTB3t6+xvX6Vh/poKBdtAMHjGDw448QEwPvv2/scTjT3sJC7k1O5j95efzJw4MPoqK4RO+I1jSr0UtStSYRGgorVsDixcbxn5deCvfcA2emESoT0Z9ERbG3qIiYhASe0oloTWuRdFDQGkQpGDcO9u6FRx6Bjz829jZ8/nnNvQ1KKaZYEtGTAgJ48dAhem3ezC86Ea1pLYoOClqjcHc3Cupt2QLh4TB5MvzlL7BnT8121RPRtkpx9Y4d/G3PHjJ0IlrTWgQdFLRGFR0Nv/8OH34IO3YYj2fMgKKimu0u9/ZmR2wss4KDWWzZET1f74jWNKvTQUFrdDY2cNddxlGgN98ML79snBu9fHnNdk62tsy27IiOdnXl7uRkhm7dym69I1rTrEYHBa3J+PvDp5/Cf/4Dbm7GwT5jx8KhQzXbdXd1ZU1MDJ9WS0TP1IloTbMKHRS0JnfZZbB1K7zyCvzyC/ToAa++CuXlp9sopbjNkoj+W0AAL1kS0T/rRLSmNSsdFLRmYW8PTzxhJJ6vvBKefBL69TPyD9X5OziwoEcPVkdHY6cUI3bsYNKePRwvLbVOxzWtndFBQWtWwcHw3XewbBkUFBijiNtvh6ysmu3+4u3NdksieklmJj02b+ZDnYjWtCang4JmFddeC7t3GyOGL74wjgL96CMwm0+3qUxE74iLI8bNjXuSk7ls61Z2nTplvY5rWht33kFBKeWqlLJtys5o7Yurq7Eyads2o/LqXXcZI4cdO2q2i3JxYXV0NAu6dyepqIi+W7Yw88ABinQiWtMaXZ1BQSllo5SapJRaoZQ6ASQCx5RSe5RSrymlwpuvm1pbdsklsHYtLFhgHOTTrx88/jhUHxAopZjcoQOJAwZwc2AgL1lKc+tEtKY1rvpGCmuAMGAG0EFEuopIADAE2AC8opS6uRn6qLUDShm7oJOS4I47jN3RPXrAN9/ULJfh5+DAp92710hET9SJaE1rNPUFheEi8ryI7BCRqpleEckRkSUiMg748lwfoJSyVUptVUottzwOUUptVErtU0p9qZRysFx3tDzeZ3m+W8N+Na018vGBDz6A9evB19eoqzRqlFGRtbq/eHuzIy6O2d268Y1lR/QHOhGtaQ1WZ1AQkXIApVSYUsrRcn+YUuohpZRX9Tbn8DCwt9rjV4A3RSQcyAXusFy/A8i1XH/T0k5rpwYNgoQEePNNY/PbJZcY5zdUHxA42tgwq1s3dsTF0dfdnXuTkxmydSs7dSJa0y7a+SSalwAmSw7hQ6ArsPB83lwp1QUYCXxkeayAy4HFliafAWMs90dbHmN5/gpLe62dsrODadOMo0BHjYKnnzZqKa1ZU7NdZSL6s+7dSS4qot+WLczQiWhNuyjnExTMIlIBjAXmisjfgY7n+f5vAU8AldNPvsBJy/sBpAOdLfc7A4cBLM/nWdpr7VznzvD118ZhPuXlcPnlRk2ljIzTbZRS3GpJRN8SGMjLlh3RP2VnW6/jmtYKnU9QKFdKTQQmA5UlzezraQ+AUmoUcEJEtjSgf7W9791KqQSlVEJmZmZjvrXWwo0YAbt2wTPPGEEiKso47a36gMDPwYFPundnTXQ0Dkrx1507dSJa0y7A+QSFKcBgYI6IHFRKhQBfnMfr/gRcp5RKBRZhTBu9DXgppSpP8+0CVB79fgRjagrL857AWV/zRORDEYkVkVh/f//z6IbWljg7G2dC79gBsbFw//0weDD88UfNdsO8vdkeF8c/qiWi/3XkiE5Ea9o5nDMoiMgeEXlIROItjw+KyDmTwCIyQ0S6iEg3YAKwWkT+hrHU9QZLs8nAUsv9ZZbHWJ5fLa35AGmtSUVFwa+/wsKFcPgwxMXBQw9BXt7pNo42NjxrSUT3c3fnvpQUnYjWtHOob/Pa90qpa5VSZ00VKaVClVLPKaVuv4jPfBJ4VCm1DyNn8LHl+seAr+X6o8D0i3hvrR1RCiZONBLR998P8+YZexu+/LLm3oYoFxdWnZGInr5/v05Ea1otVF1fxpVSHTD+OI8DcoBMwAnoBuwH5onI0lpf3ExiY2MlISHBml3QWpCEBLj3XuNI0CuvhHffhYiImm2yysp44sABPj1+nG5OTrwfEcEIX72eQWtflFJbRCS21ufOZ4bGspGsI1AMJItIUb0vaCY6KGhnMpngX/+CmTOhpASmTzeOA3Vyqtnut5MnuTc5mcSiIsb4+nKppye3d+yIr/0511BoWqvX4KDQUumgoNXl+HF47DEj5xAebowarrqqZptSs5kn9u/nnSPGWgcnpfgjNpYerq5W6LGmNZ/6goIuna21SR06wL//DStXGmdGX301jB8PR4+ebuNoY8Nwb2/cbIz/DEpEuGH3bhL1GdFaO6aDgtamXXGFsXz1+eeNg326d4e334YKy/bJSBcXlFI429jgZGNDemkp0QkJPJeaSln1wx00rZ04r6CglHJWSkU1dWc0rSk4OholMnbvhiFDjNIZAwbAxo3Gyje6cfcAACAASURBVKSDgwbxVng46YMHkzxwIGP9/JiVmkrfhATWVV/jqmntwDmDglLqWmAb8JPlcYxSallTd0zTGltoKKxYAYsXw4kTxqa3++6DhNX2ZH/WieQEewIdHFh0ySUs792bApOJIVu38kByMvmVQwtNa+POmWhWSm3B2I28VkT6Wq7tFJHezdC/eulEs3axCgpg1ixjKslsNvIOjo6wapURLABOVVTw9MGDvHPkCJ0cHJgXEcEYvYteawMammguF5Ezx9Ctd8mSpgHu7vDPf8IDDxiPzWYoLjZWK1Vys7PjrYgINvTrh4+9PWN372bcrl0c1XWUtDbsfILCbqXUJMBWKRWhlJoLrGvifmlas5g40ainZFmAxLvvwoMPQm7u6TYDPDzY0r8/L4aEsCI7m576QB+tDTufoDAVuAQoBeKBfGBaU3ZK05rL4MHGlNELL8BPPxkB4f33ITISPvrIGEEA2NvYMCM4mJ2WOkr3Jifz523b9PJVrc3Rm9c07Qw7dhjB4b//NQrtzZtnrFaqJCIsOH6cx/bvp9BkYmZwMNODgnC00Su8tdahQTkFpVSsUuobpdQfSqkdlbfG76amtQx9+sBvvxmb39LTYeBAuPNOqDy+QynFlI4d2TtgANf7+zPbsnz1f3r5qtYGnM9Xm38DCzAK411b7aZpbZZSMGkSJCXB3/8On31mTCnNm3d641uggwPxPXuyondvCi3LV+9PTiZPL1/VWrHzCQqZIrLMco5CWuWtyXumaS2Auzu8+irs3GlMJU2dCv37G1NLla7x9WV3XBzTunThg6NH6blpE9/qUwG1Vup8gsIspdRHSqmJSqnrK29N3rMmlLEwg98cfyNjYca5G2saRnmMn3+Gb76Bkydh6FDjnOjKWkpudna8GR7Ohn798LO35/rdu7leL1/VWqHzPY4zBhjB6amjUU3ZqaaUsTCDxCmJSJmQOCWRY18cs3aXtFZCKRg7FvbuNc6JXrzYOAHu9dehrMxoE+fhQUL//rwUEsKPOTn00MeAaq3M+exoThKRFln36GJWH/3m+BtSVvN37nhPRwInBuJ5mSfKRjVmF7U2bP9+eOQR+P57YyTxzjvG4T6V9hUVcU9yMqtPnmSIpycfRkbqstxai9DQHc3rlFI9G7lPVtP90+7YOBu/tnJQuA92J+OLDLYN28b6oPXse2wf+Qn5tOalulrzCAszKq8uXw7l5cZ5DTfcAGmWjFu4iwsro6P5NCqK3YWFxCQk8I/UVEp19VWtBTufkcJeIAw4iLGBTQEiIn2avnv1u9h9CpVTSN0/7U7gpEBMhSayvs/iRPwJcn7MQcoF5whnAiYEEDAxANce+tudVr+SEqNsxgsvGI9nzoTHHz994tuJsjKm7dtH/IkT9HBxYX5UFH/y9LReh7V2rUEnrymlgmu73hJWIDXF5rXy3HKyvskiIz6Dk6tPgoBbjBsBEwMImBCAU5DTud9Ea7cOHTKCwddfG1VZ33oLrq22gPvH7GzuTU7mUGkp93bqxMuhoXja2Vmvw1q7dFFBQSnlISL5Simf2p4XkZxG7ONFaeodzaXHSsn8KpOM+AwKNhYA4PEnDwInBuJ/oz8OAQ5N9tla67ZqlbF8de9eGDnSCA7h4cZzpyoqeDY1lbfT0+lgqb46Vldf1ZrRxQaF5SIySil1EKMqavUMrIhI6Dk+1An4D+AI2AGLRWSWUioEWAT4AluAW0SkTCnlCHwO9AeygfEiklrfZzRnmYviA8WcWHSCE/EnKNxVCLbgfYU3ARMD8B/rj52n/ran1VReDnPnwuzZUFpqbIKbMQMqc82b8/O5KymJ7YWFjPXzY25EBJ0dHa3aZ619aND0UQM+VAGuInJKKWUP/A48DDwKfCMii5RS/wK2i8j7Sqn7gT4icq9SagIwVkTG1/cZ1qp9dGrXKU7EGwGi5GAJylHhe40vARMD8B3li62zbbP3SWu5jh2DJ5+EL76Arl2N3MO4ccYS13KzmX+mpzM7NRUHpXg5NJR7OnXCRulVcFrTaWjto1Xnc+1MYjhleWhvuQnGgT2LLdc/A8ZY7o+2PMby/BWWwNLiuPVyI3ROKAP3D6Tfhn50urcT+evz2XPTHtYFrmPvrXvJ/jEbc7leZaJBx47w+efGLmgfH7jxRmPp6p49RvXVJ4OC2BUXR5y7O/enpDB061b26OqrmpXUGRSUUk6WfIKfUspbKeVjuXUDOp/PmyulbJVS24ATwK/AfuCkiFQWh0mv9l6dgcMAlufzMKaYWiylFB4DPYh4K4LB6YOJXhmN/03+ZH+fzc5rdrKu4zqS70vm5H9OIma9xLW9GzIEtmwxzmzYsgWio+GxxyA/H8Kcnfk1OpoF3buzt6iImIQEZh88qJevas2uvpzCwxjnJnQCjnA6p5APzBeReef9IUp5Ad8CzwALRCTccr0r8KOI9FJK7QJGiEi65bn9wEARyTrjve4G7gYICgrqn5Zm9UVQZzGXmsn5OYcT8SfIWpaFuciMYxdH/Mf7Y+dtR9pzaVXLYbX2KTMTnnrKOLMhMNCor3TzzcaU0omyMh7Zt4+FJ07Q3cWF+ZGRDPHysnaXtTakoUtSp4rI3EboxLNAMfAk0EFEKpRSg4HZInK1Uupny/31Sik74DjgL/V0sDWcp1BxqoLs77M5sfAE2T9kQ+UXP1vo8mgXQp4LwdZJ5yDaq82bjbMbNm2CP/3JqMIaE2M892N2NvclJ5NWWso9HTvySliYXr6qNYoG5RQuNiAopfwtIwSUUs7AlcBeYA1wg6XZZGCp5f4yy2Msz6+uLyC0FnZudgRODKT3971RdtVSJCZIfy2d//n8j53X7eToB0cpSS+xXkc1q4iLg/Xr4eOPjTLd/fsbQSInB/7q68vuAQN4tEsX5h87Ro9Nm/hGV1/VmlhTrj7qg5E4tsUIPl+JyHNKqVCMJak+wFbgZhEptSxh/QLoC+QAE0TkQH2f0RpGCtVlLMwg6c4kzMVmbJxt6DS1E+ZCMzkrcihJNQKCa7QrviN98R3li8cAD5Rti8y1a00gNxdmzTJyDj4+8OKLcPvtYGsLCfn53JWczLZTpxjj58c8vXxVawCrLEltDq0tKMDZJTbAON6xaG8R2cuzyV6RTd7/8sAE9n72+IzwwXeUL95Xe2PvZW/l3mvNofpxoLGxxpTSwIHG8tU309OZpZevag3U0JzC0Nqui8h/GqFvDdIag8L5KM8tJ+fnHHJW5JD9YzYV2RVgC55/8sR3lC++I31x6eFCC12xqzUCEYiPN0pmHDtmjBheegkCAmB/cTH3JiezMjeXSz08mB8VRU9dfVW7AA0NCt9Xe+gEDAC2iMjljdfFi9NWg0J1YhLyN+aTvcIYRRRuN9avO4U44TvSF5+RPngN89LJ6jaqoACefx7efNPYCf3883DffWBrK3yRkcEj+/ZRYDIxIyiImcHBONqcT+Fjrb1r1OkjyzLSt0RkXGN0riHaQ1A4U8nhEnJ+yCF7RTa5K3ON/ISLDd7DvY1cxEhfHDvruea2JjHRqKW0ciX07m1MKQ0dCpmW5av/tixf/TAyksv08lXtHBo7KChgt4hY/YyF9hgUqjMVmzi59mRVLqI0zTj60S3GDd9RxijCI04nq9sKEfj2W+Ngn0OHYNIkeO016NQJfs7J4d7kZFJLSri7Y0deCQ3Fy17noLTaNXT6aC5GeQowVhHFAKkicnOj9vIitPegUJ2IULSnWrJ6XbVk9TU++I70xfsqnaxuC4qK4OWXjQ1v9vbw7LPw8MNQbmti1sGDvJmeTqCDA3MjIrjez0/nnrSzNDQoTK72sAIjIPyvEft30XRQqFt5jpGszl6RTc6POVTkVKDsFJ5DPPEZaQQJl+46Wd2aVT8ONCrKqMh65ZWwpaCAu5KS2HrqFKN9fXk3MlIvX9VqaPD0kVLKAeiOMWJIEpGyxu3ixdFB4fyIScjfYElWL8+mcKclWR3qVJWH8Pyzp05Wt1IrVhgjhf374frrjSqsnbueXr5qZ1m+eq9evqpZNHSkcA3wAUYxOwWEAPeIyI+N3dELpYPCxSk5VEL2D0aAOLnqJOYSMzaulmT1KF98r/HFsZP+ZtmanHkc6IwZxvkNR6WYeyzLVwdblq9eopevtnsNDQqJwCgR2Wd5HAasEJHujd7TC6SDQsOZik2cXFMtWX3Ikqzu61a1J8I9zh1lo79htga1HQc6cqTw7xPG8tV8k4npQUHMDArCyVaPDNurhgaFzSISV+2xAjZVv2YtOig0LhGhcFehkYdYkWMkq81g7386We1zlY8+Za4VqH4c6DXXGMHBK7iMR/fv5/8yMohydmZ+VJRevtpONTQovA8EA19h5BRuBA4BKwFE5JtG7e0F0EGhaZVnn5GszrUkqy/zrKrP5BzprJPVLdSZx4E+/jjMnAm/l+rlq+1dQ4PCp/U8LSJye0M61xA6KDQfc4WZ/A35RumN5dnGOdWAU5hTVYDwGuqFjaPeUdvSVD8OtEsXI/fw17Em/pGWyj8PHybAwYG54eGM8/fXAb6d0AXxtEZXklYtWb3aSFbbutnifaWxs9rnGh8cO+pkdUvy++9Gob3t2+Hyy41RREnXAu60LF+9zteXdyMi6OLkZO2uak2soSOFz4CHReSk5bE38IY1RwiVdFBoGUxFJnJX5xqjiBXZlB62JKv7u1UteXWP1cnqlsBkgg8+ME59O3UKHnoInnrGzCf56TxrWb76Umgo9+nlq21aQ4PCVhHpe65r1qCDQssjIhTuLKwq4Je/Pt9IVgfa4/tXo/SGz1U+2HnoZLU1VT8ONCDAKJdx6bhi7ktJ5le9fLXNa2hQ2A4ME5Fcy2Mf4DcR6d3oPb1AOii0fOXZ5eT8VC1ZfbICZW9JVleWAY90sXY3260zjwOdO1fY3TGDaZblqw917kyIkxMTAgPx1cnoNqOhQeFWYCbwteXSjcAcEfmiUXt5EXRQaF3MFWby15/eWV20uwgA53DnqgJ+XkO9sHHQyermZDbDggVGMjonB+69F6Y9U8bjmUksy85GAc42NvwRG0uUiw7gbUFjlLnoCVSen7BaRPY0Yv8umg4KrVtxanFVHiJ3dS5SKkay+qpqyeoOOlndXKofB+rtDROeKmC+/T7KtntCzElGXGbHop498bTTU3+tnV59pLV4pkIjWV05iig7YpTXco91Nwr4jfLFvZ9OVjeH7duNjW///S+gLH8fHMzwxnY69y3hvchIrvPzs2oftYbRQUFrVUSEwh2FVQEif0M+iCVZfY2xJ8L7Sm/s3PU31qYiAhMmwFdfCaBQSpj8aBlbxu9gZ2EhN/n7805EBIEODtbuqnYRdFDQWrWyrDIjWb08m9yfc08nq4dWS1ZH6LnuxrZ+PVxxhVFsT8RyHOgcMwWjDjMnPRVXW1v+GRbG5A4d9Ka3VkYHBa3NMFeYyV+XX1XAr2iPJVkd6Xy6DPhlnjpZ3UjWr4e1a6FbN/jkk9PHgT7+ajHzAxP5PS+P4d7efBAZSaizs7W7q50nqwQFy1nOnwOBGDWTPhSRty1LWr8EugGpwE0ikmsptPc2cA1QBNwmIn/U9xk6KGjFB4urCvjlrrEkq91PJ6t9r/HFIVBPcTSGM48DnTBRiH70OC8W76NChOdDQni4c2fsbHRAbumsFRQ6Ah1F5A+llDuwBRgD3AbkiMjLSqnpgLeIPGk5t2EqRlAYCLwtIgPr+wwdFLTqTIUmcldVS1YftSSr49yrppnc+rrpZHUDVT8O1M4OHp5Rzva/JrEiP4tYd3c+iooi2s3N2t3U6tEipo+UUkuBeZbbMBE5Zgkca0UkSin1geV+vKV9UmW7ut5TBwWtLiLCqe2nqgr45W80ktUOHRyqjiP1Hq6T1Q1x4IAxali2DCIjhRuez+OjzrvJqajgia5deSY4WJ/Z0EJZPSgopboB/wF6AYdExMtyXQG5IuKllFoOvCwiv1ueWwU8KSIJZ7zX3cDdAEFBQf3T0tKavP9a61eWeTpZnfNzDqY8E8pB4fVnL+w72JO5KJPuC7oTOCnQ2l1tdX74wTgOdN8+GDnajNPUgyyxPUyk5cyGofrMhhbHqkFBKeUG/IaxC/obpdTJyqBgeT5XRLzPNyhUp0cK2sUwl59OVh9feJzyo+XGEwr8bvQjdE4oLuF6NdOFKC09fRyo2Qw3TSvit6t3kkYx93TsyCthYXrTWwtSX1Bo0oyQUsoeWAL8u9phPBmWaaPKvMMJy/UjQNdqL+9iuaZpjcrG3gavP3sR9loYFVkVp58QyPoqi00Rm0jol0DaS2kU7SuyXkdbEUdH41zoxES47jr4/GUX1O0DGJ0UyYdHj9Fz0yaWZmVZu5vaeWiyoGCZGvoY2Csi/6z21DJgsuX+ZGBpteu3KsMgIK++fIKmNYbun3bHxtn4z8DG2Yawt8MIeyMMG0cbDs48qAPEBeraFb780jgO1MVZsfTeTgx6cQjuR90Ys2sXN+7ezfHSUmt3U6tHU64+GgL8F9gJmC2XZwIbMY72DALSMJak5liCyDxgBMaS1Cn1TR2Bnj7SGkfGwgwSpyTS/dOaOYWSQyVkLs4k8+tMY1c14NbXDf8b/fG/0V9PMZ1DeblRR2nWLCguFgbfmc+Ga3fg4q54IyyMKXrTm9VYPdHcVHRQ0JpLnQHiJn8CbgzAOUxv3KrL8eMwfTp89hkEdjLjPTWNxIFpXO7txYdRUYTpTW/NTgcFTWtEOkBcnHXrjLMbtm6FyMGlHLlnF+aQQp7r1o1pXbroTW/NSAcFTWsiOkBcGJMJ5s83Tn3LyxOCJ2ZzYPxe+ndw5qOoKGLc3a3dxXZBBwVNawY6QJy/7GwjMHz4oeDhK5jv2Ufh5Ud5olsQzwYH46w3vTUpHRQ0rZnVGiD6GUlqHSBO27LFmFLasAH8Y4rIvG8PEdEm5kdF0cvVlSWZmYzz99dHgTYyHRQ0zYp0gKif2QxffAFPPAGZmYLb6BMU3JaCvZcJO6WwU4rN/fvro0AbkQ4KmtZCVAaIE1+doGBjAaADRKW8PJg9G+bOFWzcTJSPOAquFTj0z+PrG7twnb+/tbvYZuigoGktkA4Qtdu1C266xcTebbYYR+4JfT5M4dvxQfrMhkZitTIXmqbVzSnIia6PdqX/hv4MShtE2BthKHvFwRkH2Ri+kYT+CaS9nEbx/mJrd7VZ9eoFN99oi7IxjgKlXLHr1Y70XL2VF9PSKDObz/ke2sXTIwVNa2H0COL0UaBlZaAUmEyCg7eJ0rtS6DG2gA+jIhmiq69eND19pGmtVEmaJUB83f4CROVRoMOGgbMz3H+/cc0xOp/Sh5O481IPXgkNxUevTLpgOihUk1WUxdyNc6tqrjw44EH8XPyaonua1qjac4AAY5XSp5/CE08KuSeBcen43HmEt3p342+BgbqO0gXQQaGan/b9xDX/vgbB+L2vDrua4aHDifKNIsovilDvUOxsdN13rWWrK0AE3BSA/43+OIe23QCRnW2U6Z4/H+wDyii/P4XLryvn/ahIIvWy1fOig0I1IkKXN7twtOAoTnZOuNm7kVV8us67nY0d4T7hRpCwBIrKn3pEobVE7TVAbNgA990nbNumsBuQg3poH0//KYAng4Jw1HWU6qWDwhmWJi5lzJdjWDphKddFXUducS5J2UkkZSUZP7OTSMxKZF/OPspMZVWv83H2qREouvt1J8o3ijCfMBxsHRrzV9O0i9LeAkRFBbz3Hjz9jFBYLJgnHiLi9hN82CeCYd7e1u5ei6WDwhlEhI1HNjKw88B65yFNZhOpJ1PPChhJWUkcO3X6/B9bZUuId0ito4tA19NznVlFWczbNK/qdTqfoTWl9hQgjh2Dxx+HhQvBrnMJFVOTmXydPa+HheHnoL+wnUkHhSaQX5p/OlBUCxjJ2cmUVJRUtfN09CTSN5Luft1RKL7Y8QWCYKNsWDFpBSPCR1il/1r7UmuA6O9GwI1tK0CsXg333S8kJynU0Ew8HjrIm5d25TZ9oE8NOig0I7OYOZx3uGoKqnrASM9Pr9E2wCWACN8Iwn3CCfMOM376GD99nH2s9BtobV1bDxBlZfDGG/Dc80KZmDHfmsqQO/L58JJIeri6Wrt7LYIOCi1EYVkhH275kEd/eZRbet+Cg50D+3L2sS9nH0cKjtRo6+XkRbhPeM2AYfnZwU1/69EaR1WA+OoEBZvODhD5G/JrPaq0NUhNhYenCcuWKmy6FaKmpTDjOk9mBgW1+9LcOii0IHXlM4rLizmQe4D9ufurAkXl/bSTaZjEVNXWxd7lrEBROcLo6tEVW5v2/X947eLUFiBQgIByVHT/pPUFBoDvv4cHpgqH0xRceZzgaUeYPyiEK33a72hcB4VWrtxUTlpeGvtzzg4YB3IPUGoqrWprb2NPiHdIrSOMbl7dcLRztOJvorUWJWklbAjbAKZqFxV0/XtX/Eb74THIA2XTekarRUXw4ovwyquC2dGE+Y4DTLijgreiwglsh4loHRTaMLOYOZJ/pNYRxv6c/RSUFVS1tVE2dPXoWusII8w7DFcHPd+qnZaxMIOkO5MwF5tR9grnKGeKE4uRCsE+0B6/a/3wHe2L93BvbJ1ax+g0KclIRK9ZrVCRBbg+tp83RgVwZ8eO2LSjKVmrBAWl1CfAKOCEiPSyXPMBvgS6AanATSKSq4x5lLeBa4Ai4DYR+eNcn6GDQv1EhMyizKoAUSNg5O4nqyirRvsObh1qHWGE+4Tj7azXfLdHGQszauQUyk+Wk/NjDllLs8j5IQdTgQkbVxt8rvbBb4wfviN9sfdp2bWIROCrr2DqNDOZGQquPcqAR7L4KDaM3m5u1u5es7BWUBgKnAI+rxYUXgVyRORlpdR0wFtEnlRKXQNMxQgKA4G3RWTguT5DB4WGOVlykv05+2uMLPblGj/PTHx7O3mfHll4h9cYZVTfi6G1H+ZSMyfXniRraRZZS7MoO1oGtuA11Au/0cYowrlby13JlJ8Pzz4rvDMX8ChH3XOAx++0Z1ZIN1zaeCLaatNHSqluwPJqQSEJGCYix5RSHYG1IhKllPrAcj/+zHb1vb8OCk2nqLzISHyfMcLYl7OPtLw0zHK6pr2rvSthPmG1jjC6eHTRie92QMxCwZYCsr4zAkTR7iIAXKNd8Rvth98YP9xi3Frkl4dt2+Du+8xs3mADvU/S6Yk05v+1C9f4+lq7a02mJQWFkyLiZbmvgFwR8VJKLQdeFpHfLc+tAp4UkbP+4iul7gbuBggKCuqflpbWZP3XaldmKiPtZNpZI4zKxHf10iAOtg6EeIVUjTAqRxeViW9dHqRtKtpXRPbSbLKWZpH3vzwwg2OQoxEgRvvhOdQTG/uWU5/IbIYFC+DRv5vJy1MwLp0xj53i3ehQOjm2vcUZLTIoWB7nioj3hQSF6vRIoeUxmU0cKThS6whjf+5+TpWdqmpro2wI8gyqMbKovB/qHaoT321EWWYZ2cuzyfoui9xfcjGXmLHzssNnpA9+o/3wGeGDnXvLqEycnQ1PTBc++UiBXynOUw/w6u0e3Ne5E7YtcJRzsVpSUNDTR+2YiHCi8ESdK6Wyi7NrtO/o1vH0yMK75kopnfhunUyFJnJ+zTFGEd9nUZFdgXJQeF/hbeQhrvPFsaP1v5lv2AC332Nm7w4biM2h98yjfD48mBh3d2t3rVG0pKDwGpBdLdHsIyJPKKVGAg9yOtH8jogMONf766DQtuQW57I/d3+tK6WOFhyt0dbH2afWEUaYT5hOfLcS5goz+evyjUT1d1mUHDBqhrkPdMdvjB/YQurTqVbbTV1RAe++K8x4RiguATXxEA/+3cSL3YNxs2sZI5uLZa3VR/HAMMAPyABmAd8BXwFBQBrGktQcS35hHjACY0nqlHNNHYEOCu1JYVlhjR3f1VdK1Zb4rr5SqvoIQye+WyYRoXB3oTGC+C6LgoTT+2uUnSL4uWCCpwdbJdgfOwYPPWpm8SIb6FSM32OpfHSrP6P9Wm+FY715TWvTykxlpJ5MrTWPcfDkwbMS36HeoVUji0DXQHad2IWPsw+eTp5MGzRNlzNvAX5z/A0pq/m3ySnUCf9x/vjf4I97nHuzB4jVq2HKvSYOpdjCZZlc9WwOHw0JpquTU7P2ozHooKC1WyazifT89FpXSu3P2U9heWGN9r7OvvTw70GIV4hx8z79s7N7Zz3KaCbVd1PbONkQcHMAZell5K7MRSoEx66O+I/zx2+cH56XejZbyY2yMnj1dTPPPQ/lCA63HWLO43ZMC+mMXSs67U0HBU2rhYhw/NRx4ubHcaTgCD7OPoyKGEVqXioHcg9wJP9I1VneYNSVCvIMOh0ovEII9Q6teuzn4qdzGY3ozN3UAOW55WR/n03m4kxyfslBSgWHjg74jfXD/wZ/PC/zxMau6f84p6bC3VNN/LrcFoILCZ9xmPiJnYj18Gjyz24MOihoWj3OPJ61UmlFKYfyDnHw5EEO5h40fla7f2aZEFd71xoBo/ooI8QrBHfHtrFypaWoKKgge4UlQPyQg7nYjL2fPX5jjADhdblXk++FWLZMuPNBM5mHbeHK49zxjyL+GReERwtPROugoGn1ON/jWc9UUFpA6snUOoNG9T0ZYExN1RU0gj2DdQXbBjAVmsj5KYfMJZlkf5+N6ZQJOy87fEf74j/OH+8rm65oX1ERzHrBzD9fV5gdTHjec4gPHnfjpg7+LXbkqIOCpjUzESG7OPt0sDgjaKTlpdVIgCsUnT061wwW1e53cu+k8xnnyVRiIvfXXDIXZ5K9LJuKkxXYutviO8oIED5/9cHWpfH/LZOS4JZ7Ktj8mx1EFvCnZ4/zf9d3oZtzy6v/pIOCprUwJrOJowVH6xxl1JbPCPYKrnNqSuczamcuM5O7OpesJcZeiPKscmxcbPD5qw/+N/jjO9K3UXdTi0D8IjP3TRPyM22wu/YYTz9vZmavGwmBbwAADTRJREFUTti3oES0Dgqa1spcaD7DzcGNbl7dak2Ah3iH4ObQPkpC18dcYSbvP3lkLskk65ssyo6XoRwVPlf74D/OH99rfbH3bpyy3/n58OhTFXz8ni24l9Pl4XS+fMSXS708G+X9G0oHBU1rYy40n+Hn4lfn1FSwV3C7K0woJiFvfR6ZizPJWpJFaXopyk7hPdwbv3FGVVcHv4b/m2zbBhPuKicpwR56n+TGOSf5cERnvOyte+aEDgqa1o6ICFlFWTUDRrXAkXYyjXJzeVX7uvIZ3s7erDm4Bg9HD5RSPDjgwTa5sU/MQsHmAjKXZJK5JNMot2ELXn/2wv8Gf/zG+uHY4eIXAZjN8K+PTTz2JJTk2eBy01HmvmDPlFDrJaJ1UNA0rcqF5jMqDew8kEu7XkqETwThPuFE+EbQ1aNrm0qAiwintp0yAsTiTIqTikGB5xBPY7Pc9X44db24HczZ2XDX42V8u8AB/Erp88QxljwQSLhL8yeidVDQNO28VeYzDuQe4JZvbyGzKNPYg+EVwv7c/RRXFFe1tbexJ9Q7lAjfiKoT+SJ8jaAR5BmEnU3LXq9fHxGhaE8RmYuNEUThTmP3u/tAd6Pcxjh/nEMv/A/6/9YL4+8q58huB2zicnj41WJeHtoRh2ZMROugoGnaRTlzY5+IcLTgKPty9pGSk1JVY6ryflF5UdVr7W3sCfEOqSp9XhksKg9Yam0Boyi5qGqK6dQWI2fj1tcN/xuMAOES5XLe71VRAS/NLee5Z22oKFX43XqM+OddGd7R69wvbgQ6KGiadlEuZGNfZdmQygCRkp1SVWcqJTulRp0pOxs7unl1M0YWPqeDRYRPBN28upFXmse8TfOq2re0fEbxwWKyvskic3Em+RvyAXDt5cr/t3fnMVaVZxzHv493VgcYZivI1qEwUKIhsVo0tipdnBDbgBUj2FabuILRtv7hUpM2mjZCYkyUNI3S1iY0iBqoAoql2rqQFiPEBRELDlsdFGQ2hIHZn/5xzlzu3Lmzwb1z78z8PskJZ3nPuc/MmzsP73nPed/ShcHb1AXnF/Srv+Dzz+GGu5t5c10uTDhF5W9qeObm8ZSkuCNaSUFE0srdOdJ4JEgUca2LT+o+6fK0VMQilBWUcfjEYSDoCH/oyodYdMEiphZNzbgnpZqqm4IEse4ox7YcA4f8GfnREV1HXdj33NQvv9rOTUvbqdubQ84VtSx/rJ1fXpS6jmglBRHJWJ0z8nVJFLWfsH73eprbm7uU7ZzCNbZ10dnCmFo0lbys9A5j3Xy4mZoXa6hZV0P96/XQDnnlecFTTAtLGTNnTI8jura0wL3Lmvn9smw6zJl2xxFefGgsF4zt/22p/lJSEJEhp7M/Y9U1q5hePL1LH0ZnAmloaoiWN4zJhZO73ZLqnGApP3twn/JpqWmhdkMwYF/9a/V4q5M7KZfSa8MRXS8rxCLdE8S+/c61S5v4YHM+Vt7Iz5Z9yZPXjyM3iR3RSgoiMuT01Z/h7tSdqut2O6pziZ/ze9KYSQk7vacVTaMgpyClP0trQyu1L4Ujuv49GPI7e1w2ZT8KbjEVXtl9yO9Vf2vlzp87jYdyGD2vhr88nsXCmcnpiFZSEJERp/5UfbeWRef60ZNHu5SdMHpCjwkj2UOetx1vo25TOKLry7V0nOwgqyQrGPJ7YRlF3yvinJwgQZw8Cbf++hRrVuRCbgdz7qlh/YPFjM8/u34VJQURkRjHmo51SRhV9VXRTvAjjUe6lB0/avzpW1FxSWNM7tlNqtN+sp26zXXBiK4ba2k/3k6kMELp/DBBVBYRyY/wwa52rrmtmQP/OZfIjBNcd2sLTXURllxdwLzLB/5or5KCiEg/HW8+zt76vV2elOp8tPaz4591KVt2btnpJBH38t7YvIHd6ulo7qD+tWDI75r1NbTVtxEZFaH4B8GIrsXzSnhiXSv3/yJC+7EsMCCngzWvtLL4OwPrYFdSEBFJgsaWxuh83/HvYlR/Wd2lbEl+yekkURTzpFRJBcX5xb1+TkdrBw2vNwQjur5QQ+vRVs7JD4b8XtZ0Hs9sKgYMrIMbHzzFqt8NrE+kt6SQUa8Umtk84AkgAvzJ3ZenOSQRkaiCnAJmj5vN7HGzux071XqqS8Lo7Md46+BbrN6xust4UkV5RV0ep41NGCX5JZyTfQ7FlcUUVxYz4w8zaNjSQM264F2Iyz5vYR1jacXIdue7fhxIXkd5xrQUzCwC7AGuAqqBbcAN7r6rp3PUUhCRoaCprYn99fu7PSFVVVfFwWMH6fCOaNnC3MKECWN68XTK8svYWLyRp76xicPkMZ4m7njvahY0LBhQPEOlpTAHqHL3fQBm9iywAOgxKYiIDAV5WXnMKpvFrLJZ3Y41tzVzoOFAt6ektn+2nbW71tLu7dGyo3NGU3J/CQdaDoR3j4wbfzw3qbFmUlKYCHwas10NXBJfyMxuB24HmDJlyuBEJiKSIrlZucwsncnM0pndjrW2t3ZJGFV1Veyp3UP13mraaKM8t5zFty5OajyZlBT6xd1XAishuH2U5nBERFImO5JNRUkFFSUVXfZ3vu39+MLHkz4+UiYlhUPA5JjtSeE+ERGJMX/mfLbespVLJna7mXLWBm9Wh75tAyrMbKqZ5QCLgQ1pjklEJOOYGZdOujQlo6hmTEvB3dvM7C5gM8EjqU+7+0dpDktEZETJmKQA4O6bgE3pjkNEZKTKpNtHIiKSZkoKIiISpaQgIiJRSgoiIhKVMWMfnQkzOwocjNtdCBzrx+l9levteE/HEu1PtK8UqOlHjMnW399NKq6Trnrpb12lq04SxTJY19F3pWfD+btSCIx197KEZ7n7sFqAlcko19vxno4l2t/Dvu2Z/LtJxXXSVS/9rat01Uk660Xflcyrk8Gol77OG463jzYmqVxvx3s6lmh/f+MZDMmK5Uyuk656GUhdpUu66kXflZ4N5+9Kr+cN6dtHQ5WZbfcehq2V9FCdZCbVy+Abji2FoWBlugOQblQnmUn1MsjUUhARkSi1FEREJEpJQUREopQUREQkSkkhw5jZLDN70szWmtnSdMcjYGbXmNkfzew5M6tMdzwCZvY1M/uzma1NdyzDjZJCEpnZ02b2hZntjNs/z8x2m1mVmT3Q2zXc/WN3XwJcD3wrlfGOBEmqkxfd/TZgCbAolfGOBEmqk33ufktqIx2Z9PRREpnZFcAJYJW7XxDuiwB7gKuAaoIZ5m4gmEhoWdwlbnb3L8xsPrAU+Ku7PzNY8Q9HyaqT8LzHgNXu/u4ghT8sJblO1rr7dYMV+0iQUZPsDHXu/paZlcftngNUufs+ADN7Fljg7suAH/ZwnQ3ABjN7GVBSOAvJqBML5jxcDryihHD2kvU9kdTQ7aPUmwh8GrNdHe5LyMzmmtkKM3sKzUKXKgOqE+Bu4PvAdWa2JJWBjWAD/Z6UmNmTwIVm9qtUBzeSqKWQYdz9DeCNNIchMdx9BbAi3XHIae5eS9DHI0mmlkLqHQImx2xPCvdJ+qhOMo/qJEMoKaTeNqDCzKaaWQ6wGNiQ5phGOtVJ5lGdZAglhSQyszXAVmCmmVWb2S3u3gbcBWwGPgaed/eP0hnnSKI6yTyqk8ymR1JFRCRKLQUREYlSUhARkSglBRERiVJSEBGRKCUFERGJUlIQEZEoJQUZ8czsRB/Hx5rZnTHbEzSOvwxXek9Bhr1wlFNz944ejp9w91G9nF8OvNQ5zPNgio+9r58l5rys8IUwkQFRS0GGJTMrDydsWQXsBCab2b1mts3MdpjZwwnOGWVm/zSzd83sQzNbEB5aDkwzs/fN7NHw2jvDc942s/NjrvGGmV1sZgXhZDLvmNl7MdeK/8xuMSWI/fIEP8ujZrYzjHNReN5cM9tiZhuAXUn7ZcrI4u5atAy7BSgHOoBLw+1KYCVgBP8Zegm4Ijx2Ivw3CxgTrpcCVWH5cmBn3LV3huv3AA+H6+cBu8P1R4CfhutjCSaQKYiLMWFMCWKP314IvEowAc044H/hZ88FGoGp6f79axm6i1oKMpwddPe3w/XKcHkPeBf4OlARV96AR8xsB/AawXj+4/r4jOeBzpm/rgc6+xoqgQfM7H2CodDzgClx5/YWU2zs8dvfBta4e7u7HwHeBL4ZHnvH3ff3EbNIjzSfggxnjTHrBixz96d6Kf8ToAy4yN1bzewAwR/zHrn7ITOrNbPZBPM3d47xb8BCd9/dy+kJYwr7MBrjysZv96S/5UQSUktBRorNwM1mNgrAzCaa2VfiyhQCX4QJ4TvAV8P9x4HRvVz7OeA+oNDdd8R83t1hxzBmduEZxpTIFmCRmUXMrIzgltM7/ThPpE9KCjIiuPs/COa73mpmHxLc5on/Q78auDg8fhPw3/DcWuDfYcfuowkuv5Zg/P/nY/b9FsgGdpjZR+H2mcSUyAvADuAD4F/Afe5+uB/nifRJj6SKiEiUWgoiIhKlpCAiIlFKCiIiEqWkICIiUUoKIiISpaQgIiJRSgoiIhKlpCAiIlH/B1qwYU6u0cIvAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "ACVI_Higher_dimension_bg.ipynb",
      "provenance": [],
      "toc_visible": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}