test_1d_data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "import ot\n",
    "\n",
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "from src import sinkhorn\n",
    "from dataset_generators import empirical_distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('font', family = 'serif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiments on 1d data\n",
    "---\n",
    "We are generating the data on the real line $\\mathbb{R}$ from the normal distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 20\n",
    "\n",
    "x_i = generate_1d_normal(n_samples)\n",
    "x_f = generate_1d_normal(n_samples)\n",
    "\n",
    "C = cdist(x_i, x_f)\n",
    "C = jnp.asarray(C)\n",
    "\n",
    "a = jnp.ones(n_samples) / n_samples\n",
    "b = jnp.ones(n_samples) / n_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us have a look at the effect of regularization on the optimal assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4sAAAEkCAYAAACGxqRxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYOklEQVR4nO3de4ydd53f8e9vzsz47jhxHIdJDElwAk4ImcImkMB2U4o0QNku2y2r7AILqlDLVtVWaosqLaqKqq5ahLZqV+q2VNWKZdm2oUoaupR62G7KEhSyIcETTBzI/ebBiWPwdTK2zzlP/5iJavKzx/t9MhfP+PWSLMKceZ/nOXN5PJ85Y7s0TRMAAABwqoGlPgEAAADOPcYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIvMu1LKr5dSnimlHCul3FVKuWiO121mX+/o7K//vJjnCpCRvL6NllIeLKVMzf7v6Cm3fbaUcvKUa9/RUspVi/IgAF5lHq9tpZTyuVLKgdlfnyullFNu93XfMmMsrhCllMGlPoeIiFLKdRHxhYj4WERsjYipiPj9s2Q3NE2zfvbXJxf6HIHlZTle30opwxHx1Yj4ckRcGBF/GBFfnX35K24/5dq3vmmaJxf0AQDnlBV6bfu7EfGhiLghIt4aEb8YEX/vVXfj675lxFhcxkopT5dS/mkp5fsRcayUMlhKeWcp5d5SysFSykOllFtPef0rSynfKqUcKaX8n1LKvy+lfHmeT+sjEfEnTdN8q2maoxHxzyLib5VSNszzcYAVbAVc326NiMGI+LdN0xxvmub3IqJExHvm+ZyAZeQ8uLZ9PCJ+t2ma55um2RsRvxsRn5jn82URGYvL369FxN+IiE0x892g/xUR/zIiLoqIfxIRd5RStsy+7n+JiPsjYnNEfDZmvoN0WqWU189etM7069fPkF4XEQ+98n+apnkiIk5ExDVzPIZvlVL2lVLuLKVccdZHDJwvlvP17bqI+H7TNM0pL/v+7Mtf8YullJ+UUh4upfzmmc4XWHFW8rXtZ+5r9r9Pve5F+LpvWTknnv7mNfm9pmmei4gopXw0Ir7eNM3XZ2/701LKAxHxgVLK/42IGyPirzdNcyIivl1K+Z9nutOmaZ6NmYtY1vqIOPSqlx2KiDM9s/gLEXFfRKyNmQvl10opo03TdFscG1hZlvP17Wyv+5WI+E8R8UJEvCNmvjg82DTNf21xXsDyspKvba++/VBErC+llNmB6eu+ZcYzi8vfc6f89xsi4sOnfhcpIt4dEa+LiJGI+EnTNFNnaOfL0YjY+KqXbYyII6d75dkfeTjRNM3BiPiHEXFlROxYgPMClp/lfH2b83WbptnTNM1k0zS9pmnujYh/FxF/e57PFzg3rdhr22lu3xgRR195JtLXfcuPsbj8nfpjAM9FxB81TbPplF/rmqb51xHx44i4qJSy9pTX33amO539UYajc/z6yBnSh2PmDzW/cj9XRcSqiHg08XjKWV8LOB8s5+vbwxHx1lLKqdezt86+/EyP1bUPzg8r+dr2M/c1+99nuu5FuPad84zFleXLMfNnYMZKKZ1SyupSyq2llMubpnkmIh6IiM+WUoZLKTfHzN9QdVpN0zz7qr+l79W//vgM6R/PnsPPl1LWRcS/iIg7m6apvjtVSrmuzPz1y51SyvqY+UPQeyPikdf4dgBWnmV1fYuIb0ZELyJ+q5SyqpTyD2ZffndERCnll0opF5YZN0XEb8XM3zAInF9W1LUtIr4UEf+olHJZKWUkIv5xRHwxwtd9y5WxuILM/vz7L0XEb0fE/pj5btWn4/+/nz8SETdHxIGY+Tnx2yPi+Dyfw8MR8amYufC8GDM/w/73X7m9lPK/Sym/Pft/t86ew+GIeDIiroiIDzZNc3I+zwlY/pbb9W32zxd9KCJ+IyIORsTfiYgPzb48IuK2iHg8Zn5060sR8bmmaf5wPs8XOPetwGvbFyLiTyJid0T8IGb+8p4vzN7m675lqPzsX2bE+aSUcntE/LBpmn++1OcCMJ9c34CVyLWNxeaZxfNIKeXGUsobSykDpZT3xcx3su5a4tMCeM1c34CVyLWNpeafzji/XBoRd8bMv9XzfET8ZtM0u5b2lADmhesbsBK5trGk/BgqAAAAFT+GCgAAQGXOH0P9+b/5+fTTjmvGJ9InMXD1FemmeWZvuomIOH7zm9PNqm/P9c/DnF7/eIu/qKrFs7xlaDjdTL/3hrO/0quse+SFdBMRESe76aRZv/bsr/QqvceeSjedLZvTTRnIf3/l5esvTzcREUPfeCAfDXTyTb+XTv60/99X1L+JdPW/+jfpT74rPvOdhTiV+XPT9fnm/t3zfx5L6KefuDndXPjFc/z9eg4bvHRrq667r+XvL4tgpV3rDuy9LH2tu23bLQtxKsyjgbX5r5v6U1MLcCYsV3Nd6zyzCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoDM5145rxifQdHvrw29LNpru+n276U1PpJiJi+M93p5uT73pLuln1zIF00xw8nG6O3HpNuln3tV3pplmzOt1ERMTQnB9ipz/WwUP55p3591F56oV0092Xb4b3v5RuIiLK2rXppnnTlfkD7Xk836wwV3zmO+nmjd/Nf048ceN0umnt/vy1LgY6+abfyzctnBj7uXRz4Rfz71faa3N9ZHHdtu2WdDM+OZFuxkZG0w3ttf2aGP4yPLMIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAICKsQgAAEDFWAQAAKBiLAIAAFAxFgEAAKgMznXjwNVXpO9w013fTzfTf/W6dLPqzx5KNxER5do3ppvhh55KN91Dh9NNNP10sv4bP0g35Zor001v3XC6iYgYON5NNyc2r003Q/fsTjfNldvSzcDR9emmvO6SdBMR0Xv0ifyxdv8o3TTd/PuIiCdunE43L4/nP/fWjOWvPxERUUq+6ffaHWsRDI8/sCjHGVi3rlXXn5pKN0d+9R3pZsPt96UbeC3GRkbTzfjkxKIcB1h4nlkEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAICKsQgAAEDFWAQAAKBiLAIAAFAxFgEAAKgYiwAAAFQG57qxeWZv+g77U1PpZtWfPZRudj5zf7qJiHj/9jkf8mn1WjymxdI/dizdDDz9fLopx4+nm4iIfr9JN0OdTrppTp5IN73Hn043L33ypnSz9St70k1bTbe7aMcib83YU+lmfHKi1bHGRkZbdee7NtfUtjbcft+iHQsWU5vrj2sdnJs8swgAAEDFWAQAAKBiLAIAAFAxFgEAAKgYiwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVAbnuvH4zW9O3+Hwn+9ON+XaN6ab92+f89TP6OVbr0s3a775cLrpv/xyuommSSdl1ap0c/xdO9LN6icPpJuIiHKym26a9WvTTe+Rx9PN4NYt6WbrV59INyevvyrdREQM3LOrRdTJN/1evllpbro+39yfv9ZFKelkbGQ0f5yIePyP/kq62f6xFh9z57AT77sx3Qzv/O4CnMn5YWD16lZdf3p6ns+E5ajttW58cmLRjgXnI88sAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqg3PduOrbD6fv8OS73pJuhh96Kt30pqbSTUTEmm/mH9OBD9+Qbi7e+US6aY4eSzfT796RblbfsyfdxIb1+SYiYng43+zbn046V1+Zbpp0EdF7LP+xOnjocIsjRTSrVqWbzmWvSze9Z59PNyvO/bvzzUAn3/R7+aal7R/blW7+4zPfTjefesO7000bnY0b083wzu8uwJlwJv3p6aU+Bc5DYyOj6WZ8cmJRjgMrgWcWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAACVwblu7B8/nr7DVc8cSDfdQ4fTTVv9l19ONxfvfCLdPPuJ7elmy8SJdLN2z750023xNign8uc2E+a/H9H0evnmmm3pprPrR+km+vlz6x/v548TEZ1rr0k305euTzer++3O77zX4mPhXPepN7w73Xx97/fSzQcue1u66R1evN8ngJVtbGQ03YxPTizKceBc45lFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAICKsQgAAEBlcM5bmyZ9h83Bw/mzaPr5pq02j+nosXSzZeJEutn714bSzfa/OJJuWr0Ner38cSIiSv5YbT4ehn7803TT63bTTSst3t4REQevvzDdbPre/nQzfdWWdAOv+MBlb0s345MT6WZsZDTdAMyXNtcg1zpWAs8sAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqg3PdWIaG03d45NZr0s36b/wg3fSPHUs3ERFl1ap0M/3uHelm7Z596Wb7XxxJNz99f/7cLrhzV7oZWLs23URExOCcH2Kn1UxPp5ve5g3pptO7NN10n9+bbsrgULqJiLjgrol002y/It0M3bM73RBxYuzn0s3w+AMLcCbzp7NxY7rpHT6cbsZGRtPNlns3pZv9txxMNwDzpc21bnxyYlGOA39ZnlkEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAICKsQgAAEDFWAQAAKBiLAIAAFAxFgEAAKgYiwAAAFQG57px+r03pO9w3dd2pZtyzZXpZuDp59NNRMTxd+1IN6vv2ZNuui+/nG6iadLJBXfm397T731rulm358V0ExERJ7v55qIL0kl/1w/TTbNlc7oZ3HpJupm+9vJ0ExExePeD+WjPY/mm38s3K8xPP3Fzurnwi99ZgDOZPyfed2O6Gd753QU4k/mx/5aD6ea5z9ySbrb9zr3phhmdTflrd0RE7+CheT4TzmRg7dp005+aWoAz4UzGRkbTzX97Ln/dum1b/vrI+ckziwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAyuBcN6575IX0HTZrVqeb3rrhdFOOH083ERGrnzyQjzasTyflxIl00/R66WZg7dp0s27Pi+nm2JsvSTcREUNT3XRTuk266TwzlG762/KPqex5Mt2sfiz/eRQR0S0l3ZShOT+lT6s5nv+4W2ku/OJ3FuU4A+vWpZv+sWOtjjW887utupVk2+/cm252PJj/HIqIePS216ebMjWdbrp7J9NNDHTyTT9/XegdPJQ/DouqPzW11KfAArht2y3pZnxyotWxxkZGW3UsX55ZBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqg3PeerKbv8ehue/ydAaO54/T7zfpJiKitHlMw8MtDtRih5cWj2kw//Zu834dmmrxdouIExuH0s2ayWP5A3U66aS7If9+Her3002cPJlvWiqlpJt2n0m00Z+aWupT4Cweve31rbpfuHN3urn7ho2tjpXV2XxRuunt378AZwKcK8ZGRlt145MTi3Yszg2eWQQAAKBiLAIAAFAxFgEAAKgYiwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVAbnurFZvzZ9h83BQ+nmxOb8cYY6nXQT0e4xxb79+eP0evnjNP18Mj2dP85FF6ST0m3yx4mINZPH0s3xi9ekm1XdbroZfuFoumlafNy1+piLiNiXf5s3vfzHEIvnyK++I91suP2+BTgTzqRMtbimRsTdN2xMN5N3vCndjPzynnTT25//PQzgdMZGRtPN+OTEohyHheGZRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAZXCuG3uPPZW+w+adb0k3Q/fszh/n5Il0ExHRe+TxdNO5+sp001yzLd0M/fin6aa3eUO66e/6YbrpPDOUbmbCTjpZ1e2mm5/82tvSzcV37Uk3/WPH0k08+Wy+iYjO5ovSzaH3XJ1uNu7Mvx1oZ8Pt9y31KXAW3b2Ti3askV/Of+79yiMvpps7dlySbs515e3XpZvmwYcX4EyAsxkbGU0345MTi3Iczs4ziwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAyuBcN3a2bE7fYXnqhXTTXLkt3fQefzrdREQMbt2SbpoWx+ns+lG66XW7+eP0Lk03TYv3a3/bJekmIqK7YTjdDL9wNN1cfNeedHNwbEe62XjXrnQzsGF9uomI6B34SbrZuDP/dugfzb+9YcUa6LTKOpsvSje9/fvTzR078tfizz99X7r59BXvTDeLqXnw4aU+BWABjY2MppvxyYlFOc75xjOLAAAAVIxFAAAAKsYiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAIDK4Fw3loH8luzueyHdDBxdn25e+uRN6SYiYutXn0g3vceeyh+o38s3LXSf35tuBrdekm7KnifTTUTEUL+fbppOJ930jx1LNxvv2pVuNt29Lt0c/vjGdBMRES8dSCf9I0faHQuY0fLa3du/f55PZP58+op3ppvxyYl0MzYymm4A5kuba5Br3dl5ZhEAAICKsQgAAEDFWAQAAKBiLAIAAFAxFgEAAKgYiwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUBmc68aXr788fYfD+19KN+V1l6SbrV/Zk24iIk5ef1W6GTx0ON30j/fTTTRNOimDQ+lm+tr8+3X1Yy+km4iIOHkynTTr1+aP8+Sz6WRgw/p0c/jjG9PNgZsvTTcREZsefyoflZJvWnzcrTSDl25NN919LT8nFsnA6tXppj89vQBnsnQ6my5IN72DhxbgTJafsZHRdPMbP3qu1bG+9KZtrTqA16rNtW58cmLRjnUu8MwiAAAAFWMRAACAirEIAABAxVgEAACgYiwCAABQMRYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAICKsQgAAEDFWAQAAKBiLAIAAFAZnOvGoW88kL7DsnZtuuk9+kS6aWvgnl3pplm1Kt10rr0m3Ry8/sJ0c8FdE+lm8O4H0023lHTT2r4mnXQ2X5Ruegd+km7ipQPpZNPjT+WPExGP/oeb0s2Ozz6dbp7/6PZ0s9J0972w1Kcw7/rT00t9Ckuud/DQUp/CeeVLb9rWqvva3vzvSR+87O35Aw108g3Aq4yNjLbqxicnFu1Y88kziwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACAyuCctw500nfYvOnKdFN2/yh/nG433UREq8fUuex16Wb60vXpZtP39qebZvsV6Sb2PJZOytDcHypn7EpJN02vn24OvefqdLNx55500z9yJN1Ei7dBRMSOzz6dbl784BvTzeV/8Ei6ic/nE+C1K2+/Lt00Dz68AGcyfz542dvTza888mK6uWPHJekGYL6MjYymm/HJiUU5zlw8swgAAEDFWAQAAKBiLAIAAFAxFgEAAKgYiwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqAzOeWu/l7/HPY+nk6bbzR+nrRaPqffs8+lmdb+fbqav2pJuhu7ZnW7avA2a4y0+FiKiaVXlbdy5J930jx5dgDM5jabdW+H5j25PN5f/wSPpZuB/rE43wNJoHnx4qU9h/g100skdOy5JN1/b+2C6AVhKYyOj6WZ8cmJez8EziwAAAFSMRQAAACrGIgAAABVjEQAAgIqxCAAAQMVYBAAAoGIsAgAAUDEWAQAAqBiLAAAAVIxFAAAAKsYiAAAAFWMRAACASmmaZqnPAQAAgHOMZxYBAACoGIsAAABUjEUAAAAqxiIAAAAVYxEAAICKsQgAAEDl/wEkycOAaHy7jQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize = (16, 6))\n",
    "\n",
    "for idx, reg in enumerate([5e-1, 5e-2, 5e-3]):\n",
    "    G = sinkhorn.sinkhorn_knopp(C, reg = reg, a = a, b = b)\n",
    "    ax[idx].imshow(G)\n",
    "    ax[idx].axis('off')\n",
    "    \n",
    "    ax[idx].set_title(f'reg = {reg}')\n",
    "    \n",
    "plt.savefig(\"results/1d_with_different_regularization.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a larger dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 1024\n",
    "\n",
    "x_i = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "x_f = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "\n",
    "C = cdist(x_i, x_f)\n",
    "C /= C.max()\n",
    "\n",
    "a = np.ones(n_samples) / n_samples\n",
    "b = np.ones(n_samples) / n_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us measure the dependence of the time elapsed on the regularization parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.42 s ± 167 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.85 s ± 730 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.21 s ± 89.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "2.04 s ± 548 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.82 s ± 244 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.4 s ± 237 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.6 s ± 122 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.31 s ± 165 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "688 ms ± 35 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "364 ms ± 41.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "205 ms ± 31.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "142 ms ± 14.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "91.9 ms ± 3.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "98.9 ms ± 8.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "83.4 ms ± 6.25 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "77.6 ms ± 2.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "75.2 ms ± 2.95 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "85 ms ± 6.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "80.5 ms ± 4.68 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "82.3 ms ± 4.25 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "regs = np.logspace(-10, 0, num = 20, base = 2)\n",
    "\n",
    "tm1, tm1_err = [], []\n",
    "\n",
    "for reg in regs:\n",
    "    t = %timeit -o sinkhorn.sinkhorn_knopp(C, reg = reg, a = a, b = b)\n",
    "    tm1.append(t.average)\n",
    "    tm1_err.append(t.stdev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9c8c4070d0>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGCCAYAAAAMp6wpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABG1klEQVR4nO3dd3hUZdrH8e9NJyABREAlJCoWcAVWomJ5BQuyimBb3dWoiKzYu66uKNaoa+8FZcHCimXVFXsDZVWUsCq2FRsgIEgXCD3P+8czIWGYSWaSmTlTfp/rOlcypzznnvGYuXmqOecQERERyUUNgg5AREREJChKhERERCRnKRESERGRnKVESERERHKWEiERERHJWUqEREREJGc1CjqAdNWuXTtXVFQUdBgiIiKSAFOnTl3onNsqfL8SoSiKioooKysLOgwRERFJADObGWm/msbCmNlAMxu5bNmyoEMRERGRJFMiFMY5N945Nyw/Pz/oUERERCTJlAiJiIhIzlIiJCIiIjlLiZCIiIjkLCVCIiIikrM0fF5EJEutXr2aBQsWsHr1atavXx90OCJJ0bhxY9q3b0+rVq3qdL0SIRGRLLRs2TLmz5/PVlttRceOHWnUqBFmFnRYIgnlnGPVqlXMmTMHoE7JkJrGRESy0MKFC+nUqRNt2rShcePGSoIkK5kZeXl5bLvttvz66691KkOJkIhIFlq7di3NmzcPOgyRlGjevDnr1q2r07VKhEREspRqgSRX1OdZVyIkIiIx69vXbyLZQomQpA39gRURkVRTIpRFlEiIiCTXQw89RLdu3TAzxowZA0BZWRk9e/akSZMmnHLKKYHGF4sRI0bQpUsXzIyJEydGPe8vf/kLnTt3xsyYMWNGyuJLNSVCIiIiUdx1110UFRVRXl4OwBlnnMGrr766yTnFxcV89tlnbLPNNimLa+XKlRQWFnL33XfHfe11113Ho48+Wut5jz76KNddd11dwssoSoRERESiaNu2LZ07d6Zhw4ZBh7KJRo0a0blzZ7bccsugQ8l4mlBRRERiMnYsTJ4Ma9ZAURGUlkJJSdBRJdfJJ5/MySefHHQYm2natCmTJk0KOoysoBohERGp1dixMGyYT4IAZs70r8eOTc39V61aRffu3TEztt56643Jydy5c+nZsyeNGzemR48eTJgwgYEDB7L77rvTo0cP9tprr82asg477DA6duyImTFlyhT69+9PUVERvXv35quvvtp4XmlpaUx9aaJZvXo1l19+Ob169WL33Xene/funHHGGSxdurTWa+fMmcOf/vQnunfvzu9//3t69+7NzTffDMCiRYvo2bMnLVu2pG+1jqGxvq9Ixo4dS1FREXl5efTs2ZMJEyZscvyHH35g0KBB7Ljjjuy22268//77m5Xx0ksvsccee7DjjjtSWFjIkCFDNpnksHqfo3fffZeBAwfStWtXzIwXX3yRnj170rZtW4qKinjttdc48MAD6dSpE/369WP27Nm1fmZ15pzTFmHr1auXyzR9+vgtU2V6/CLp5Ouvv456rPL/terb/ff7YytXRj6+5ZbOweZb06b++Lhx/vpZsyJf/9JL/vj//le/99W1a1c3YMCATfb98MMPrvJv9umnn+6uuOIKV1FR4Zxz7oMPPnDNmzd3U6ZM2eSaq6++2gHuvPPOcxs2bHDr1q1z+++/v9t77703OW/ChAkOcBMmTNi476effnKAGz169CbnFhYWusGDB298/csvv7j27du7WbNmOeecW7VqlTv22GPdwIEDa32fBx10kDvttNM2vo+XX37Z+a/sKn369HF9wv5o1vV9zZ8/3/Xq1cu9++67m5w3evRoB7g//vGPrry83Dnn3AknnOAKCgrcunXrNp43btw417BhQ/fiiy9ufK/9+/d3Xbt2dStWrNisvAEDBrglS5Y455wbNGiQe+GFF5xzzg0ePNi1atXKXXXVVc4555YvX+522mknd/zxx9f6mdX0zDvnHFDmInzfq0ZIRERqtWhR5P2VNUSpMnjwYN544w3mzZu3cd+YMWM21hBdccUVXHnllRsn2Ntnn33o3r07o0aNiljekCFDaNCgAY0aNWLgwIF89NFHrF27NiGxtmvXjg8//JCCggIAmjVrxqmnnsr48eOZP39+jddOnjyZoqKije9jwIABXHHFFTHfO573tWjRIgYMGMCNN97IAQccEPGckpKSjTOVH3XUUfz888/8+OOPgK9QufTSSznggAM44ogjNr7XW265hW+++YaHH354s/IGDx5M69atAXjkkUc46KCDNh5bvnw5F1xwAQAtW7akX79+daqRi5X6CImI5JiavlPy8iIfLyryzWHhCgs3Pb+goObyd945phCjOvHEE7niiisYO3YsF198Mc45xo0bx4cffghAixYtuPLKK5k4cSLr1q2jQYMGfP/99+Tn50csb6eddtr4e9u2bQGYP3/+xuSlPho1asT06dM566yzmDNnDo0aNWLFihUA/Pjjj3To0CHqtfvvvz/XXnstP//8MyeddBK9e/emtLQ05nvH+r6WLFnCwQcfzO9//3sOOeSQuMvbaaed+Pbbb/n555858cQTN7mme/fuNGvWjLfeeouLLrpok2Ndu3bd+Hv79u03OdauXbuN96i8X22JY32oRkhERGpVWuqTpOry8vz+VNp22205+OCDeeyxxwCYOHEi3bp1o127dlRUVDBw4EBef/11/v3vfzNt2jQ+++wziouLWROl6iqv2ptq0MB/JW7YsCEhsb722msMGDCAQw89lC+++ILPPvts47D1aPFUeu655xg+fDivvvoq++67L9ttt11MQ94rxfq+zjrrLDp06MATTzzBp59+WqfyFi5cCECbNm02u65NmzYbj1fXsmXLmO5Veb+Kioqo59eXEiEREalVSQmMHAlNm/rXhYX+dRCjxgYPHswXX3zBp59+ypgxYxg8eDAA33//PR999BFDhw6lU6dOqQ8szBNPPEHLli254IIL4l4LKy8vjxEjRjBjxgzeeecdCgsLOe2003j77bcTGuM111zDiy++yPbbb8/QoUNZv3593GW0a9cOgMWLF292bMmSJRuPpyslQiIiEpOSEujdG/r0gRkzghs6f9RRR9GqVSvuv/9+3nvvPQ4//HCgqpYlPOmo3p8oldasWbOx9iTeWI4//njAv5cDDzyQF198EYBp06YlNMadd96ZZs2aMWrUKD777DNuueWWOpVRUFDAlClTNtn/xRdfsHr1avr165eocJNCiZCIiGSU5s2bc9xxxzFq1CgGDRpE48aNAdhll13YfvvtGT16NEuWLAHg2Wef5dtvvw0kzgEDBrBs2bKNTVrLly/nrrvuiunacePG8fzzz298/Z///IeGDRuy//77JyNU9ttvP84++2yuu+46/ve//8V1rZlx6623MmHChI0J2+rVq7nsssvYZZddOP3005MQcQJFGkqmTcPng5Dp8Yukk9qGEtdVuvx/OmnSJAe4srKyTfZ/+eWX7oADDnAdOnRwffr0cRdccIHr1auXa9GihevRo4dbs2aNO+GEE1yHDh0c4Hr06OGmTZvmbrrpJldQUOAA17VrV/fss8+6G264we2www4OcDvssIO79tpr3YMPPui6du3qAFdQUOBOP/10N2XKFNejRw/XuHFj16ZNm41D+SsqKtz111/vCgoK3O9+9zt3yCGHuBEjRmws7+677476/v7+97+7PfbYw3Xv3t11797d7bnnnhuHpi9cuND16NHDtWjRYuP7mj9/fszv64477tjkfd1xxx3u3XffdTvttJMDXIcOHdxJJ53kLrnkkk2uff31191jjz22ybXV38OLL77oevXq5bp06eIKCgrc4MGD3fz58zceDy/vlFNO2eQ99+3b17Vp08Y1btzY9ejRwy1cuNCde+65m7ynDz74IOpnVtfh8+aPSbji4mJXVlYWdBhxqZxXK4mjDJMq0+MXSSfffPPNJiNzEkX/n0q6qu2ZN7Opzrni8P0aPi8iIjFTAiTZRn2EREREJGcpEZK0ULmY43vv+YnbUrV+kYiI5DYlQhK4oBdzFBGR3KVESAI3fDiUl2+6r7zc70+kvn2rOnqKiIiAEiFJA7NmRd8fYaJSEYmRRgVLrqjPs65ESALXuXPk/c7BNtvA8cfD229DEpeaEck6TZo0YdWqVUGHIZISq1at2jixZryUCIUxs4FmNnLZsmVBh5IzSkshbBb6jYs5nnYavPEG9OsHO+wA778fTIwimaZdu3bMnj2bxYsXs27dOtUOSVZyzlFeXs6cOXM2W8U+VppHKIxzbjwwvri4+LSgY8kV7dr52p5GjWD9er+YY2lp1TpGt94KL7wA//gHbLed3zdpEsybB4MGVS0CKSJV8vPzadq0KQsWLGDRokV1WkxTJBM0btyYDh060KpVqzpdr0RIArdoEeyxBzRr5muGwidsa9bMN4+F1iAE/KrXTz4JW24JJ50EQ4fC736X0rBF0l6zZs0oKCgIOgyRtKamsSyRyfPwnHACfPzx5s1jNRkzBl591Y8Cu/9+2G03X46IiEg8lAhlgUydh2fFCh9jRQWYxXdtw4Zw6KHw3HMwZw7cfjv07++PrVoFZ5wB//mP73ANmZ0oiohI8mjR1SgyadHVoiKf/IQrLIQZM1IdTeyuuQauvRbKyqBXr8Qt5vjJJ3DQQT7R2nln6NkTXnrJJ0iV8vJ881plPyQREclu0RZdVY1QFqhpHp50NXeu7wR93HE+CUqkPfeEX37xnavbtYOnn940CYLkTNgoIiKZR4lQFojWFzLa/DzpYMQIWLcObropOeW3bAlDhvjmsWjNbumcKIqISGooEcpwX3zhazfCh5A3bOiHoKejadN8bc2558L22yf/ftESQg2mERERJUIZ7P334f/+D5o08f1tKpOhNm1gw4aqztPp5rfffPNVqpqmSkt9n6BwLVv6WEREJHcpEcpQL7wAhxwCW28NH34Il18OvXtDnz6wcCH8+c/QoUPQUUa2335+BFfbtqm5X0mJ7xhdmSgWFvpRZdOn+89r3rzUxCEiIulHiVAG+vhj+OMf4fe/931gCgs3Pd6gATz1FAwYEEx80WzYAPfeCytXpv7eJSVVieKMGfDggzB+PHz3nYbSi4jkMs0snYH23BPuvBP+8pfITT6VKip8Z+RWrXx/nKCNGQPnnecXUj3mmKCjgT/8AT7/vKqf0tq1vpkxnSVqigEREfFUI5QhNmyAyy6DH37wo6DOO6/mJAj8eR9/DJdeCl9/nZo4o1mxAq66CvbeG44+OthYqtthB/85TZ/u5xx6882gIxIRkVRSIpQBVq3yTWG33OInBoyVme8b07IlnHyyH64elNtv93P73H57/LNIp8IWW0B+vm9OVFOZiEjuUCKU5pYs8UtH/PvfcM89cOGF8V3fsSM8/DBMnQo33picGGvzyy8+iTv2WF8jlI623tovv7HffnDiiXDHHUFHtDktEyIiknhKhNLYL7/A/vv7L79x4+rez+eYY3xn4ZtvhvnzExtjLFau9AlGsiZPTJT8fHj9dV/7dvHF8Pjj8V3ft29VH55Ey9T15ERE0p0SoTTWqpWvqXj9db8URX3ce6+vSQhiSH2XLvDGG74/Trpr2tQnnXfc4ROidDF8uJ84szotEyIiUn9KhNLQ1KmwfDm0aOETiAMPrH+Zbdr40Wbgh4ynyh13wOzZqbtfIjRs6Jsg8/Jg2TI4/3zf2TtIkRbVBS0TIiJSX0qEUiiWppNXXvGzRV98sX+d6I7FzzwDu+zia4eS7c03/ft45pnk3ysWEyfGP+x80iS47z6fjC5YkIyoops50yfEAFtuGfmcdF5PTkQkEygRSiOjR8MRR0C3bnDDDcm5x4ABft6cU05J7vISGzb4YfvbbQdnn528+yTb4Yf7Wby/+AL23Rd++in59/z4Y/jTn/x/p1Gj/L677958uoSmTdN3PTkRkUyhRCgNOOc7Ep96qq95mDAB2rdPzr1atIDHHvNNKhddlJx7gO9oPG2a76AdviBsphk0CN5+2y9dss8+PilKNOfgX//yyVbv3r5J9JJLqiaeDF8mpHNnnySVlPjJNV9/PfExiYjkAiVCaWDxYt+Z+YQT4OWX/Zw2ybTPPvDXv/ov0pdfTnz5K1f6Trx77eWHzGeDfff1zWTbbQetWyeu3Mq5ncz8MzBvnp8mYfZs+PvfoaCg6tzqy4TMnOlfr13rZ+w+9FA466xgli8REclkSoQCtHatXwZjyy19c8gTT6RuiYdrrvFz+ixZkviy16zxX8zpOnliXe26K3zwgU9ONmzwv9fVrFm+6XDbbf00CeBHq02f7qdJaNkytnKaNPHPzkUXwUMP+fXnPv647nGJiOQaJUIBWbbMr3V16aX+dUGBXyw1VZo29Qu2nnRS4stu29bXNu27b3zX1aUzc6pVJnb33ec7td9/f3zXT5kCxx/v+//ceSccdFDV3EAdO/oRa/Fq1swnne++68vq2xd+/TX+ckREcpEWXU2RylmB16ypSnrmzoUhQ4KLqUED3zflscd836FENGPddZefPLG4uP5lpbNhw3zicc45vjnruutqv2buXN+0tcUWvgbnnHPiH/VVU6LYt6/vlzVpUlUfs3nzfIIVKy3qKiK5RjVCKRA+K/Ds2b5p5MILk1MjE4+KCt8Jd9gwmDOnfmV9+aUfLv/kk4mJLZ01b+47N//lL36E34EHbr78xfLlvr/PWWf5a7bZxi+V8vPPfsmRZAx9z8/3I93A9//abjsfQ0VF4u8lIpINlAilQKRZgSE95tdp2NDXCK1dC0OH+hqiuvrrX31tx1VXJS6+dNaokU8ijzzS16BUX/7ilFP8LN7nn+9raVav9scOPzz5neErFRfDwQf7GA45xCdgIiKyKSVCKRBt9t90mRV4xx3h1lv9kO2HH65bGW+9Ba+9BldeGX3yv2xkBp9+uvn+9ev9sY8+8n2xmjVLfWwdO8JLL/lkbfJk2G03ePrp1MchIpLOlAilQLQmkHSaFfjMM6FfP9+0FW9H2w0b/Jw3RUV1Xxg2k0VLaFet8n2CgmQGp50Gn3/uR72tWhVsPCIi6UadpVOgtNT3wanePJaXl16zApvBP/7hazfincxxwwa/QOmuu2b+5Il10blz5LXA0inR3WEHeP/9qpGJY8f6mrs//CHYuEREgqYaoRQInxW4sNC/LikJNq5wnTrBwIH+98WLY7+uSRPfL+joo5MTV7orLd18+Yt0S3TB9wcz8x2n771XkzCKiIASoZSpPivwjBnplwRV98ILPlmLZSmJkSP96KlclimJbqUGDXzn7osv3nQSxsopHqqPfBMRyXZKhGQz++3nazROOsmPJotm3jw/H84//5m62NJVJiW64Dtv33Zb1SSMe+/t+xJVH/k2bJiSIRHJfkqEZDNbbQWPPOI72F57bfTzrrnGf3HefHPKQpMEq5yEsW3bzTtSl5f7qR9ERLKZEiGJaNAgP+v1zTf75pJwX3/tk6Wzz/bD7yVz5edH7xOWLlM8iIgkixIhiequu3xfkbIy/7pv36olGC69NLcmT8x2mTDFg4hIMmj4vETVqpVfNqN58033OwcnnADHHJNbkydms0yY4kFEJBmUCKVQJi5kWZkEXX45fPghrFvn168qLYVTTw02Nkmcys7dQ4f6fl8NG/oRZene6VtEpL7UNCa1euIJv0jounX+tUYUZafKkW+77uonyWzSJOiIRESST4mQ1OqqqzZfjFUjirLXllv6Zs/WrYOOREQk+dQ0JrVK90VjJbHM4Lnngo5CRCQ1VCMktdKIoty0ZAk8+2zQUYiIJJcSIalVpqylJYl1331w3HHw7bdBRyIikjxKhKRWmbaWliTGGWf4pTjuuCPoSEREkkd9hCQmJSV+JmnIzGkAUiHbPpettoLBg2HMGLj+emjfPuiIREQSTzVCIrLRxImbJnQXXujnFbr//qAiEhFJLiVCIhLVzjv7dee++y7oSEREkkNNY1kk2U0z2db0I7F55pmq/mEiItlGNUIiUqPKJOjnn/2M0yIi2USJkIjU6uOP/Rpz48cHHYmISGJlfSJkZgPM7F4zu9nMzgg6HpFM1KuXn0DzttuCjkREJLGyOhEys5bAncCFzrnLgcFmtkPAYYlknEaN/AiyDz6Ajz4KOhoRkcRJ+0TIzJqEanPWm1lRhOODzGyKmb1vZh+YWXG1w72Bn5xz60OvpwB/SEHYIllnyBC/EOvttwcdiYhI4qR1IhRKfN4DtgYaRjjeC/gnMNg5tz9wE/CGmXUMnbIVsLzaJb+F9olInFq2hDPPhJdfhoULg45GRCQx0joRAloCJwGjoxz/G/CGc+5rAOfcy8B84OzQ8QXAFtXObxXaJyJ1cNFFfk6hdu2CjkREJDHiToTMrKGZtTOzbcysdRJi2sg596Vz7vsaTjkYKAvbNwXoF/p9MrCdmVXOl7QH8HpioxTJHe3aQUGB/72iIthYREQSIaZEyMwOMrNRZjYdWIWvdfkZWGRmK8zsYzO7PpUdkc2sLZAP/BJ2aB6wPYBzbgVwEXCPmd0KPOac+6GGMoeZWZmZlS1YoIojkUjWrIFDDoGbbgo6EhGR+qtxZulQYjMG2Bffv+Zb4GtgJbAeaAK0AToBfwUuN7MHgEucc+uSFzYALUI/14TtXwPkVb4INZe9HEuBzrmRwEiA4uJil4AYRbJO06bQsCHcey9cfLFfoV5EJFNFTYTMrCvwPr6zcm9ginMuanJgZnnAEcAVwGtm9odqo7WSYWXoZ/jk/02B8iTeVyTnXXIJHHwwjB0LQ4cGHY2ISN1FbBozs4bAU8DZzrk/Ouc+qSkJAnDOlTvnngJ64JvOrk14tJvebzGwFOgYdqgjELX5S0Tq78ADoWdPP5RefYVEJJNF6yPUGDjaOfdMvAU65yrwI72eqE9gMXobKA7bVxzaLyJJYuZrhb75Bl57LehoRETqLmLTmHNuNfBjXQsNJUP/q+v1cbgZmGhmXZ1z35jZYfg5h+5Pwb1Fctpxx8GiRbDvvkFHIiJSdzV2lq7OzFoBPYFy51xZaN9WwDlAa+Bp59yHiQzOzJoAb4bKBxhnZnOdc0cDOOemmlkJ8LiZrcJPutjfOTcvkXGIyOYaN4bzzgs6ChGR+ok5EQLOxM/cPB44wsya4TtT7wQYcJaZ9XPOTUxUcM65tUDfWs55CXgpUfcUkfg89xxMmwbXXRd0JCIi8YtnQsUjgNOdc0eEXp8A7AyMwM/efBUwPLHhiUi6mzwZbrwRZs4MOhIRkfjFkwh1ds49Uu11CTAbuNE5txK4FeiayOCCYGYDzWzksmXLgg5FJCOcf77vPH333UFHIiISvzqtNWZm7YH9gScrh9U75zbgm8gymnNuvHNuWH5+ftChiGSEggL405/gkUdg6dKgoxERiU88idAsMzs09PvloWv/WXnQzLoBqkYRyUEXXwwrVvhkSEQkk8TTWfpW4CUzWwpsCbzknPsSwMwuAc4CNKOISA76/e/hjDNg++2DjkREJD4xJ0LOuRfM7BBgIPArcC+AmRmwK34pjoeTEaSIpL8HHww6AhGR+MVTI4RzbgIwIWyfA4YkMigRyUy//eaH0w8Z4jtQi4iku2hrjTUIzRNUZ2bWsj7Xi0jmef55vwjr21rkRkQyRLTO0k2B/5pZnSbPN7MbgNvqHJWIZKTjj4eOHf1irCIimSBiIuScWwWcju8cXWpm7WIpzMx2N7Px+MkXL0lcmCKSCZo29ctuvPGGn21aRCTdWWgaoMgHzfbDryLfCfgCv5DqPGAlsB6/Sn2b0PEeQAHwBnCic25xUiNPsuLiYldWVhZ0GCIZZ/FiP7fQscfCmDFBRyMi4pnZVOdccfj+GjtLO+f+E5of6BTgz8DRQJMIpy4G3gZOc869Vf9wg2NmA4GBXbp0CToUkYzUtq3vJ/Ttt1BRAQ3qNG2riEhq1FgjtNnJZg2BIvw8Qk2AcmBuNq72rhohkbpbvx4axTUmVUQkuepUIxQutIzGD6FNRCSiyiRo7lxo1QpaagypiKQpVVqLSFLMmAFFRfDoo0FHIiISnRIhEUmKoiLYay+46y7fVCYiko6UCIlI0lxyCcycCf/6V9CRiIhEpkRIRJJm4EDYcUe49VaIY1yGiEjKKBESkaRp0AAuvhg+/RS++SboaERENpfQRKi+65OJSPY5+WT4/nvo1g369vWbiEi6SHSN0I8JLk9EMlzz5rDddjB2LEyeDO+95ztSjx0bdGQiIjXMI2RmI+pQnmYLEZHNjB0LgwfDhg3+9cyZMGyY/72kJLi4RESizixtZhURdjvAouz3vzjXMDGhBaPaEhunfffdd0GHI5IViop88hOusNDPNyQikmzRZpauqWlsIbBdta0Y+Bi4CtgP6BLa9gNG4Bdk3T+xYaeec268c25Yfn5+0KGIZI1Zs+LbLyKSKjUtsfG0c27jv+HM7Drgb86598LO+xH40Mz+A1wIfJD4MEUkk3XuHLlGqHPn1MciIlJd1Boh59y5Ybv6REiCqp//HrBnogITkexRWgp5eZvuy8vz+0VEghTPqLG2ZrZttINm1gloVf+QRCTblJTAyJHQtKl/3amTf62O0iIStHgSofeBV83scDNrXLnTzJqEOhi/HDpHRGQzJSXQuzf06QM//6wkSETSQ019hMKdD0wC/g1gZktD+1uHfs4DjkpUYCKSfSZO9D9ffx2eew4eeQQs0jhUEZEUiblGyDn3A7AbcAfwLZAX2r4FbgO6O+d+SkaQIpJdZsyAUaNg0qSgIxGRXBd1HqFcV1xc7MrKyoIOQyQrrVrlR4z17g3jxwcdjYjkgrrMI1RTYR3NrFfody3cKiJxad4czjkHXn5Zi7GKSLDiSmLM7GAz+wyYA7wb2n2AmX1uZv0THZyIZK+zz4ZmzeD224OORERyWcyJkJntB7wKbAm8BqwPHZoCjAXGmVnfBMcnIlmqXTsYMQL23TfoSEQkl8XcR8jM3gCmAxc659ab2Vzn3DbVjg8IHTs4OaGmlvoIiYiIZI9E9BHqBVzqnFsf6aBz7hVgxzrGlzbMbKCZjVy2bFnQoYjkhJUr4YEHYMWKoCMRkVwUTyJkwJqoB80a4YfTZzQtuiqSWl984fsLjRoVdCQikoviSYR+BGqaC/Ys4Lv6hSMiuaZ3b99P6M47YX3E+mYRkeSJJxG6ExhjZo+b2UlAk1Az0plm9lbouMZ/iEjcLr3Ur07/3HNBRyIiuSauCRXN7HLgenwCZYAL/VwP/M05lzWJkDpLi6RORQV07QotW0JZmZbdEJHEi9ZZOp61xnDO3WxmTwHHAF1Cu6cDzzvnZtU/TBHJRQ0awMUXw+OPw9Kl0KZN0BGJSK7QEhtRqEZIJLUqKnxCJCKSDPWuETKzVkBPoNw5VxbatxVwDn4F+qedcx8mJFoRyTmVSdD8+bBuHXTqFGw8IpIb4vn315nAROAqADNrBrwfen0u8J5mlhaR+li9Grp1g6uuCjoSEckV8SRCRwCnO+eOCL0+AdgZnwhtEfo5PLHhiUguadYMTjgBxo6FOXOCjkZEckE8iVBn59wj1V6XALOBG51zK4Fbga6JDE5Ecs+FF8KGDXDPPUFHIiK5oE5dE82sPbA/8KQL9bZ2zm3AD6UXEamz7beHY46Bhx6C334LOhoRyXbxJEKzzOzQ0O+Xh679Z+VBM+sGaIEuEam3Sy+F5cvh3XeDjkREsl088wjdCrxkZkuBLYGXnHNfApjZJfglNl5LeIQiknP22ANmzIDOnYOORESyXcyJkHPuBTM7BBgI/ArcC2BmBuwKvAc8nIwgRST3VCZBq1f7TtQiIskQ78zSE4AJYfscMCSRQQXJzAYCA7t06VLruSKSXBdfDBMnatkNEUmeuDtLm1lbM/uzmf0ttP3ZzNomI7ggOOfGO+eG5efnBx2KSM7bdVf473/hnXeCjkREslVciVBo0dXZwFigNLSNBWab2WWJD09EcllJCXTsCLfdFnQkIpKt4lli43zgRuBT4FVgXuhQR+Bw4EYzW+Wc0+wfIpIQTZvCeefBFVfAtGnQvXvQEYlItol50VUz+xG40zl3b5Tj5wIXOOd2SGB8gdGiqyLpYckSKCiAY4+F0aODjkZEMlW9F10FWkVLgkLuA0bEHZmISA3atIFnn4VevYKORESyUTx9hOaZWV4Nx1sAm6wOZGaaBURE6u3QQ6F9+6CjEJFsFE8idC/woJk1Dz8Q2nc3cH3Yocn1iE1EZKOPP4b+/WGZ5q8XkQSKp2msN3AIMMfMPgHmh/Z3APbAT7JoZjag2jUagy4iCdGoEbz5JjzyCFxySdDRiEi2iKezdEUdynfOuYZ1uC5w6iwtkn4OPBCmT4cff4QmTYKORkQySbTO0vE0jc13zjWIZ8PXEomIJMQll8CcOfD000FHIiLZIp5E6IE6lF+Xa0REIjr0UD/b9K23QoyV2SIiNYpn0dXwjtCbMbNDnHNvxnONiEiszOCGG3yt0IYNvt+QiEh9JPrPyBhgmwSXKSKy0ZFHBh2BiGSTuBIhM+sNXATsCmw2jB7YKhFBiYjUZPVqeOwx2Gcf2G23oKMRkUwWcx8hMzsAmATsCywEtgVmhrYKoAj4JPEhiohsavVq33H65puDjkREMl08naVHAPcAnZxzfYDFzrkDQlsX4Brg3STEKCKyidat4bTT/OixmTODjkZEMlk8iVAP4ApXNfFQ+JiNUmAAIiIpcMEF/ufdd1ft69vXbyIisYonEVrlnFtT7fUGM2ta7XUFkPFri5nZQDMbuUzz+Iuktc6d4c9/9jNNL10adDQikqniSYQWmNn/VXv9E3B6tddnAL8lJKoAOefGO+eG5edrdRCRdHfJJdCtG/zyS9CRiEimimfU2CvAK2Z2u3PuWmA0MMrMhuGbyboBf09CjCIiEfXs6RdjFRGpq3gSobuAD4Gloddj8MPohwJr8LNIX5e40EREYrNoEfz8c9BRiEgmimdm6QX4WqHK1w64JLSJiARm4ED47TfYcks/+7SISKzi6SNUKzP7WyLLExGJxZlnwldfwUcfwXvvQVERjB0bdFQikgkSmggB5ya4PBGRWlVU+Jqgdev865kzYdgwJUMiUruoiZCZbYh3AzqkMHYREQCuvnrz1ejLy2H48GDiEZHMUVMfoTXA03GUZcCx9QtHRCR+s2bFt19EpFJNidAy59yQeAozs/71jEdEJG6dO0deaqNzxk/xKiLJVlMfoT51KK8u14iI1EtpKeTlbbovL8/vFxGpSdREyDk3PdJ+M2tsZr3N7KjQ6xZm1qima0REkqmkBEaOhKbVFv257DK/X0SkJnGNGjOzU4DZwAf4CRUBegNzzWxoQiMTEYlDSQn07g377uvnE5o8OeiIRCQTxJwImdkg4B/ADOAhIDRQlfeBc4DrzezoRAcoIhKPRo38GmSvvaZkSERqF0+N0F+Bvznn9nLOnQ2sBXDOrXPOPQMcg2aZFpE0cM450K4dXHNN0JGISLqLJxHqCtwe7aBz7iOgU70jEhGpp5Yt4a9/hYkTI48mExGpFE8iVOMKPmbWDGhev3BERBLjrLPghx+gsDDoSEQkncWTCH0BXFzD8WuBz+sXjohIYrRoAdtu639fsSLYWEQkfcW8+jxQCrxmZgOAt4DmZnYuvjnsSGB7oF/CIxQRidHEiZvv+/OfYeFCePvtlIcjIhkg5hoh59ybwEnA7/C1P/nAXcClQDvgBOfcxMSHKCJSd3vuCe+8A5MmBR2JiKQjc+ErFdZ2gVkLoD/QJbRrOvCWc25lgmMLVHFxsSsrKws6DBGpp/Jy2GEH6NoV3n036GhEJChmNtU5Vxy+P56mMQBCCc/zCYlKRCTJ8vLg8svhggt801nfvgEHJCJpJa6ZpUVEMtGwYbD11nDXXUFHIiLpJu4aIRGRTNO8Obz8Muy0U9CRiEi6UY1QGDMbaGYjly1bFnQoIpJAu+/uJ1rcsAHi7BopIllMiVAY59x459yw/Pz8oEMRkQT79lvfafqdd4KORETShRIhEckZhYWwahVcfbVqhUTEUyIkIjmjWTO44gr48EN4662goxGRdBBXImRmHczsLjP7n5nND+3b08zuNrNtkxOiiEjinHoqFBSoVkhEvJgTITPrDHwGnIdfVqNp6NB8YEdgipntkOgARUQSqWlTGD4cJk+GCROCjkZEghZPjdANwHfArs65lkA5gHNupnPuMOA+/NIbIiJpbcgQeOYZ6NMn6EhEJGjxzCN0MNDLOfdLlOO3At/WPyQRkeRq0gSOPTboKEQkHcRTI9SkhiQI59w6oEX9QxIRSY1HH4VDD1VfIZFcFk8itMTM9oh20MwOAhbXPyQRkdQwg9dfh1deCToSEQlKPInQk8C/zexUMysAMLMtzKyrmQ0HngUeS0aQIiLJcPLJsP32cM01qhUSyVXxJEI34keNPQrMADoAS4EvgeuBSfh+QiIiGaFxY7jySpg6FcaPDzoaEQmCuTj+GWRmBpQAfwIqh8pPB552zj2V+PCCU1xc7MrKyoIOQ0SSbP162GUX2GIL+O9/fXOZiGQfM5vqnCsO3x/X6vPOZ01PhjYRkYzXqBE88ADk5SkJEslFcSVCtTGzx5xzgxNZpohIsh1ySNARiEhQ4k6EzKwjfmbpvAiHB9Q7IhGRAKxYAZddBv36wZFHBh2NiKRKzImQmW0HPAHsHe0UQOMuRCQjNW8O77wD778PgwZBAy1JLZIT4qkRehToiF9qYzawNuy4AXcmKC4RkZRq2BBGjICSEvjXvzTztEiuiHnUmJktBXZyzv1awzkvOueOTExowdKoMZHcs2ED7Labrw2aNk21QiLZJNqosXj+N/8OWFTLOUfFFZWISBqprBX66it47rmgoxGRVIgnEboC+JuZ1XTNnHrGIyISqGOPhcsvh169Nj/Wt6/fRCR7xNxHyDn3lpltDcw0s8+BhUBF2Gn5iQxORCTVGjaEm24KOgoRSZV4Ro0NBv6B7xS9bZTTNGpMRLLCtGnw0ENw770+ORKR7BRP09gI4B5gO6Cpc65B+AZE7UgtIpJJpk+HBx+EceOCjkREkimeRKiVc+5C59xM59y6KOdclIigRESCdvTR0L07XHedX49MRLJTPInQJ6FZpWuipjERyQoNGsDVV/uaoaeyaklpEakunkTodOBWMzvUzNpFOeeOBMQkIpIWjjwSevRQrZBINotnZumZoZ8nAJiWaRaRLNegAdx4I3z0EayL1iFARDJaPInQGuDpGo4boEnpRSSrHHaY30QkO8WTCC1zzg2p6QQz61/PeERE0o5zcMkl8OGHvmaoqAhKS/26ZCKS2eJJhPrEcE63ugYiIpKu/vlPuPNOnxABzJwJw4b535UMiWS2mDtLO+emx3BaaT1iERFJS8OHVyVBlcrL/X4RyWzxzCzdOYbTjgHOqns4IiLpZ9as+PaLSOaImgiZ2QSgLbC7c24DMAPNEyQiOahzZ98cFmm/iGS2mmqEjE2bzlYCt9VyfsbPLG1mA4GBXbp0CToUEUkTpaW+T1B5edW+vDy/X0Qym7nwhu9oJ5r94pzbur7nZIri4mJXVlYWdBgikibGjoWhQ2HNGigs9EnQfvv5WiFNqyaS/sxsqnOuOHx/PDNL71ZD4Xm1nSMikslKSqB3b+jTB2bMgG7dYJdd4OGHg45MROojnkTouhqO/d3MFgB71jMeEZGM0L07HHAAnHeen3laRDJTPInQkTUcuw7fP+j2ekUjIpIhGjaEJ5+EggL44x9h3rygIxKRuognEYrKObcAeBI/ykxEJCe0bQvPPw9LlsBxx2k9MpFMVOM8Qmb2j2ov881sFH50WKRydgbmJjA2EZG016MHPPooaGyFSGaqbULFU6r97oBoa42VA9+gyRRFJAedcILfADZs8M1mIpIZamwac841qNyAX6u/DttaOuf2cM6py6CI5Kxp02DXXf1PEckM8fQR+lvSohARyQLt28Py5XDUUb7fkIikv3gWXR2TxDhERDJex47w3HPw889+3qGKiqAjEpHaJGTUmIiIeHvvDffeC6+9BtdeG3Q0IlKbmFefFxHJdRMnxnbesGHwySfw8cfqPC2S7pQIiYgkmBk88AA0aqQkSCTdqWlMRCQJmjb1SdDcuXDqqbBiRdARiUgkSoRERJLof/+Dxx7zyZBzQUcjIuGUCImIJNGBB8JNN8Gzz8LtWo1RJO0oERIRSbJLL4VjjoHLLoN33w06GhGpTomQiEiSmcHo0bDzznDllWoiE0knGjUmIpICW2wBr7wCbdr4xEhE0oNqhEREUmS77aB1a1i9Gp54YvPjffv6TURSR4mQiEiK/eMfcPLJ8MgjQUciIkqERERS7PTToX9/OOccP/u0iARHiZCISIo1bAj//Cdsu60fTTZ/ftARieQuJUIiIgFo2xaefx4WLYLTTgs6GpHcpVFjIiIB6dkTxo3zw+rHjoXJk2HNGigqgtJSKCkJOkKR7KcaIRGRAB1xBEyd6lesX7PG75s5078eOzbY2ERygRIhEZGADR8O5eWb7isv9/tFJLmUCImIBGzWrPj2i0jiKBESEQlY587x7ReRxFEiJCISsNJSyMvbdF+zZn6/iCSXRo2JiASscnTY0KG+w3TjxnD++Ro1JpIKSoRERNJASUnVkhsTJmhhVpFUUdOYiEiaMYP16+GWW+Czz4KORiS7KRESEUlDv/0Gd9wBp5wCa9cGHY1I9lIiJCKShtq2hYcfhs8/hxtvDDoakeylREhEJE0dcQSceKIfPaYmMpHkMOdc0DGkpeLiYldWVhZ0GCKS4xYvhl13ha23hrIyaKB/vorUiZlNdc4Vh+/XqDERkTTWti08+aSfZ0hJkEjiKRESEUlzBx1U9fuqVdC8eXCxiGQb/ftCRCRDXHcd9O5dtUq9iNSfEiERkQyx++4wbRpcf33QkYhkDyVCIiIZ4vDDYfBguPlm33FaROpPiZCISAa56y7o0MFPtKgmMpH6UyIkIpJBWrf2a5L99BNMnRp0NCKZT6PGREQyzGGH+USoffugIxHJfKoREhHJQO3bg3PwwguwenXQ0YhkLiVCIiIZaupUOPpouPbaoCMRyVxKhEREMlRxMQwdCrfcAp98EnQ0IplJiZCISAa7/XbYZhs/ikxNZCLxUyIkIpLB8vPh0Ufhm2/g6quDjkYk8+REImRm25vZC2Y2LuhYREQSrX9/uPxyv/xGkPr29ZtIJsmV4fN7Aa8DBwQdiIhIMtx0U9ARiGSmwGuEzKyJmd1sZuvNrCjC8UFmNsXM3jezD8ysON57OOeeAjQHq4hkNed8x+krrww6EpHMEWiNUCjxeQqYDjSMcLwX8E9gT+fc12Z2OPCGme3qnJsXOqeMyO/jMOfc3KQFLyKSZszghx/8zNMDBsDeewcdkUj6C7pGqCVwEjA6yvG/AW84574GcM69DMwHzq48wTlX7JzrGWFTEiQiOefWW6GgwI8iW7UqdfcdOxYmT4b33oOiIv9aJBMEmgg55750zn1fwykHA+FrLE8B+iUvKhGRzNWqFYwaBdOnw1VXpeaeY8fCsGFVi8DOnOlfKxmSTBB0jVBUZtYWyAd+CTs0D9g+zrIGAAOBbmZ2Xg3nDTOzMjMrW7BgQbwhi4ikhYMPhtNPh7vvhlmzknuv336DCy+E8vJN95eXw/Dhyb23SCKkbSIEtAj9DO/kvAbIi6cg59wrzrljnHPdnXP31HDeyFBTW/FWW20VZ7giIunj1lvhP/+Bzp2TU/6IEbDXXtC2LUT7d2OykzCRREjnRGhl6GfTsP1NgbB/e4iISHVbbOETFYAff6x7OatWwTvv+JFo555btf/996FRIz9/Ufv2ka9t2hS+/bbu9xZJhbSdR8g5t9jMlgIdww51BH5IfUQiIpnnmWfghBNg4kTYb7/Yrxs7FkaO9B2g166Fhg1h3339EH0zePddaBD6p3TXrr5PUPXmscaN/c/ddoMLLvCJVKtWiXpXIomTzjVCAG8D4fMGFYf2i4hILQ47zI8iO+YYaNbMJzHVR3WtXQsffAA33OD7Fi1e7PfPnQsrV8J558Err/j9773nr4eqJAigpMQnTU1D9feFhTB6NMyYASed5JvpLr00Ve9YJD7mnAs6BsysLzAB2M45N6Pa/l7ARPw8Qt+Y2WHAk0C3ynmEkqW4uNiVlYUPWBMRyTzDh8ONN266r1kz6NLFN5tV1uT06AFPPgm/+11VzU88KpfXmDhx0/2ffOIXhu3UyTeVrVgBvXrV5Z2I1J2ZTXXObTYpc9ATKjYB3gRah3aNM7O5zrmjAZxzU82sBHjczFbhJ13sn+wkSEQkm0Qaxr56tU9KTj8dDjgA+vSBLbesOh5vElSTPfes+n3ECHj2WfjLX6C0FDQuRYIWaCLknFsL9K3lnJeAl1ISkIhIFoo2emv9erj33tTGMnKkrxm65x6fEF13HZx5pu94LfUTrUZOapbufYRERKSeog2hT/TQ+okTa/8Szs+H22+HadNgjz18H6S77kpsHCLxUCIUxswGmtnIZcuWBR2KiEhClJZCXtjsa3l5fn9QunaFN96Af/8bzjjD75s6tea5h/r2rar1EEkUJUJhnHPjnXPD8vPzgw5FRCQhKkd1FRb6vj+Fhf51SUmwcZnBoEHQsqXvnH3aabDLLnD99b4Pk8ROa73VnRIhEZEcUFLih7NXVPifQSdB4czghRdgwADfobpbN3jxRZ8gSc201lv9KBESEZG0UFjoO1C//TY0bw5HHQXjx/tjqvGIbvjw7FjrLaimTyVCIiKSVg46CD77DB57zNcQjR0LQ4eqxiOS8nL/eUQycybcdptPLBcurP+9kpmoBJnoKhESEZG007gxnHyyX9rjb3+rSoIqZWKNRyItWuRn6+7UKfo5DRv6c/r18/M1bbcdbNjgj02bBt9/75tKgxZ0055mbhARkbQ2e3bk/ZU1IWvX+sVhs32Mi3Pw66/QoYNf4mTkSOjf33cwv/32TZvH8vKqjn/+ua9hW7jQJ0fg13+bMMF3VO/RA3r29GvR/fnPke9dWWOzZo2vsSktjb+fWeVs5b/9BmVlPplbuNAnutGa9lLRly0tlthIR1piQ0QkPRQVRW7+6dgRfvnFLwDbrx/svrufIbtvX/+l3rp1igNNkuXL4fHH4b77oEULmDLFJxTLl8MWW/hzqjcfFhbWnqh8/rlPRiqTpM8+g332gddf98cHDvRl9+zp15m75x6fbFbKy4O//x32398nMwsXwq67+m3OHLjssqr9lQnPPffAkCHw8cfQu3ft79sssTVWabnEhoiISG1KSzdf3T4vz/d/Af/Ff+WVvn/Jfff52hEz/+XevTvMm+cXhG3TJpDw6+yHH3zyMGaMr0UpLoZzz62qWalMgsAnPY884n+PZWbpHj38Vqmiwt+j8vdGjfxivE89Ffn68nIfS3U33OATIYCPPvJLtrRv70cAtmvn544C/3rCBL+vXTvYa6/I80clesLPaJQIhTGzgcDALl26BB2KiIhQVbMRrcZjhx3g2mv976tX+xqH99+v+uK9+WafUPToUdXh9//+D9q2rbpHvDUqyVJR4Zc+adLEJzQPPgjHHeeTjr32St59GzSoqkFr0MBPZQC+Nqhdu+jTGDz3nE942rWr6q+07bY+iYtmiy027XR9442RE91UTfipprEo1DQmIpJe6rqW1tSp8Morvsboww99srT11r4Jx8x/4ZaWbt70k8pJJ5cuhdGj4f77fdJz/vk+nmXLfBNgLJK11li0psnCQj8nVSKkIhGN1jSmRCgKJUIiIuklEV/0a9b4Pjbz58Mxx/h9jRpVjaaqrvKLft48P+qqsqNxXUWK/6uvfHPe44/7GpF99vHNfIceWr97JVLlqK5InbHTbWLOmqiPkIiI5LymTX1H6krORe+QO2uWT5AKCvx5W2/tm386dfLNVcce66+dPNnv23prP+w/kuqjrgoLfXNQSQmcc47vT3P88b4maPfdE/+e66u2pslMp0RIREQyQqKbfMA3jXXuHLnpp3Nnnwjdd58fwl+5ffkl7LmnP2fePNh336qyOnTwSdGll/pkadkyuOIKGDWqap6cWbP8umrg+wBVdhpOZ/F2xs4kSoRERCSnRRuVVlrqOy2ffnr0a1u3htde2zRRmj3b1zwBfPstPPDA5tetWuXnyUlUH5tUyLYEqJISIRERyWmVTTzDh/vams6dY2/6ycuDP/wh+vGePX1NUaTuuJGGjEvqKRESEZGcV1KSnD4vTZrU3PQmwdNaYyIiIklUWuprjqpL5Tw5UjMlQiIiIklUUuKHmhcW+maywsLMG3qezdQ0JiIikmTJanqT+lONkIiIiOQsJUJhzGygmY1ctmxZ0KGIiIhIkikRCuOcG++cG5afnx90KCIiIpJkSoREREQkZykREhERkZylREhERERylhIhERERyVlKhERERCRnKRESERGRnKVESERERHKWEiERERHJWUqEREREJGcpERIREZGcZc65oGNIS2a2DPiultPygZoWJavpeDtgYR1CC1pt7zkd71OfsuK9NtbzYzkv156vVD1bib5XNj5f2fZsQe49X+n6bNV2PJnPV6FzbqvN9jrntEXYgJH1Paem40BZ0O8xWZ9Lut2nPmXFe22s5+v5Su5/81TeKxufr2x7thL93zyV96prWen6bNV2PIjnS01j0Y1PwDmxlJFpUvWeEnmf+pQV77Wxnq/na3OpfD96vmo+J9ueLci95ytdn6147pUSahoLiJmVOeeKg45DspOeL0kWPVuSTEE8X6oRCs7IoAOQrKbnS5JFz5YkU8qfL9UIiYiISM5SjZCIiIjkLCVCdWBmg8xsipm9b2YfmFmN7Zlm1srMxoSu+a+Z/d3MGoWds7WZ/dvMPgqdc0mUsk40syVmdkoC35KkiWQ8W6Hz+pvZbDO7Jko5u5jZu2Y2ycymmtlJCXpLkkbifb5C1+xhZl+b2Zgox2v926XnKzeYWRMzu9nM1ptZUQznJ+S7sd7PV6qHqWX6BvQCVgDdQq8PBxYBHWu45nngidDvTYAPgRurHW8AlAHXh17nAz8Aw6qd0xh4DngQcMApQX8W2tL/2Qrtvzf07MwHrolQRktgFjA09LoTsBjoH/Rnoi3w5+sy4G3gK2BMhOOx/O3S85UDG1AEfAQ8FvqOKorhmkR8N9b7+Qr8w8u0LfSF8q+wfV9X/oeKcP7vQg/FbtX2HQeUAy1Drw8H1gOtq53z19B/3Mp+XHnAIaHflQhl4ZaMZyu076jQzxlREqFzgAVAg2r7HgAmBf2ZaAvu+QodHwQYMDFKIhTL3y49Xzmwhf4edQH6xpIIJfC7sd7Pl5rG4ncwPkOtbgrQr4bzVwNfhp3fHNiv2jk/OOeWhp1TAOwM4Jwrd869Wa/IJd0l49nCOfdCDPf9r3OuIqycfcwsL4a4JTPE+3zhnHvJhb5Zaiizxr9d6PnKCc65L51z38dxSUK+G0nA86VEKA5m1hZfNfdL2KF5wPZRLtsemB/2x2RetWOVPyOVWf0cyWJJfLZiEe35a4Cv7pYMV8fnKxax/O3S8yWRJOq7sd7PlxKh+LQI/VwTtn8Nvukq2jWRzqfaNbGcI9ktWc9WrPfW85fd6vJ8xVqu/r5JXSTq2an386VEKD4rQz+bhu1vim/XjHZNpPOpdk0s50h2S9azFeu99fxlt7o8X7GWq79vUheJenbq/XwpEYqDc24xsBToGHaoI74neyQ/Au3NzMLOp9o1P0Yps/o5ksWS+GzFItrzV4HvYC0Zro7PVyxi+dul50siSdR3Y72fLyVC8XsbCJ97ozi0P5K38J2/dg07fxXwQbVzuphZ67BzfnbOfVvfgCVjJOPZisVbwO5mVv3vQTHwoXNO/2LPHvE+X7GI5W+Xni+JJFHfjfV/voIecpdpG34ujuVA19Drw/BzFnQMvb4B3wu+WbVrngceC/3eGPgPkedKuDb0uhXwHdXmSgiLQcPns3BLxrMVVv4Map5HaEjo9bb4+WU0z0sWbXV5vqpdO5Ga5xGK+rdLz1dubUQZPp+s78ZEPF+Bf2iZuOHn1pgCvI/PXPeoduw2fFVd82r7WuEnmZoC/Be4BWgUVuY2wEv4Cak+BS6NcN8HQn+QHPC/0O+dg/48tKX9szUi9KysxidDE4E9w87ZBZgATAqVc3LQn4W2tHi+hoWel6X4kTgTgUFhZcbyt0vPV5Zv+AkRJwKfhb6jJgPP1/J8Jeq7sV7PlxZdFRERkZylPkIiIiKSs5QIiYiISM5SIiQiIiI5S4mQiIiI5CwlQiIiIpKzlAiJiIhIzlIiJCIiIjlLiZCIiIjkLCVCIiIBMDNnZmOCjkMk1ykREhERkZylREhERERylhIhkSxiZgeb2TwzW2tmM8xsVzN728zmV2+KMbOtzWxU6NzFZva9md1kZnkRyuxlZu+Z2XIzm2Vm/zCzo0PlLTOzr2qIZ0joHhvMbGLouk/N7Fczm2tmN5hZo7Brtjezu8zsm2rxvWpmv6/D+6y1rCjl7RZ6z0vM7CszGxA6r8TMvjSzpWb2gpltFeE91/jZVn4modP/FDpvnpk9FGc5MX0G8TLvVDP7wMwWmdlqM5tpZi+bWeO6lCmS1oJesVabNm2J3/CrQC/Er9rcCTDgSWAM0B6YiV99fNvQ+b2BufgVnBtUK2dn4DfgXWBL/D+ejg1d74BTYoxnBrAAeAPYKrTvj8B6YGTYueeEYtkt9Lo18BSwHNg+1vcZb1lh5T0B5AONgeeA1cBg/GrsDYCdgMXAP8Ouj+ezdZVxRogjnnJq/Azq8Ow8HIrtOeDM0Hu+Ify9atOWLVvgAWjTpi3xW+jL0QF7Vtu3G3AI8Gjo2E5h15wV2n90tX3Phvb1CDu3soxTYoxnRijpKQzb/3yonN2q7TsOOD3svJah62+L9X3GW1ZYeb+rtq9vaN/HYec+HkqQLMLnEstnW1MiFE85NX4GcT43+cAG4OGgn2Ft2lK1qWlMJHutds59UvnCOfcF8Da+JmaGc2562PlTQj//AGBmDYFDgXnOuc/Dzp1Qh3h+dM7NDNv3aujnwGpxPuOce7j6Sc65FfjakF0jlLvZ+3TOvVnHsgDKnXNfVns9N/Rzath5s4GmwFYAZtaAGD/bmtSxnKifQZzW4WsAe5nZnmbW3szy61COSMZoVPspIpKhFkTYtxX+X/151fqpVDJgZeicynNbAN9FKCf82ljMr6Gc7TYGYdYU3yRzPFBAVV/GrYAfI5QR6X3WtSyARWGv14Z+Loyyv7JfVTyfbU3qUk7UzyAezrlyMxsEPA18HNr9LL5mTSQrKRESyV4VNRz72jnXM1WBxGkscARwEvC8c24tgJnNiHJ+Te8z3rJqKq+m+1SXqM82nnJija1GZnYMMBp4EXgZ3w9qRiLKFklXSoREcssCYCmwbaSDodFUK0NNMgvwtQ8dI5waaV9tOtRQzk+h+7cBjgZed86Nq8M9NkpkWTGK57NNRTlxCX1ej+H7LZ2TyLJF0pn6CInkEOdcBb6po52Z/V/1Y2bWHN+HaLfQuRvwfXg6mlmPsKL61uH225tZ57B9A0I/x4d+rsV3/HVhsTUhciJVk0SWVat4PtuQckL/GA31xbnPzFrWoZxE2Q3fFJrQBEsk3SkREsk9w/HNHfeaWRcAM2uHHzL+JX4YdqWrgBXAHWa2pZk1CDWfFNfhvnOAe0L3qmyGGQg8EurIjXNuJfAacIiZHRk6Lw+4F2gWz80SWVYc4vlsvwF2DnWO7o/vx7SyDuXUysz+Y2bLajlteuj+N4bmXjrNzM4xszvN7Kl47ieSUYIetqZNm7bEbUAPfAfktfhh0POAJyKc1wF4CJ+c/ILvEP13ID/Cub2A9/Bz7/wE3IOvyXHAyTHGNQM/zPsAfCfc+fjRWDcAjcLObYufy2YO8CvwX+AiYFbofc2L433GUtauUcr7P+AyfFOVwyeE/w2V+9/Qaxc6fmO8ny2wD/B56PpvgYHx/DeK9TMInTuTsOH/Uc7bFz/X02L8FAMLgEnA4KCfbW3akrWZc5vUGouI1MrM/gSMAw53zr0Sw/kz8MPB+yY5NAljZnsBk/HzKY0MOh6RdKOmMRGJKrTUxCMRDu2Nr4mYnOKQJA6hJUBGAy8Akf47iuQ8JUIiUpMWwBAz+1NoDaqGof42fwFucc6Fz7kj6WU5vk/UsU7V/yIRqWlMRKIys7b4fjKHAe3wkwfOBB50zj0Yw/VDgJvwEwBuwPc9OdY5NylpQYuIxEGJkIiIiOQsNY2JiIhIzlIiJCIiIjlLiZCIiIjkLCVCIiIikrOUCImIiEjOUiIkIiIiOUuJkIiIiOQsJUIiIiKSs/4fgOJ44QG1SfIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (9, 6))\n",
    "\n",
    "plt.errorbar(regs, tm1, yerr = tm1_err, fmt = 'bo--', label = 'vanilla sinkhorn')\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.yticks(fontsize = 14)\n",
    "plt.xticks([1e-3, 1e-2, 1e-1, 1e0], fontsize = 14)\n",
    "plt.gca().get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())\n",
    "\n",
    "plt.xlabel(r\"reg parameter, $\\varepsilon$\", fontsize = 18)\n",
    "plt.ylabel(\"time elapsed (s)\", fontsize = 18)\n",
    "\n",
    "plt.legend(fontsize = 16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.7 ms ± 289 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "6.27 ms ± 20 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "112 ms ± 29.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "152 ms ± 8.41 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "233 ms ± 64.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "891 ms ± 33.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "5.31 s ± 333 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "ns = 2 ** np.arange(5, 12)\n",
    "\n",
    "tm1, tm1_err = [], []\n",
    "\n",
    "for n_samples in ns:\n",
    "    x_i = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "    x_f = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "\n",
    "    C = cdist(x_i, x_f)\n",
    "    C /= C.max()\n",
    "\n",
    "    a = np.ones(n_samples) / n_samples\n",
    "    b = np.ones(n_samples) / n_samples\n",
    "    \n",
    "    t = %timeit -o sinkhorn.sinkhorn_knopp(C, reg = 1e-2, a = a, b = b)\n",
    "    tm1.append(t.average)\n",
    "    tm1_err.append(t.stdev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9c487f8850>"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAF/CAYAAABKac3QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABIDklEQVR4nO3deZzW8/rH8ddVjbQoJZWlBYlQoeFkOSokJOE4lkIcKsLh2A5lOZwTjn1fwpEfOZZzLIXsJbsmS3UI2UqpRA1ap+b6/fH5Tk13M9PcM/c933vu+/18PO7HPd/Pd7vuu7u5r/ms5u6IiIiI5Ko6cQcgIiIiEiclQyIiIpLTlAyJiIhITlMyJCIiIjlNyZCIiIjkNCVDIiIiktPqxR1ApmrRooW3b98+7jBEREQkBaZMmbLQ3Tcva5+SoXK0b9+egoKCuMMQERGRFDCz78rbp2YyERERyWlKhkRERCSnKRkSERGRnKZkKIGZ9TOzUYWFhXGHIiIiIjVAyVACdx/n7kOaNm0adygiIiJSA5QMiYiISE7T0Ppq+OWXX1iwYAFFRUVxhyKSFnl5ebRs2ZImTZrEHYqISNooGaqiX375hfnz57PVVlvRoEEDzCzukERSyt1ZtmwZc+bMAVBCJCJZS81kVbRgwQK22morGjZsqERIspKZ0bBhQ7baaisWLFgQdzgiImmjZKiKioqKaNCgQdxhiKRdgwYN1BQsIllNyVCCZIbWq0ZIcoE+5yKS7ZQMJYhzaH3PnuEhIiIiNUfJkIiIiMRqzBjYeGMwg/btw3ZNUjIkGeOee+5hp512wswYPXo0AAUFBey6665stNFGnHzyybHGVxmXX345HTp0wMyYOHFiuceddtpptG3bFjPj22+/rbH4REQyzZgxMGQIrFgRtr/7LmzXZEKkZEhic8stt9C+fXuWLl0KwOmnn84LL7ywzjH5+fl8/PHHbLnlljUW15IlS2jXrh233npr0udeddVV3H///Rs87v777+eqq66qSngiIlllxAiIvgbWWLo0lNcUJUMSm+bNm9O2bVvq1q0bdyjrqFevHm3btmWzzTaLOxQRkaw3a1Zy5emgSRczxJgx8N57oZqwfXsYORIGDow7qvQ66aSTOOmkk+IOYz3169fnzTffjDsMEZGc0LIlzJ+/fnnbtjUXg2qGMkDc7aXLli2jS5cumBlbbLHFmgRl7ty57LrrruTl5dG1a1cmTJhAv3792H333enatSu/+93v1mvWOvTQQ2ndujVmxuTJk+nTpw/t27ene/fu/O9//1tz3MiRIyvVt6Y8y5cv5+KLL6Zbt27svvvudOnShdNPP53Fixdv8Nw5c+Zw7LHH0qVLF3bbbTe6d+/OtddeC8BPP/3ErrvuSuPGjelZamhfZV9XWcaMGUP79u1p2LAhu+66KxMmTFhn/1dffcXhhx/O9ttvT+fOnZk0adJ61xg7dix77LEH22+/Pe3ateOUU05ZZyLE0n2QXn/9dfr160enTp0wM5555hl23XVXmjdvTvv27Rk/fjz7778/W2+9Nb179+b777/f4HsmIpIOH3wAhYWh43RpDRuGSoEa4+56lHoA/YBRHTp08Ip8+umnZZb36LH+4847w74lS8rev9lm7rD+o379sP+xx8L5s2aVff7YsWH/jBkVhrxBnTp18r59+65T9tVXX3m3bt3c3X3o0KE+fPhwLy4udnf3t99+2xs0aOCTJ09e55wrrrjCAf/zn//sq1ev9qKiIt9vv/18r732Wue4CRMmOOATJkxYU/bNN9844A8++OA6x7Zr184HDRq0ZvuHH37wli1b+qxZs9zdfdmyZf7HP/7R+/Xrt8HXecABB/jgwYPXvI7nnnvOw3+FtXr06OE9evRIyeuaP3++d+vWzV9//fV1jnvwwQcd8KOPPtqXLl3q7u4DBgzwNm3aeFFR0ZrjHnvsMa9bt64/88wza15rnz59vFOnTv7bb7+td72+ffv6okWL3N398MMP96efftrd3QcNGuRNmjTxyy67zN3df/31V+/YsaMff/zxG3zPyvu8i4hUx7vvuu+2m/utt7q3a+duFp4feST19wIKvJzvftUMJfAY5hn66aeyy0tqimrKoEGDeOmll5g3b96astGjR6+pKRo+fDiXXnrpmkn49t57b7p06cIDDzxQ5vVOOeUU6tSpQ7169ejXrx/vvvsuK1euTEmsLVq04J133qFNmzYAbLzxxvzpT39i3LhxzC+rvrWU9957j/bt2695HX379mX48OGVvncyr+unn36ib9++XH311fTq1avMYwYOHLhmNvMjjzyS2bNn8/XXXwPhj5ULL7yQXr160b9//zWv9brrruOzzz7j3nvvXe96gwYNYtNNNwXgvvvu44ADDliz79dff+Xcc88FoHHjxvTu3btKNXMiItVR8r3XvTsUFMCf/wzffgvFxeG5pruJqM9QilX0vdKwYdn727cPTWOJ2rVb9/g2bSq+/g47VCrEcp1wwgkMHz6cMWPGcP755+PuPPbYY7zzzjsANGrUiEsvvZSJEydSVFREnTp1mDlzJuUljh07dlzzc/PmzQGYP3/+mgSmOurVq8cXX3zBsGHDmDNnDvXq1eO3334D4Ouvv6ZVq1blnrvffvtx5ZVXMnv2bE488US6d+/OyCTqYyv7uhYtWsSBBx7IbrvtxkEHHZT09Tp27Mjnn3/O7NmzOeGEE9Y5p0uXLmy88ca88sornHfeeevs69Sp05qfW7Zsuc6+Fi1arLlHyf02lDyKiKTS//4XJhi+8koYNgzqZEC1TAaEICNHhkSptBpvLwW22morDjzwQB566CEAJk6cyE477USLFi0oLi6mX79+vPjiizz77LNMnTqVjz/+mPz8fFaUU4XVsNSLqhN92levXp2SWMePH0/fvn055JBDmDZtGh9//PGaIe3lxVPiP//5DyNGjOCFF15gn332YZtttqnUcPgSlX1dw4YNo1WrVjz88MN89NFHVbrewoULAWjWrNl65zVr1mzN/tIaN25cqXuV3K+4uLjc40VEUmnGDDjgAMjLg969445mLSVDGWDgQBg1CurXD9vt2oXtOEaTDRo0iGnTpvHRRx8xevRoBg0aBMDMmTN59913OfXUU9l6661rPrAEDz/8MI0bN+bcc89Neu2shg0bcvnll/Ptt9/y2muv0a5dOwYPHsyrr76a0hj/9re/8cwzz7Dtttty6qmnsmrVqqSv0aJFCwB+/vnn9fYtWrRozX4RkUz35Zew//7h59dfh+23jzee0pQMZYiBA0PbaY8e8bSXljjyyCNp0qQJd955J2+88QaHHXYYsLa2JTHxKN2/qCatWLFiTS1KsrEcf/zxQHgt+++/P8888wwAU6dOTWmMO+ywAxtvvDEPPPAAH3/8Mdddd12VrtGmTRsmT568Tvm0adNYvnw5vTPpTysRkXL89luoESoqgtdegx13jDuidSkZknU0aNCAY445hgceeIDDDz+cvLw8AHbccUe23XZbHnzwQRYtWgTAk08+yeeffx5LnH379qWwsHBN89avv/7KLbfcUqlzH3vsMZ566qk122+99RZ169Zlv/32S0eo7Lvvvpx55plcddVVzJgxI6lzzYzrr7+eCRMmrEnali9fzl//+ld23HFHhg4dmoaIRURSq3Hj0Efo1Vdh553jjqYM5Q0zy/VHyXDy8qRjqHHJUPm4vfnmmw54QUHBOuXTp0/3Xr16eatWrbxHjx5+7rnnerdu3bxRo0betWtXX7FihQ8YMMBbtWrlgHft2tWnTp3q11xzjbdp08YB79Spkz/55JP+j3/8w7fbbjsHfLvttvMrr7zS7777bu/UqZMD3qZNGx86dKhPnjzZu3bt6nl5ed6sWbM1w/yLi4v973//u7dp08Z32WUXP+igg/zyyy9fc71bb7213Nf3z3/+0/fYYw/v0qWLd+nSxffcc881w9YXLlzoXbt29UaNGq15XfPnz6/067rpppvWeV033XSTv/76696xY0cHvFWrVn7iiSf6BRdcsM65L774oj/00EPrnFv6NTzzzDPerVs379Chg7dp08YHDRrk8+fPX7M/8Xonn3zyOq+5Z8+e3qxZM8/Ly/OuXbv6woUL/eyzz17nNb399tvlvmcaWi8iVTF7tvukSXFHEVDB0HoL+yVRfn6+FxQUlLv/s88+W2fUTiqUzPGnkc6SadLxeReR7DZ3buj6sWQJfPUVRDOIxMbMprh7fln7NLQ+gZn1A/p16NChxu+tJEhERLLBvHmhs/S8efDKK/EnQhuiPkMJPIZJF0VERLLFggWhs/Ts2TB+fBgclOlUMyQiIiIpc8cd8M038MILsO++cUdTOaoZEhERkZS54gp4//21/WBrAyVDIiIiUi2FhXDccaFprG5d6Nw57oiSo2SoGjQST3KBPuciUpFffoGDD4anngrrjtVGSoaqKC8vj2XLlsUdhkjaLVu2bM3kmyIipf32G/TtC5Mnw+OPh6SoNlIyVEUtW7Zkzpw5LF26VH85S1Zyd5YuXcqcOXNo2bJl3OGISIZZuhT69YN33oF//xuOPDLuiKpOo8mqqEmTJgDMnTuXoqKimKMRSY+8vDxatWq15vMuIlJi2TL49Vd4+GH44x/jjqZ6lAxVQ5MmTfQlISIiOSVat5vNNoP33oN6WZBJZMFLEBERkZqwciUcfXT4eezY7EiEQH2GREREpBKKisLw+eeeg8MOA7O4I0odJUMiIiJSoVWrYOBAePppuO02GDo07ohSS8lQAjPrZ2ajCgsL4w5FREQkI5x9Njz5JNx4Y/g525iGhZctPz/fCwoK4g5DREQkdh9+CG++CeecE3ckVWdmU9w9v6x9qhkSERGR9RQXw/PPh5933712J0IbomRIRERE1uEOw4aFjtITJ8YdTfopGRIREZE13OHPf4Z774WLL4YePeKOKP2UDImIiAgQEqHzz4c77oDzzoOrr86uIfTlUTIkIiIiABQUwC23hBFjN9yQG4kQaAZqERERieyxB7z9NnTvnjuJEKhmSEREJOddey28/HL4ea+9cisRAiVDIiIiOe3aa+GSS+A//4k7kvgoGRIREclRN94YEqEBA+Duu+OOJj5KhkRERHLQbbfBBRfAMcfAQw9B3bpxRxQfJUMiIiI5xh2mT4cjj4RHHoF6OT6cKsdfvoiISG5ZtgwaNIB77oHVqyEvL+6I4qeaIRERkRzx0EOw884wezbUqaNEqISSIRERkRzw6KNwyimw3XbQokXc0WQWJUMiIiJZ7okn4MQTwzpjzz4bmslkLSVDIiIiWeyVV8LQ+b33hnHjoGHDuCPKPEqGEphZPzMbVVhYGHcoIiIi1bbHHjBkCLzwAjRuHHc0mUnJUAJ3H+fuQ5o2bRp3KCIiIlX23nth5Nimm8Jdd8Emm8QdUeZSMiQiIpJlXnkFevaEiy+OO5LaQcmQiIhIFpkwAQ4/HHbYAS6/PO5oagclQyIiIlnizTfhsMPC8PlXX4XNNos7otpByZCIiEgWKCqCQYOgbVt47TXYfPO4I6o9tByHiIhIFsjLC0PnmzWDVq3ijqZ2Uc2QiIhILfbhh3D11WHx1Z13hi23jDui2kfJkIiISC31ySfQuzeMGgWLF8cdTe2lZEhERKQWmj4dDjwwzCj9+uuheUyqRsmQiIhILTNjBhxwQOgnNGECbLtt3BHVbkqGREREaplp09YmQh06xB1N7adkSEREpJYoKgrPf/wjfP55mFhRqk/JkIiISC3w7bdhtNhzz4XtRo1iDSeraJ4hERGRDDd7Nuy/PyxaBFttFXc02UfJkIiISAabMwd69YKffgozS++2W9wRZR8lQyIiIhlq0aIwamzBAnj5ZcjPjzui7KQ+QyIiIhlq003DCvTjx0P37nFHk71UMyQiIpJhFi6EX34J8wddd13c0WQ/JUMiIiIZ5OefwxIbv/4Kn30W5hOS9FIyJCIikiEWL4aDDoJPP4WxY5UI1RQlQyIiIhngl1/g4INh6lR46ino0yfuiHKHkiEREZEMMGIETJkCTz4Jhx0WdzS5RcmQiIhIBrj6aujfP6xELzUr64fWm1lfM7vdzK41s9PjjkdERKTEsmUwfDgsWQKbbKJEKC5ZXTNkZo2Bm4Gd3H2Vmb1rZq+4+1dxxyYiIrlt+XI48sgwmeK++8Khh8YdUe7K+JohM9soqtVZZWbty9h/uJlNNrNJZva2mZWen7M78I27r4q2JwMH10DYIiIi5VqxAo4+Gl56Ce6/X4lQ3DK6ZihKfv4NfAHULWN/N+BRYE93/9TMDgNeMrOd3X0esDnwa6lTfonKREREYlFUBMceC88/D/fcA3/6U9wRSabXDDUGTgQeLGf/JcBL7v4pgLs/B8wHzoz2/whsUur4JlGZiIhILL7/Ht5/H+64A4YOjTsagQxPhtx9urvPrOCQA4GChLLJQO/o5/eAbcyspAZsD+DF1EYpIiKyYatXgztss02YWfrMMzd8jtSMjG4mq4iZNQeaAj8k7JoHHALg7r+Z2XnAbWa2BHioos7TZjYEGALQtm3btMQtIiK5Z/VqOOUUaNkSrr8+LMAqmSOja4Y2oFH0vCKhfAXQsGTD3Z9z92HufqG731PRBd19lLvnu3v+5pura5GIiFRfcTEMHgwPPwzNmoFZ3BFJoqSTITOra2YtzGxLM9s0DTFV1pLouX5CeX1gaQ3HIiIissaYMdC+PdSpA02bwoMPwhVXhFmmJfNUqpnMzA4ABgC/B9pTamSXmS0D/ge8DIyuqTl83P1nM1sMtE7Y1RrQPEIiIhKLMWNgyBBYGv1Z/ttvUK8ebL99vHFJ+SqsGTKz7czsTeAV4A/AIuAF4DHgEeAJ4C1Ck9VFwAwzu9XMamqd3VeB/ISy/Ki8Ssysn5mNKiwsrFZgIiKSm0aMWJsIlVi1SrVCmazcmiEz6wRMAt4gTF442d29guMbAv2B4cB4Mzu41GSH6XItMNHMOrn7Z2Z2KLAFcGdVL+ju44Bx+fn5g1MVpIiI5I5Zs5Irl/iVmQyZWV3CZIdnuvsTlbmQuy8F/m1mjwMPA1cC1cqDzWwjQvPbplHRY2Y2192Piu45xcwGAv8XNdfVBfpEEy6KiIjUqHHjyt+nQcqZq7yaoTzgKHf/OtkLunuxmZ0IdKxWZOFaK4GeGzhmLDC2uvcSERGpqlWr4LLL4NproV07WLAgLMJaomFDGDkyvvikYmX2GXL35VVJhEqdX+zuM6oeloiISO2wahX06RMSoSFDYMYMuO++kBSZhedRo2DgwLgjlfJUetJFM2sC7AosdfeCqGxz4CxCM9bj7v5OGmIUERHJWPXqwQEHwKBBcNJJoWzgQCU/tYlV0Cd63QPN/gpcA4xz9/5mtjHwEaE5zIDVQG93n5imWGuEmfUD+nXo0GHwl19+GXc4IiKSgdzhhhsgPx969Yo7GqkMM5vi7okj0IHkJl3sDwx19/7R9gBgB+BywmKol1HNDtOZwN3HufuQpk2bxh2KiIhkoMWL4aij4KKL4IlKDTGSTJdMMtTW3e8rtT0Q+B642t2XANcDnVIZnIiISCb5+ONQG/Tcc3DzzXDXXXFHJKlQpYVazawlsB9wfcncQ+6+2kwrroiISHaaNg322gs22wzeeAP23jvuiCRVkqkZmmVmh0Q/Xxyd+2jJTjPbCdC0zSIikpV22SXMIv3hh0qEsk0yydD1wFgz+xE4l9CRejqAmV0APAdMSHmENUzLcYiISImZM2H//eGbb8Iw+UsvhZYt445KUq3SyZC7Pw0cRJhdejihzxBR09jOhGU77k1DjDVKHahFRATg6aehWzf45BOYPTvuaCSdkuoz5O4TSKj9ifoMnZLKoEREROJSVASXXAI33gh77glPPqmlNLJdmTVDZlYnmkeoysyscXXOFxERicM//xkSoTPPhEmTlAjlgvJqhuoDU8xssLu/nexFzewfQAvg9OoEJyIiUlNWroSNNoK//AU6d4b+/Td8jmSH8tYmWwYMJXSYHmlmLSpzMTPb3czGESZovCB1YYqIiKRHcTFcfTXssQcsWQKNGikRyjXl9hly9zfNrD+hw/RFZjYNmAHMA5YAqwir2zcDtga6Am2Al4Ae7v5bmmMXERGplp9/DuuJPf88HH98WGZDck+FHajd/a1o/qCTgeOAo4CNyjj0Z+BVYLC7v5LqIGtSqbXJ4g5FRETSqKAAjj4a5s6FO+6AYcPC8HnJPZVeqBXAzOoC7YHNCEnRUmCuu89LS3Qxys/P94KCgrjDEBGRNHCHffeF778Po8X23DPuiCTdKlqoNdmh9auBr6KHiIhIrbJkCaxaBU2bwmOPQcOGYXkNyW3JzEAtIiJSa82YAb/7HZx6athu00aJkARKhkREJOs98UQYLTZ/PpyuSV8kgZIhERHJWitXwjnnwLHHQpcu8NFHcOCBcUclmUbJkIiIZK2ff4bHHw8TKU6cCFtvHXdEkomS6kAtIiJSG7z3XmgWa90a/vc/9Q2SiqW0Zqi665llAjPrZ2ajCgsL4w5FRESStHo1XHEF7L033HVXKFMiJBuS6mayr1N8vRrn7uPcfUjTpk3jDkVERJLw449w6KFw1VVw4olrR42JbEi5zWRmdnkVrqeV6kVEpMa9/36YTfrHH2HUKDjtNM0mLZVXUZ+hv5VR5kBZHy+t5iIiIrFxDwusPvss7L573NFIbVNRMrQQ2KPU9mbAHcDzwATCgq0ArYFewABgSBpiFBERWc+vv8LYsTBwIHTvHjpK160bd1RSG1WUDD3u7t+VbJjZVcAl7v5GwnFfA++Y2VvAX4C3Ux+miIjIWtOnh2axmTPDumLbb69ESKqu3A7U7n52QlGPMhKh0se/AWipOxERSatHHgnLaixeDK++GhIhkepIZjRZczPbqrydZrY10KT6IYmIiJTtL38JI8Xy88Ns0j17xh2RZINkkqFJwAtmdpiZ5ZUUmtlGZtYPeC46RkREJC26doWLLoLXXoMttog7GskWycxAfQ7wJvAsgJktjso3jZ7nAUemKjARERGA558PnaWPOw5OPjnuaCQbVbpmyN2/AjoDNwGfAw2jx+fADUAXd/8mHUHWJM1ALSKSGVatghEj4LDD4I47oLg47ogkW5m7pggqS35+vhcUFMQdhohITpo/H44/HiZMgMGD4bbbYONav+CTxMnMprh7fln7qrRQq5m1BrZy9ylmVsfdla+LiEhK/PxzmDhx0SIYPRoGDYo7Isl2Sa1NZmYHmtnHwBzg9ai4l5l9YmZ9Uh2ciIjknubN4ZxzwsrzSoSkJlQ6GTKzfYEXCDNRjwdWRbsmA2OAx8ysZ4rjExGRHFBYGDpIf/BB2L7oIujSJd6YJHckUzN0GXAvsI27HwasAHD3X9z9OuAE4NLUhygiItnsk0/CvEH//W+YWVqkpiWTDHUDLnT3VWXtdPfnAc0DKiIilfbgg2FdsaVLYeJE+NOf4o5IclEyyZAR1QaVudOsHmGovYiIyAY9/XRIfvbZJ8wmvc8+cUckuSqZZOhrYGAF+4cBX1YvHBERyXarV4fnww+H+++Hl16Cli3jjUlyWzLJ0M3AaDP7PzM7EdgomqDwDDN7Jdp/Y1qiFBGRrPD007DzzvDDD2GV+VNP1WrzEr9KzzPk7o+aWVvg74QaIgOeiZ5XARe5+3/TEaSIiNRuRUUwfDjccEPoLF1UFHdEImslNemiu19rZv8G/gB0iIq/AJ5y91mpDk5ERGq/uXPDsPk334Rhw+Cmm6B+/bijElkr6Rmo3f07wvpkWcnM+gH9OnTosMFjRURkw0aMgClTYMwYGDAg7mhE1lfptcnMrAmwK7DU3Quiss2Bswgr1z/u7u+kJ8yap7XJRETW1bNneJ44ccPHFhfD4sVhNulFi0IfoZ12SmNwIhtQ0dpkyXSgPgOYSJh8ETPbGJgUbZ8NvKEZqEVEZNEi6N8f+vSBlSuhWTMlQpLZkkmG+gND3b1/tD0A2IGQDG0SPY9IbXgiIlKbFBSERVZfeimsK5aXF3dEIhuWTDLU1t3vK7U9EPgeuNrdlwDXA51SGZyIiNQO7nDvvWHixNWrYdIkOOssMIs7MpENS2rV+hJm1hLYD3jEo05H7r6aMMxeRERyzIoVcOut0KsXfPhhWGJDpLZIZjTZLDM7xN3HAxcTEqlHS3aa2U5AYYrjExGRDPbFF7D11tCwIbz+ephJuk6V/swWiU8yH9nrgbFm9iNwLjDO3acDmNkFwHPAhJRHKCIiGemJJ6BbN/jrX8N269ZKhKR2SmYG6qfN7CCgH7AAuB3AzAzYGXgDuDcdQYqISLzGjIH33gvNYe3aQadOoZP0XnvBRRfFHZ1I9SQ7A/UEEmp/oj5Dp6QyKBERyRxjxsCQISERApg1Kzz69IGxY2GjjeKNT6S6kp6B2syaAwcB20RF3wAvu/vPqQxMREQyw/DhsHTp+uUzZigRkuyQVDJkZhcDlwP1WTtyzIEVZnalu/8zxfGJiEgNcw+JzqRJ8MYboRaoLOWVi9Q2lU6GzOwc4GrgI+AFYF60qzVwGHC1mS1z99tSHqWIiKRNcTFMnw6dO4d5gYYNg3vuCfu22CKMFCurZqht25qNUyRdkqkZOgc4x91vL2PfZWZ2NmGUmZIhEZEMtmpVmAuopObnrbfCOmJffAHbbx9WmM/Phx49YLvt4NFHQ5+h0glRw4YwcmRsL0EkpZJJhpqUkwiVuIPQhCYiIhlkxQqYPBnatw9zAj3zDPzxj2Ffx45w9NEh8WnZMpT16BEeJQYODM+nnrp2NNnIkWvLRWq7ZJKheWbW0N3LqCwFoBEwp3SBmbV1d7Uqi4jUoKIiePPNUOszaVIYEr98OdxwA5x/Puy/Pzz+OPz+96EZrDIGDoT7ogWZKrNqvUhtkkwydDtwt5md7u7LSu8wswbArcDfE855D9iyeiGKiEhFfvkF3n4b6tWD3r3DSvF9+oS+QLvtBmecEWp6fv/7cHzz5nDMMfHGLJJJkkmGuhOG1M8xsw+A+VF5K2APwkSMZmZ9S53TNCVRiojIOl58EV55JdT+fPRRSHx69QrJUKNGofZml12gqX4Li2xQMsnQoFI/H1TG/mbADgllnnREMTOzfkC/Dh06xB2KiAgA8+aF5q4vvoBLLw1lt94KEyaEBVEvvRT22y/MBl1in33iiVWkNrJo0fkNH2j2g7tXsnW56udkivz8fC8oKIg7DBHJUW++CQ8/HGp+vvgilG2yCfzwQ6j5mTMHWrSA+vVrLqaePcOz+gxJbWRmU9w9v6x9ySypd1cV7l2Vc0REcoY7zJwJDzwAgwbB7NmhfNq0sBBqx45w3XXw/vvw008hEQLYaquaTYREslkyC7Umdo5ej5kd5O4vJ3OOiEgumjkTLrssNH/NnRvKWrSAU06BNm3CMPahQ6Fu3XjjFMkFSa9NtgGj0egxEZE1Vq8OtTwlExwefnioAdp441DWo0fo77PffmEleIsWOlKtj0jNSXZtsu7AecDOQIMyDtk8FUGJiNR2xcVw1FGhf01hYShr3z7M8QNh8sPvv1+b/IhIfJJZm6wX8DJhCP1MwlD7d6LdbYBtgXdTHaCISCZbvjzM7lxS89O4MTz1FNSpA3l5YT6fkjl+EtfyUiIkkhmSqRm6nLDu2AXu7tFIsV4lO83sciAv1QGKiGSSFSvWNmGddx7cdVcog7DQ6SGHrD32ySdrPr500igyyVbJJENdgYN97Vj8xDH5I4HJwGWpCExEJBMUFobZnSdNCo+PP4b588Mw9513hjPPDDU/++4bZnYWkdonmWRombuvKLW92szqlyorBtqWcZ6ISMYZM6bshUd/+il0bm7UKMzzc/LJof9PXh7ssQecc044Z5NNwvkiUvslM8/Qj2b2+1Lb3wBDS22fDvySkqhERNJozBgYMmRt89Z334Wkp02bMLz9+edD+Z57huHvr70GixeHGqJrrgnHiEj2SKZm6HngeTO70d2vBB4EHjCzIYQms52Af6YhRhGRaisqCjM5T5sWVm5funTd/atWwcKFoYZot91C2Q47wN/+VuOhikgNSyYZuoUwemxxtD2aMMT+VGAFYbbpq1IXmohI8txhyZIwqss9NGV9+CF89llYzb0iK1bA8OE1E6eIZI5kZqD+kVA7VLLtwAXRQ0QkFlOmQEFBqPGZOhWmT4fdd4dXXw1D12fPhi23hD59wmivzp3DxIezZq1/rcSh7yKSG1I6A7WZXeLu16TymiIipZu4pk2DRYvCkHaAiy8Oic8mm4RE55hjYO+91577yivrX+/qq0OfodJNZQ0bhiYyEck9qV6O42xAyZBITGr7quLuYVbm6dPh4INDzc6ll8L1169t4qpXLwxpX706rNt1661h5FfbtpWfxHDgwPA8YkSoIWrbdu1oMhHJPeUmQ2a2uiYDEZHc9MEH8NBDa2t9Fi8O5d99F5KUbt3g3HNDrU+XLqFTc+l1u3baqWr3HThQyY+IBBXVDK0AHk/iWgb8sXrhiEi2KWnimjp1bcIzbRo89hh07w5ffw2PPAK77ALHHru2X0/LluH8I48MDxGRdKkoGSp091OSuZiZ9almPCJSS5U0cZUkO716hXl63n03zNAMoYlrxx1Dn56NNw5lRx8dkiCt0yUicakoGepRhetV5RwRqWUKC2HZMmjdGn7+Gfr3D/18Spq4AK67LiRDu+0Wan46dw6J0EYbrXuteqnuuSgikqRyfw25+xdllZtZHtAN2MLdnzazRsAKd19V3jkiUrs99hh88snaWp9Zs+C00+C++2DTTUOCc9xxoU9P586hyWvTTcO5m2yivjkiktmS+pvMzE4mzDLdAvgNeBroDvw7Glb/QMojFJG0K2niKt2vp3lzuP32sP+yy+Dbb0PNzj77wOmnw377hX116oTlKkREaqtKJ0NmdjjwL8LK9P8Bjo12TQLOAm4xs0Xu/lTKoxSRDRozBt57L8yi3L59+UPFCwtDsjNnTuirA9C3L4wfv/aYNm2gd++126+9FprEEpu4RESygYWJpCtxoNlbwDh3/2e0Pdfdtyy1fy/gRnffu7xr1Cb5+fleUFAQdxgilVKy8GjiJIKjRoWE6LHHwgrs06aFGZkhJDa//RZWY3/88dD3J7GJS0QkW5jZFHfPL2tfMs1knYCe5e1093fNbOskY6sxZrYtcCOhf9NxcccjkkojRqy/8OjSpWF25oEDQxPY7Nnw+9+vHbreufPazsslNUQiIrkomWSowoGvZrYx0KB64aTV74AXgV5xByKSamWtswWhKQzgggvCQ0RE1lcniWOnAedXsP9K4JOqBGFmG5nZtWa2yszal7H/cDObbGaTzOxtMyuzmqsi7v5vwkSSIlmnvAVGtfCoiMiGJVMzNBIYb2Z9gVeABmZ2NrA1cASwLdC7/NPLFiU//wa+AOqWsb8b8Ciwp7t/amaHAS+Z2c7uPi86pqCc13Kou89NNiaR2mbkSBg8OMz9U0ILj4qIVE6lkyF3f9nMTgTuAPaNim8hNJ8tAga4+8QqxNAYOJGQVJ1Uxv5LgJfc/dMojufMbD5wJnBZVJZ0TZFItli2DI4/Pvx86qlhNFm7dlp4VESkspKaZ8jdHzWzZ4E+QIeo+AvgFXdfUpUA3H06QAWdrw8kzG1U2mRCLdRlVbmnSDb5y19gxgx49dUwCSLU3lXrRUTikPRE+FHSUyNzCZlZc6Ap8EPCrnnAIUleqy/QD9jezP7s7reVccwQYAhAW3W2kFrggw/C8PlzztGyFiIiVZXpvz4bRc+JHZ9XAA2TuZC7Pw88v4FjRgGjIMwzlMz1RWra6tVwxhlhMsQrr4w7GhGR2ivTk6GSprf6CeX1gYRZVURyy733wocfwqOPQpMmcUcjIlJ7JTO0vsa5+8/AYqB1wq7WwFc1HpBIhnCHBx6AXr3CAqkiIlJ1mV4zBPAqkDhaLJ8a6rckkonM4K23YNGi8LOIiFRdRtcMRa4F+phZJwAzOxTYArgzHTczs35mNqqwsDAdlxeptq+/DsPpGzSALbfc8PEiIlKx2GuGzGwj4GVg06josWgR2KMA3H2KmQ0E/s/MlhEmZuxTMuFiqrn7OGBcfn7+4HRcX6Q6Vq2CI4+E5s1hwoS4oxERyQ5JJUNm1oowCeLBQDN3b2VmewIDgevcfU6yAbj7SipYADY6ZiwwNtlri2SbO+6AqVPhKTUSi4ikTKWTITNrC7wPtCKM5FoV7ZoPbA9MNrPfu7s6Noukwdy5cPnlcMghcMQRZR+jyRZFRJKXTJ+hfwBfAju7e2Oioe3u/p27H0pYpkOznYikyQUXwMqVcPvt6jQtIpJKyTSTHQh0c/fE2aBLXA98Xv2Q4mVm/YB+HTp02OCxIjVl2TL47ju4+GLYbru4oxERyS7mXrmJls1sobu3KLU91923TDhmvru3SnGMscjPz/eCgoK4wxBZo7g4zDqdlxd3JCIitY+ZTSlvYfdkmskWmdkeFdzkAODnZIMTkYo9/zzMnw916igREhFJh2SSoUeAZ83sT2bWBsDMNjGzTmY2AngSeCgdQYrkqlmz4Jhj4Pzz445ERCR7JdNn6Grgd8D9QEnb2uLo2YBxhH5DIpIi554blt4YOTLuSEREslelkyF3LzKzvoQ5hY4FSrpxfgE87u7/TkN8Ijlr/Hh4+mm4+mpo1y7uaEREslelO1DnilKjyQZ/+eWXcYcjOWr5cthlF6hXL0yyuNFGcUckIlK7paoDdWVuVOv7DLn7OHcf0rRp07hDkRy2dCnsuSfceacSIRGRdEt6bTIzaw1sDTQsY3ffakckIjRvDo8+GncUIiK5IZnlOLYBHgb2Ku8Q1nasFpEqcIdLL4XjjoPOneOORkQkNyRTM3Q/0JqwLMf3wMqE/QbcnKK4RHLSs8+GDtMtWigZEhGpKckkQ92Aju6+oLwDzOyIakckkqOWLIFzzglJ0Nlnxx2NiEjuSCYZ+hL4aQPHHFmNWDKC1iaTuIwcGSZZfPPNMIpMRERqRjKjyYYDl5hZRefMqWY8sdNoMonD55/DDTfAoEGw775xRyMikluSmXTxFTPbAvjOzD4BFgLFCYcpgxCpgjZt4JJL4Mwz445ERCT3JDOabBDwL0JH6a3KOUyjyUSqoGFDuPLKuKMQEclNyTSTXQ7cBmwD1Hf3OokPoNzO1SKyvl9+gf32C/2EREQkHsl002zi7n/ZwDHnVScYkVxzxRXw1ltQv37ckYiI5K5kaoY+iGafroiayUQqaepUuP12GDw4LL0hIiLxSCYZGgpcb2aHmFmLco65KQUxiWS94mIYNgw23TRMsigiIvFJppnsu+h5AICZpT6aDKB5hqQmjB0Lb78N998Pm20WdzQiIrnN3CvXsmVmS4HHKzoE+KO7N0pFYHHLz8/3goKCuMOQLFVcDM88A0ccAXWSqZ8VEZEqMbMp7p5f1r5kaoYK3f2UDdyoT1KRieSgpUvDUPqjjoo7EhERgeT6DPWoxDE7VTUQkVwwZQq0bRtGkImISGaodDLk7l9U4rCR1YhFJKutXg1nnAF5eVqRXkQkkyQzA3XbShz2B2BY1cMRyV733w+TJ8Mjj4CWvhMRyRzlJkNmNgFoDuzu7quBb9E8QiJV8uOPYe2xnj1hwIC4oxERkdIqqhky1m1GWwLcsIHjNQO1SBmeegp+/RXuvBOydFYKEZFaK5mh9T+4+xbVPaa20NB6SbWvvoLttos7ChGR3FTR0PpkRpOV2+XTzBpu6BiRXLRqFXz+efhZiZCISGZKJhm6qoJ9/zSzH4Fav8KSmfUzs1GFhYVxhyJZ4K67YJddYNq0uCMREZHyJNNMNtfdtyxn3+bAwcBwd++Uwvhio2Yyqa4ffoAdd4Tu3eHFF9VXSEQkTqlqJiuXu/8IPEIYfSYiwIUXwvLlcMcdSoRERDJZhfMMmdm/Sm02NbMHCKPGyrrODsDcFMYmUmtNnAhjxsCll8L228cdjYiIVGRDky6eXOpnB8pbm2wp8BmacFEEgOnToWPHMLeQiIhktgqTIXdf04yWTcPmRdLtrLNg8GCoXz/uSEREZEOS6TOkv3FFNmD2bHj55fCzEiERkdohmYVaR6cxDpGs8Je/wBFHhOU3RESkdkjJaDIRgZdegv/+F0aMgM03jzsaERGpLCVDIimwfHnoJ9SxI1xwQdzRiIhIMjY0mkxEKuH662HmzNBfSH2FRERqF9UMiaRAmzYwdCj07h13JCIikizVDImkwMknh4eIiNQ+qhlKoIVaJRnjx8O990JxcdyRiIhIVVV6odZco4VaZUOWLoWddoLGjeGjjyAvL+6IRESkPBUt1KpmMpEqGjkSvvsurEOmREhEpPZSM5lIFXz+eRhBdsIJ0KNH3NGIiEh1KBkSSZJ7mFOoQYOQEImISO2mZjKRJJmFiRV/+glat447GhERqS4lQyJV0KdP3BGIiEiqqJlMJAnDh4eHBmGKiGQPJUMilTR9Olx3XViR3izuaEREJFWUDIlUgjuceSY0bQrXXBN3NCIikkrqMyRSCY88ApMmwahR0KJF3NGIiEgqqWZIZAOKikI/od/9Dk49Ne5oREQk1VQzJLIBeXnw4ovh5zr680FEJOsoGRKpwNKl0LAh7Lxz3JGIiEi66O9ckXIUF8P++4fZpkVEJHspGRIpxwMPwPvvw557xh2JiIikk5IhkTIsXAgXXwy//z2ceGLc0YiISDopGUpgZv3MbFRhYWHcoUiMLrkECgvhzjs1waKISLZTMpTA3ce5+5CmTZvGHYrE5Oef4b//hXPOgc6d445GRETSTaPJRBI0bw4zZkCDBnFHIiIiNUE1QyKlfPZZGEXWsiVssknc0YiISE1QMiQSmT8f9toLLrgg7khERKQmKRkSiVx4YZhkcejQuCMREZGapGRIhLAI68MPw0UXwQ47xB2NiIjUJCVDkvOKimDYMGjXLizIKiIiuUXJUA0aMwbatw+LfbZvH7YlfrNnw/LlcNttYR0yERHJLRpaX0PGjIEhQ0KfFIDvvgvbAAMHxheXwLbbwqefwkYbxR2JiIjEQTVDNWTEiLWJUImlS0O5xOff/w7/DkqERERyl5KhGjJrVnLlkn4vvwwDBsDdd8cdiYiIxEnJUA1p27bscjO45ZYaDUWAFSvgrLNgu+3gzDPjjkZEROKkZKiGjBy5fufcjTeGXXaBL79cW7Z6dc3GlatuuCG873fcEf4dREQkdykZqiEDB8KoUWH4tll4vv9++OSTtTVDkyaFOW4efTQsCSHp8e23ITn9wx/g4IPjjkZEROKmZKgGDRwYvoiLi8NzySiyvLzwXLcuNG4cynfdFcaNA/eYgs1iRUXQowfcfHPckYiISCZQMpRB9tkHPvwwjHBatgwOPzw8JLW23x7Gj4c2beKOREREMoGSoQxTpw4cd1yY92bUKOjfP5SvXg0ffRRvbLXdsmVwzjkwd27ckYiISCZRMpSh8vJg8GA47bSw/fjjsPvucNRRIVGS5F1zTZhl+osv4o5EREQyiZKhWqJfP7jySnj1VejcGQYNCv2OpHK+/BL++c/QH6tnz7ijERGRTKJkqJbYZBO4/HL4+ms47zx44omQIKmD9Ya5w9lnhyH0N9wQdzQiIpJplAzVMi1awPXXw8yZ8MADYZj+b7+FWqNFi+KOLjM9/TS89BL8/e/QunXc0YiISKZRMlRLbbUV7Lln+PnVV0MytM02cPXVITmStfbZB4YPh2HD4o5EREQykZKhLHDEEfDxx7DffmHh1+22g9tv18SNJVq1CpMs1qsXdyQiIpKJsj4ZMrNdzOxRM7vAzO41s1PjjikdunSBsWPh7behU6fQNGQW9uVqv6JPPw2dpb/6Ku5IREQkk2V9MgRsDtzn7jcAZwHXm1mzmGNKm733hgkT4NlnQzL0/ffQrRv897+5lRS5h2axqVOhSZO4oxERkUwWezJkZhuZ2bVmtsrM2pex/3Azm2xmk8zsbTPLT+b67j7B3SeUKioCVlUz7IxmFkafAcybB8uXw9FHwx57wMsv50ZS9Oij8MYbYW6hzTePOxoREclksSZDUfLzBrAFULeM/d2AR4FB7r4fcA3wkpm1LnVMgZl9XMZjyzJueQZwtbv/mo7Xk4ny82HaNBg9GhYuhD59oHdvWJXF6WBhIZx/fkj+SiatFBERKU/cXUobAycCWwMnlbH/EuAld/8UwN2fM7P5wJnAZVFZpWqKzOwPwCbuPjIVgdcmdeuGSRqPOw7uuw9mzVrbmfi776Bdu3jjS7WbboIFC+C558JrFxERqUisNUPuPt3dZ1ZwyIFAQULZZKB3MvcxswFAe3cfaWZdzaxjkqFmhfr14ayz4LrrwnZBAWy7LQwYEOYtyhYXXxw6k+cn1aAqIiK5KvY+Q+Uxs+ZAU+CHhF3zgG2TuE4v4G6gn5lNBMYAZTWhYWZDoma3gh9//LFKcdcm220Hf/1r6Gy9444wdCjMmRN3VFVXXBwWY23QAA47LO5oRESktsjYZAhoFD2vSChfATSs7EWiDtRN3b1n9NjF3SeWc+wod8939/zNc6DXbbNmYZLGr76CM86ABx8Mi8GuSHzHa4nRo2GnnWD27LgjERGR2iSTk6El0XP9hPL6wNIajiWrtW4dJmn84gsYNSo0p7nDPffAr7Wkq/nPP4darq23DrNzi4iIVFbGJkPu/jOwGEhcTao1oGn00qB9e+jfP/z83nuhtmjbbUOH5GXLYg1tg4YPD2uz3XUX1MnYT7WIiGSiTP/aeBVI7AabH5VLGu21F3zwQWg2O/982H77MBItE4fkT54carT+/Gfo3DnuaEREpLbJ9GToWqCPmXUCMLNDCXMS3ZmuG5pZPzMbVVhYmK5b1Bp77BFWe58wAdq2hX/8A1avjjuq9f3rX6Gp729/izsSERGpjcxjnI7YzDYCXgY2BboC7wNz3f2oUsccTphTaBlhYsZz3X1yumPLz8/3goLEUf25yx3mzg39cUpmtD7jDDj00LVroMWluDjMl7TNNvHGISIimcvMppQ3N2Gsky66+0qg5waOGQuMrZGApFxmazsmf/stzJgRhq/vs08YkbbffjUf08KFodmudWslQiIiUnWZ3kwmGWjHHeGzz8Jos2++gR494OCDwzIYNemCC0Ifod9+q9n7iohIdlEyJFWSlxcmaZw5E66/PizvUbI6fE0kRW+9BQ89BKeeCo0bp/9+IiKSvWLtM5SJzKwf0K9Dhw6Dv/zyy7jDqXUWLICOHeEPf4DLL0/PumerVoVRbosXhxqqRo02eIqIiOS4ivoMqWYogbuPc/chTZs2jTuUWqlePfjTn2DMmJAUnXMOzJ+f2nvcfjtMmwa33qpESEREqk/JkKRU8+ZhksYvv4RBg+DOO6FDh1BjlCpffBFGsR1xROquKSIiuUvJkKRFmzZhIsRPP4WRI6Fly1D+3HOwtJqLqdx9Nzz9dPxD+kVEJDsoGZK06tgxzAwNYS6g/v1hu+3CshkrVyZ3rQ8+gKlTw88bbZTaOEVEJHcpGZIa064dvPFGWNrjzDPDEP2HH67crNYrV4Zmt2OOCZMsioiIpIqSoQRajiO99t03JETjx8Omm8KQIZXrYH3TTWGix5tv1kKsIiKSWhpaXw4tx5F+xcUwfTp06RKW+/jzn0Mz2oEHrnvcrFnQqRP06QNPPRVPrCIiUrtpaL1kpDp1QiIEoXZo7Fjo3RsOOADefz8Mz2/fPjSvLVsGPXvGGa2IiGSrWNcmEynRunUYMn/PPWH0WffuULfu2v5E7nDJJbDZZjBwYLyxiohIdlHNkGSM+vXDJI1ffw1Nm67fsXrpUhgxIp7YREQkeykZkozTuDH88kvZ+2bNqtlYREQk+ykZkozUtm1y5SIiIlWlZCiBhtZnhpEjoWHDdcsaNgzlIiIiqaRkKIEWas0MAweG5TzatQvLbrRrF7bVeVpERFJNo8kkYw0cqORHRETSTzVDIiIiktOUDImIiEhOUzIkIiIiOU3JkIiIiOQ0JUMiIiKS05QMJdA8QyIiIrlFyVACzTMkIiKSW5QMiYiISE5TMiQiIiI5TcmQiIiI5DRz97hjyEhm9iPwXTUu0RQorxd2C2BhNa5dnXtn6j1Scb2qXqMq59XEe1xV6f581YRMfn/LEne86b5/6eun4vMV5//3qp4b979xtit5f9P5+6udu29e5h531yMND2BUBfsK4rp3pt4jFder6jWqcl5NvMfVeB/S+vmqodeQse9vJsab7vuXvn4qPl9x/n+v6rlx/xtn+6Pk/Y3r95eaydJnXJbfO9X3SMX1qnqNqpwX579vLqht72/c8ab7/tn0/72q58b9b5ztYn1/1UwWAzMrcPf8uOOQ7KTPl6STPl+STnF9vlQzFI9RcQcgWU2fL0knfb4knWL5fKlmSERERHKaaoZEREQkpykZSgMz28/MnjKziWb2pplNM7NzSu0/zMxeMLPXzOw9MxtvZl3ijFkym5ltZGbXmtkqM2tfqtzM7ITos/SamX1gZk+UPiY6rp2Z/dfMJpvZJDN7xcw61/TrkMxjZieb2Yzo91XpR9NSx5T5+Yv2VeozKNmvst9tZnaamU0xs7ei30XbVXDN/5jZek1YZtbFzF4ys3fN7O3oO7ddVWNXMpQeA4CP3b2nu/8eOAG4ycwOi/aPBh5x9wPcvTvwCfCambWKJ1zJZNGXyhvAFkDdhN2NCJ+nv7n7AcBeQBHwkpk1KHXcI8Bq4Hfuvh/wHPCCmdVPb/RSS1wb/b4q/SiEDX7+oPKfQcl+o9nAd5uZ9QeuBvq6+77As8DLZrZx4sWi78wDyig3wu+wae6+l7vvA3wL/KeqgSsZSo/bgJtLNtz9E2Ax0CEqmuTuj5Y6/kbCRFMH1VSAUqs0Bk4EHixj32rgMXd/E8DdVwO3AB2B0iMydgUmuHtxtP0KsDWwU3pClixS0ecPKv8ZlOxXme+2y4CH3X1etH1vdMzA0hcys0bASODaMu7THGgDvFaq7FUg38yaVSVwJUNp4O6fuvuvAGZWx8wGAyuAJ6P9RyWcsix61l/psh53n+7uM8vZt8zdT0goLuvz9F/gyFJ/qQ8Eiqn9M1VLmlX0+Yv2V/YzKFluQ99tUaLSDSgodU4R8DHQO+HcvwN3A/PLuM9PwETgWDOrZ2b1gOOAJdEjaUqG0sjMLgXmAn8BDnX3OeUcuhewHBhbU7FJVtsLmAdMKlV2KjAbmGtm3wHnA5e4++wY4pPMc5iZvR714fiPme1RzeuV9RmU3JP43bZN9PxDwnHzgG1LNsxsN2BPKh5mfziwGfB99DgCGOruK6sSqJKhNHL3fxDa2UcCb5jZ3onHRG2flwGXuvuCGg5RskzU7n4RcFbCL4XRQHtC1XJ7Qj+2z2o4PMlM84GZhD/Y9gWeBt41s+5VuVgFn0HJIeV8tzWKnlckHL4CaBidVwe4CxhWqlk/8dp1CQlWIeF3WhtCpcM3VY1XyVCaeTCG8BdSWW2fVwPfufuNNRuZZJvol8+/gMfd/b+lynclJD9XuftvHiYXew543Mz2jSVYyRjuPt7dL3b35dH2GOBd4OJkr1XeZ1ByUlnfbSVNWInNp/WBpdHPZwNvufvUCq7dD+gJjHD3oqip7TXg9YpGplVEyVAamNlGZRR/CuyccNy5hA6sp9RAWJL9bgYWufulCeUdo+dvSwqiL775wB9qJjSpZb5i7YCPZJT3GZQcUsF3W0nNTeuE8taEzxyEztb7lEzxQJSUR9vjo2M6Evo8zip1je+APOAwqkDJUHpMif5CKm1LYE2fITM7DTgUOMbdV5nZtmZ2YE0GKdnDzP4ONAPOira7mVm3aHfJ526LUsfXBTZn7V9jkqPM7Boza5hQvBXrftFU5joVfQYlR1T03ebuiwidp/NLHZ8HdCWMBsPd+7r73iVTPBC1qETbh0SnzSHkL6Wno2kVlVXpd5qSofTYhOgXAoRfCsDRwAPR9nHACEJfos5mlk/oSa8mC0mamV1IqDa+C+gWfZ76ASWTKn4A/A84L0qCAM4kVE0/VcPhSubZi9DBHgAz6wHsT/g8VUolPoOSAyr53fYP4MRScw8NBn4CxiRxq+cJNdsXliq7EPgFeLFKsWttstQzswHAaUADwhwcDQiJ0N3u7mZWBNQr49Qr3f1vNRao1ApRs+vLwKaEv6DeB+a6+1FmthVhJEVZTnH30dE12gHXE0ZzrIr2/93dX0hj6FILmNnBhH4amxD+QK4H3ODu/4n2l/v5i/ZX6jMo2a+y321R7dEwQi3OcuD0sqZviJrJWgM7ECb+fMXdR0b7OgPXEWoj6wC/AcPd/b0qxa5kSERERHKZmslEREQkpykZEhERkZymZEhERERympIhERERyWlKhkRERCSnKRkSERGRnKZkSERqjJm9aGY/m5mb2clxx1MVZtbIzO40s1lmNt/MvjSzITHEcZWZzTWzxKUNRCRJSoZEpMa4+8HAUXHHUU1XEtZc6uXurYAHgb1jiKM50JT1F71Mipl9G01uJ5KzypopUkREyncAMM3dSxaWvBFIXNurJpwNXOTuWl9OpJpUMyQikpxmwLKSDXdfES1AWaM8UCIkkgJKhkSyhJkdaGbzzGxl1PTR2cwmmtlCM/vKzP6ScPzDZfXfifr1FEblPRPKS47/k5n9zcy+MbPFZvZ/UV+aFmb2lJktiPrSnFRByHXN7IrouF/MbIqZHVTG62pqZjeZ2ezo/rOjPjublTrm8ui1u5mNNrMjo+v9Upn+SWbWxMxuiN63+Wb2nZndbGZNSx1zipnNA9oAe0f3m1eyIncZ19wh2r8siqGzmb0S9fMpNLMxZtayjPO6mdlzURwLzKwgWu+w9DFvlvNv979S/3YHm9lt0Wv5ycyeNLMWpY49sJzX83ypY440s3fNbE4U97tmdrGZbVrR+ylS67i7HnrokUUPYCKwEHgS2JzwR88lgANHJhzbMyo/OaH85Ki8ZznHfwT0Bwz4HVAE3AncBmwJ1I1+LgZ2KOcanwFnAXlAI8Kq1UVAj1LHNgA+BGYAO0ZlO0bnfgo0Sri2R+X3RtdsDHye+PoSzim5x+fA9lHZ9tH2x0DDhOO/BSYm8e8xOorrA6BrVLYzYXHTqUD9hPdmefReNoje35MICz6PSPLf7j3gwKhsF2AR8GwZ8ZX5eoBehEV9D4u26wBDyvpc6KFHbX+oZkgkO20GjHT3H929GLiJ8MV2RIqu/5W7P+vB+8DbwFDgZXef6+6rgVsIX+Z/KOcac9z9DncvcvclwDmEL9qbSh1zPrAbcLa7zwCIns8HOkX3TNQKONfdl7j7b8C5hESkPCX3GOHuX0b3+BIYQVil/aKK34pKu83dP4mu/z/CitudCQkGZlYHGAUsAc5z92XR+/t/wHjgSjPbLon7vefur0b3mw48BxxiYRX6yjgcKHT356JrFLv7qOg6ap6TrKJkSCQ7LXP3j0s23H0F8COh1iYVEpOLuYTaoNLl30fPW5dzjddLb7j7wuj83c2sJM5jgZWE2q7SJkfPB5dx3cnuXrpPz3h3/7ScGEruAfB8Qvk4QnJ2TAXnJuP1hO0Xoud+0fNuhBqpV6J/r9LGEt7f8hLLsryTsD2bUAu3eSXPXwA0j5oLW5UUuns/d68ouRSpdZQMiWSnhWWUrSR8Gabj+isTy929pKy8kVbzyyibFz1vEz13IIx6nV2qT8s8YBqhBmWzMq6xYAOxJ+oALC6dQMGaBPLnaH8qJL7esl4rwA9lnDs34ZjKKO/fqLKfgdsJCeG5wBwze93MTjWzBknEIFIraGi9SHYqrub5G/pDqczrR01yqfaruyczsWA6YqiNqvU+RE2Mh5tZF2AAcAJwP3ChmfVw97KSWZFaSTVDIrltVfSc+IfReqOc0qBVGWUlSc830fMXQFMza5R4oJntaGZdUxDHl8CmiTUeZlafMLHhlym4B6z/ehNfa8l9ymrK3DLhmLQzszpmZu4+1d0vBtoBlwM7EDq+i2QNJUMiuW1O9NwmoXyvGrh3z9Ib0bDvPYEp7l7SLPR49HxUwrF1gP8AZQ5rT9IT0XPfhPJ+hA7gT5AaPRO2S+43Nnr+mJDsHBglYqUdThhR9lSKYiltCVEybGb1zOwOM2sL/Au4oOSgqFP8zdHmpmmIQyQ2SoZEcpi7fwP8DxhoZltHtQHHA21r4PZdzGyYmeWZWUPgVkLycX6pY24CCoCRZtYNwMwaE/qz1CU021TXjYSpAv5hZttH9+gA/AP4BLg+BfcAODFqcsLMdgYuJPR9ug/WNDEOJkwJcJOZbWzBScAhwBW+dtbrVPoM2NbMNiYkwUNYO6nkmVGsmFk9Qo1QMalLEEUyQ9xj+/XQQ4/UPAjDwOcROsqujn4+kDBfzLyobGX0886lztsOeBlYDMwELgX+RBhJ9TPwfHTcw9G2A4XAs4SEZB7hy9OjnwdEj3lR2bLo5+bAi6WucToh2fka+AWYAhxUxuvahDAM/VtCJ+SvgXuALUodc0YZ9zsjifeuCSEpKrnHLEItSNNSx5xSxvv4fCWuPTqKqx2hFuj76L1+FGhZxvH5hJFtCwgjAAuAgQnHvJnwbzE1Kp8YbZf8290XlU8FfovKfwSuKnWtHQkjz34kNNmdFpXvCtxNmOPpB0In7teAA+L+rOuhR6of5u6IiEh6mNloYJC7W9yxiEjZ1EwmIiIiOU3JkIiIiOQ0JUMiImlQslAr0QzX0YSRf405LBEpg/oMiYiISE5TzZCIiIjkNCVDIiIiktOUDImIiEhOUzIkIiIiOU3JkIiIiOQ0JUMiIiKS0/4f0OcdGfzXfEIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (9, 6))\n",
    "\n",
    "plt.errorbar(ns, tm1, yerr = tm1_err, fmt = 'bo--', label = 'vanilla sinkhorn')\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.yticks(fontsize = 14)\n",
    "plt.xticks([32, 128, 512, 2048], fontsize = 14)\n",
    "plt.gca().get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())\n",
    "\n",
    "plt.xlabel(r\"number of points\", fontsize = 18)\n",
    "plt.ylabel(\"time elapsed (s)\", fontsize = 18)\n",
    "\n",
    "plt.legend(fontsize = 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Measuring performance on stabilized version\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 256\n",
    "\n",
    "x_i = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "x_f = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "\n",
    "C = cdist(x_i, x_f)\n",
    "C /= C.max()\n",
    "\n",
    "a = np.ones(n_samples) / n_samples\n",
    "b = np.ones(n_samples) / n_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "603 ms ± 26.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "799 ms ± 241 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "971 ms ± 160 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "807 ms ± 101 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "719 ms ± 51.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "672 ms ± 52 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "676 ms ± 90.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "583 ms ± 63.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "404 ms ± 54.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "245 ms ± 94.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "91.7 ms ± 2.73 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "58.1 ms ± 3.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "36.9 ms ± 3.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "26.1 ms ± 2.41 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "20.9 ms ± 2.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "15.3 ms ± 491 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "17.6 ms ± 1.6 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "17.4 ms ± 1.41 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "17.9 ms ± 909 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "18.4 ms ± 1.76 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "regs = np.logspace(-10, 0, num = 20, base = 2)\n",
    "\n",
    "tm1, tm1_err = [], []\n",
    "\n",
    "for reg in regs:\n",
    "    t = %timeit -o sinkhorn.sinkhorn_stabilized(C, reg = reg, a = a, b = b)\n",
    "    tm1.append(t.average)\n",
    "    tm1_err.append(t.stdev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9c486ac970>"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGCCAYAAAAMp6wpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABGwklEQVR4nO3deXxU1f3/8deH3YCgyCZrUJQKIiApdWkLKoIb4toqQZFiwfqrXa3Ll7pQwb3WBTewigvVqlULuFSogq2VFlRUXKoFAQHZRAEJsn5+f5xJSIaZMJNMcpOZ9/PxuI9kzj333s8Ml+STc89i7o6IiIhILqoTdQAiIiIiUVEiJCIiIjlLiZCIiIjkLCVCIiIikrOUCImIiEjOUiIkIiIiOate1AHUVC1atPD8/PyowxAREZEMePPNN9e6e8v4ciVCSeTn5zNv3ryowxAREZEMMLMlicr1aExERERylhIhERERyVlKhERERCRnKRESERGRnKVESERERHKWEiERERHJWRo+LyISkQ0bNrB69Wq2bdsWdSgitVb9+vVp1aoVTZs2rdDxSoRERCKwYcMGVq1aRbt27dhrr70ws6hDEql13J3NmzezfPlygAolQ3o0JiISgdWrV9OuXTvy8vKUBIlUkJmRl5dHu3btWL16dYXOoURIRCQC27ZtY6+99oo6DJGssNdee1X4EbMSIRGRiKglSCQzKvN/SYmQiEgW6d8/bCKSGiVCkjL9gBURkWyjREhERGq0efPm0atXLxo0aMAFF1yQ9vEnnXQSbdq0ifRR5EMPPUS3bt0wMyZPnhxZHFOnTqVXr16YGddee23Sevfdd1+NiLc6KBESEZFqce211zJr1qy0jysoKGD+/Pm0bds2pfqHH344v/71r0tev/DCC1x00UW71bv99tvJz8+nqKgo7ZjSNWLECF544YWMne+MM87gzDPPTPu4U089lfnz5++x3kUXXZTReGuyrJ9HyMxOBk4ANgGL3f2+iEMSEclJY8eOBaB/FT9j79ixI61atdpjvebNm9OxY0fq1q1bpfFUhbZt26qzfYZkdSJkZk2APwDd3H27mb1hZjPcfWHUsYmIZNqUKTBnDmzZAvn5MH48FBZGHVX1e+6551Kqd/7553P++edXbTBVZMKECVGHkDVq/KMxM2tgZjea2XYzy0+w/1Qzm2tmr5nZ62ZWUGr3EcCn7r499nouoXVIRCSrTJkCo0aFJAhgyZLwesqU6rn+hx9+yEknnUSvXr3o1asX/fr144EHHgDg3//+N7169QJC35PiOh988AEAEydO5Oijj6agoIAePXpw8skn89FHHyW91o033kjfvn1p2bIlAwYM4H//+1/Jvu985zs0b96c/Pz8cuMdP348Xbp0wcxKHtctXryYJk2alMTXq1cvDj74YMyszGOodevW8eMf/5hOnTpx8MEH07dvX1588cXdrnHTTTfRrl07DjnkEE466aQycZZn27ZtXH755Rx66KH07t2bXr168atf/Yo1a9YAcPrpp+/W5+nqq68ueT9PPfUU55xzDt26deOggw7iL3/5S7nX++9//8uhhx6KmXHYYYdxyy23lNm/adMmRo0aRc+ePcnPz0+YhC1YsIBTTjmF/Px8OnfuzMCBA3nrrbdK9pfuc3TPPfcwevRo+vTpQ926dfnFL37BhRdeSMeOHTEz/v73v3Pqqady0EEH0aNHD1577bWUPrcKc/cauwH5wBvAw4AD+XH7+wBfE1p8AE4BvgDaxF6fCzxdqv444NpUrt2nTx+vbfr1C1ttPb9ILvnggw+S7iv+v1Z6u/vusG/TpsT799vPHXbfGjYM+594Ihy/dGni46dODfs/+qhi7+fAAw/06667ruT1hAkTvFOnTmXqAH7NNdfsduy3vvUtf+mll0pe//73v/d27dr5hg0bytTr1KmTt2zZ0idNmuTu7t98840PHDjQO3fu7Fu2bCmpN3z48N2ufc0113j4lbfLq6++6oC/+uqr7u7+6aefer+4H3LDhg3zRo0alfx7ffPNN967d2//zne+4+vXr3d39yeffNLr1q3rr7zySslx9957r9etW9dffvlld3dfs2aNH3vssQ74Qw89tNtnUNp1113nPXr08KKiInd3X7Rokbds2bIkzj29n+9///u+du1ad3e/8sorvXHjxr5u3boydUv/W2zfvt1/+MMf+g033FCmzqeffuqA9+jRwxctWuTu7vfff7+bmX9U6kb55JNPvGnTpv7rX//ad+7cWRJf48aN/d13393tfF27dvX33nvP3cO/9c9//nN3d3/ooYcc8LPOOqvkvQ8dOtQ7dOjg27ZtK/czcy///1TsPc/zBL/va3qLUBPgPOChJPuvBP7m7h8AuPt0YBXw/2L71wB7l6rfNFYmIpJVvvgicXlxC1FVWrt2LQsXLuSAAw4oKbvwwgv50Y9+lNLxzz77LIMGDSp5ffHFF7N8+fKEnXWbN2/OhRdeCEDDhg0ZN24cn376KY888kgl3wW0aNGC3/zmNyWvn3rqKR577DFuvPFGDjnkEAAeffRR3n77bcaNG1eyrtXZZ59NQUFBSR+onTt3Mm7cOI4//niOP/74knOn+nnMmTOHNm3alMw83rlzZ2666Sbat2+f0vGnn346++23HxA6VW/atIm5c+cmrLtz506GDx9Ot27duOKKKxLWOfbYY+ncuXPJ+dyd2bNnl+wvHn123XXXlbRSjRkzhsaNGzNmzJjdznfcccdx6KGHAjB69OjdrltYWFjy3k8//XQ+++wzFi1alNJ7r4ga3UfI3RcAmFmyf/0BwE1xZXOB44GrgDlAZzOr5+Hx2LeBO6oo3Kymvgci1ae8gVV5eYn35+eHx2HxOnUqW79Dh/LP37VrSiGWsd9++9G9e3dGjx7N/PnzGTp0KL169eLqq69O6fidO3cydOhQFixYQJ06u/4+T/TLr1u3bmVe9+7dGzPjjTfeKEmQKqpJkyacfPLJAKxYsYKLLrqIAQMG8LOf/aykzsyZMzEzjjrqqDLHHnrooTzyyCNs27aNFStWsHz5coYPH75bnVR8//vf5/LLL+fkk09m1KhRDBw4kBEjRqT8Pg4++OCS75s3bw7AqlWrdqvn7owcOZI5c+bw6KOPVvh8M2fOpHv37mWWjKlfvz69e/dm5syZuHuZx3jFSSVA48aNady4cUrXK12eSTW9RSgpM2sONAM+j9u1EjgAwN2/Bn4F3GlmtwAPezkdpc1slJnNM7N5xc9iJfq+ByKyZ+PHhySptLy8UF7VivvZjB49msmTJ9O7d2+6devGX//61z0eu3z5cr73ve+xc+dO5syZw/z580uGd29J0JwVv7p4vXr12HvvvVmxYkVG3guEBGHEiBG4O5MnTy7zS3zt2rUliVDpvkSvvvoqzZs358svv2TlypUA7LPPPmXO26xZs5Suf9lll/Hwww+zevVqTjvtNFq3bs2ll16a8PNIJK/UjVCcWO7YsWO3ehMmTOCLL75g4cKF3H333RU+39q1a9l33313O6558+Zs3rx5t+kJmjRpkpH4M6XWJkJAcQoZf2dsAUo+RXef7u4Xu/tvfA9D5919orsXuHtBy5YtMxxu7TVmDMRPs1FUFMpFpGYoLISJE6Fhw/C6U6fwurpablu0aMGtt97KihUreO6556hbty5nnnlmuZ2eAZ5//nnWrVvHZZddVuYXYDIbNmwo83r79u1s3Lgx5TmGUjFhwgRefvll7r33Xtq1a1dmX4sWLahTpw5vvfVWSdI2f/58Fi5cyMqVK2nVqhX7778/AF9++WWZY7/66quUYzj//POZO3cu77//Pj/4wQ/4/e9/z3XXXVfp91bawIEDmTp1Kueddx5XXnklSxI1KaagRYsWrFu3brfydevWsddee6X07xql2pwIbYp9bRhX3hCo+tmxcsjSpemVi0g0CgvhiCOgXz9YvLj6kqDVq1fz85//HAgtNEOGDOHhhx9mx44dJSPDiveFPqthlNn8+fNLWjlKt7oUt6gk8v7775d5/dZbb+HuHHnkkRl5Lx999BGXX345Q4cO5Yc//GFJ+fXXXw/A8ccfz/bt23eL4+2332b06NEAtG/fnvbt25cZNZUo9mSuvPJKFi9eDIRHgQ888AA9evTg3XffrejbSqhr7Dno7bffTuPGjRk1alSFzjNgwADef//9Mi0/27dvZ/78+QwYMKDGz3dUaxMhd18HfAW0idvVBtA8QRnw8cfha8eOifd36BC+fvZZGJ8iIrmpqKiI++67jzlz5pSU/fOf/6Rx48b07du3pKxz584sW7YMgHHjxjF16lSOP/546tevzx133MGOHTtwd8aNG5f0WsuWLWPSpElAeHR21VVX0blz54zMB7Rt2zaGDRtGixYtdntUNHHiRADOO+88+vTpw69//Ws2btwIhJaPSy65pCSxqFOnDldffTUzZ85kxowZAHzxxRfcddddKcXxxhtvcMcdd7Bz504AlixZwrJlyzj22GMr/R4Tad68OXfffTcvv/xyhZbTuOaaazAzrrrqqpJE9/rrr2fjxo2Mr45ns5WVaChZTduA/iQePv8UpYbHx8reB66r7DVzefj8e++5n3CCe7167h9/7P7YY+55eWWH5OblhfKvv3Zv1cr96KPdZ8+u/LVFcsWehvpWVBTTXBQVFfnVV1/tvXr18l69evmhhx7q/fr189dee61Mveeee84POOAA79Gjhx9zzDG+evVqd3d/9tlnvXv37t6pUyfv37+/T5gwwQFv3bq1n3feeT537lzv2bOn169f38877zy//PLL/fDDD/cWLVr4cccd55988knJNfr27ev77ruv169f33v27Olvv/22n3jiid66dWsHvGfPnv7888/7uHHj/MADD3TADzzwQB87dqw//vjjDnjTpk29U6dOu23FvvzySx89erR36NDBDzvsMO/Tp4/fXTy/QSk333yzt23b1rt27er9+/f3Z5991gHv0KGDjxw5Munn+dxzz/nAgQO9e/fu3rNnT+/evbvfcMMNJUPTTzvttN3ez2233Vbm/Tz44IM+Y8YMP+SQQ0quefnll/trr73mPXv2LPl8hw4d6hs2bPCePXt6nTp1vHHjxt6zZ09/4oknyhw7btw4f//998sce95555XE/N577/lJJ53kHTt29E6dOvmAAQN83rx5Jfvjz9ezZ0/fvn17yf5LL73UO3To4IAfcsgh/tJLL/nDDz9c5j3dcccd5d6HFR0+H3mSk8pWTiLUB9gIHBJ7fRKwjtg8QpXZcjERWrHC/cIL3evUcd9nH/dbb3X/5puw77HHwnwk4N6pU3jt7r51q/t997m3bRv2DRrkXureF5EksikREqkJKpoI1ejh82bWAHgZ2CdW9ISZrXD3MwDc/U0zKwQeMbPNQF1gkLsnf8AsCW3cCN27w9dfw89+Br/9LcSmoQBCX4NYa3SZobf168Po0XD++XD33XDDDVBQAG+/DbGJZFNWvPxQBdZkFJEY/f8RSU+NToTcfSuhNai8OlOBqZm6ppkNBgZ36dIlU6essXbsgBkz4IQTYO+94Q9/gKOPhoq89b32gksvDcPqn3oKevYM5U8+CX37hjlOREREappa21m6qrj7NHcflep8D7XV3/4GvXvDiSfCv/8dyoYPr1gSVFrTpjByJJjBpk2htejgg+GSS6CcgSAiIiKRUCKUY959FwYNCq1AX3+9q8WmKjRuDO+9ByNGwL33woEHwv/9H8RNrSEiIhIZJUI55JtvYMAAmDsXbrsNPvwQzj47tN5Ulfbt4f774aOP4LTT4JZb4PP4ucBFREQiokQoy339Ndx1V+gP1KgRPP00/O9/8Mtf7pqBtjp06RKW5Fi8GIqXCrr4YrjzzrB0R/FaZrNnh/5EWr5DckEYyCIilVWZ/0s1urO0pC5+UdTrrgvfX3VV6JvTvTsceyx8//vRxlk8W/2WLfDf/4ZHZmPHhoRt69awr3gtM9DCrpK96tevz+bNm2v88gMitcHmzZupX79+hY5Vi1AWSLQo6vDh8OMfwwEHwL/+FZKgmqRhQ5g5M4xaK50EFdNaZpLtWrVqxfLlyykqKlLLkEgFuTtFRUUsX76cVq1aVegcpv+AiRUUFPi8efOiDiMl+fkh+YnXogWsXl21fYAyoU6dxEt0mMHmzdX7CE+kOm3YsIHVq1ezbdu2qEMRqbXq169Pq1ataNq0abn1zOxNdy+IL9ejsTi1cR6hZIuffvFFzU+CIKxlliiRc4fWreHUU+Gss2DgwNDPSSRbNG3adI8/vEWkaunRWJzaNo+QOyQLNdliqTXN+PEQ300iLw9+8xs44wyYPh2GDIHBg3ft3769emMUEZHspBahWmzTJvjRj+Crr6Bu3TAyrFheXkgwaoPiDtEjR4Z+Tp06hdiLy++/H155JTxCg/B+u3SB448Pw/9PPDHMbC0iIpIutQjVUosWwZFHhuHwN90EDz+8qy9Np04wcWLtGnFVWAhHHAH9+oUh9qVjr18/TAJ5/PHhdVFRaCmaMQPOPBNatoQf/jDMi1Se/v13rWdWW2XDexARqUmUCNVSO3aEFqEXX4TLLis/kcg2bduGRG/lypAMDRsGr74aWsUgTCPw5JPh8ymmeYpERCQRPRqrRdxh6tTQefigg8JszRWcNiEr1KsXZsoeMCCsfF+cCP3xj/DAA+Fx2Uknwf77h7LS0wtoniIREQG1CNUamzbBOeeEZSqmTg1luZwExStOggDuuy+0EP3oR/D66zBhQhiGX5rmKRIREVAiVCvE9wc69dSoI6rZ6tYN/WgmTIBly5JPIbBkCfziF6F/1aJFlb+u+u+IiNQ+SoRquJkzoaAg/EIv7g9UG+YGqinq1k0+jUCDBqGv0QUXwOTJoWzDhvDY7N57Q5+ioqLqilRERKKgRCiOmQ02s4nr16+POhQg9Avq1CmsGD9wYNTRVK1Zs8KWacnmKXrwQdi4ET74ICxHAqGj+dNPhwVhjzwS9t47LBI7c2bYv2kTrFu3+zWqozO2OnyLiGSeOkvHcfdpwLSCgoIfRxXDpk3hF++QIWHI+Lx5ZfvASHr2NE/RIYfsqnvYYWFG7qVL4e234a23wtf99gv7p08PfbU6dYLevcO2fn1oQcp0Z+ydO8PEkTt2wDPP7L6enDp8i4hUntYaSyKqtcYWLQodoj/8EP73v/ALN1XF/VOqolUlG2Ti8/n4Y3juuV0J0iefJF4nDcK8TocfHpKZOXPChJBjx4aWnO3bd20NGoSWKIARI3btLz5v69ZhaZFEy5B06rTrWBERSU5rjdUCM2aEiQEBnn8+vSRIqsfBB4d+WsW+/hqaNk2cDG3ZAk2ahGH+O3eGRKh9+9Dnq169XVvp9dNOOCEM9y+9f++94ec/TxzPkiWhBau4xUpERNKjFqEkqrtF6Lbbwtpa3bvDs8/CgQemfw61CEUjP7/qW2uSXQPCfXPzzSEZU0d6EZHEkrUIqbN0DVGnTlhh/Y03KpYESXSSdcbO5Fpvya5x441hCgAILYpHHw2PPw5bt2bu2iIi2UyJUIQWLQqLiUJ49PHEE9C4cbQxSfoKC8Mw/Kpc6y3ZNS6/PCw5AuFR3Jo1MHRomDLgmmtgxYr0rqO5kEQk1ygRqiZTpoTHG3XqhK9XXBH6iowcCdu2hUcaeqxRe1XHWm97usbgwWHZlRdfDPfWddeFFqKdOzMfi4hItlAiVA2mTAlDnZcsCf04liwJM0Q3bhyGyWupDMmUOnVCh+vp08OItkmTQtm2baH8gQc0SaSISGlKhKrBmDGJf/mYqT9QNqmqCSEr6sADw4K0EGYm//zzMHFk+/Zw6aWZWVZERKS2UyJUDZYuTVy+bFn1xiG5q3NnmD8fXnstTNJ5++3QpUuYD6mYZq4WkVykeYTimNlgYHCXLl0yds6OHRMPfU62BpZIVTCD730vbMuXw5//DL16hX2FhfDUU+ERGmjmahHJHWoRiuPu09x9VLNmzTJ2zuoYXi2Sjnbt4Fe/Cv2HICzhUZwEFSsqCo91RUSymVqEqkHxX9RjxoTHZB07ll3rSiRVVdUHqXgNs3jJHuuKiGQLJULVpLBQiY/UXHp8KyK5So/GRESPb0UkZykRyiI1bfi21B7JZq6uVy+MNhMRyVZ6NCYiQEiGJk0K38+aBZs3Q7du8PXXYUh9t26RhiciUiXUIiQiCe21F7z8cmgVOu64MFO1iEi2USIkIkkddBD8/e+wfXtIhhJ1qBYRqc2UCIlIubp1gxkzYOPGMN+QiEg2UR8hESmRrLN9r17w/vvQtm11RiMiUvXUIiQiKSlOgubPD+uVrVsXaTgiIhmhRCiOmQ02s4nr16+POhSRGmntWvjHP2DQINB/ExGp7ZQIxamKtcZEssmAAfD006Fl6OSTYdOmqCMSEak4JUIikrZTToHHH4c33oBTTw1zDomI1EZKhESkQs46Cx5+OMw3JCJSWykREpEKGzYMpk0LydCGDbBtW9QRiYikR4mQiFSKGWzZAsccAxdcADt2RB2RiEjqlAiJSKU1bAhnnw1/+hOMHg07dyau179/2EREagpNqCgiGXHFFVBUBNddB40awV13hdYiEZGaTImQiGTM2LEhGfr976FdO7jyyqgjEhEpnxIhEckYM7jlFmjWDIYOjToaEZE9Ux8hEckoM7jqKujUKfQV+tvfoo5IRCQ5JUIiUmUeeABOOCH0FxIRqYmUCIlIlRkxAk47DX72M7jwQpgzB2bPhvx8mDIl6uhERJQIiUgVql8fnngCDjsM/vjHMN8QwJIlMGqUkiERiZ4SIRGpUg0bwpdf7l5eVARjxlR/PCIipSkREpEqt2xZ4vKlS6s3DhGReEqE4pjZYDObuH79+qhDEckaHTumVy4iUl2UCMVx92nuPqpZs2ZRhyKSNcaPh7y8smUNGoRyEZEoaUJFEalyhYXh68iRocN0w4ZhvqHvfjfauERE1CIkItWisBCOOAL69YMPPwwjys4/X6vVi0i0lAiJSLXr3BnuvBNeew1uuy3qaEQkl+nRmIhEYvjwkAi1bRt1JCKSy9JOhMysLrAv0AAocvevMh2UiGQ/M3jwwaijEJFcl9KjMTM7zsz+aGYfA5uBVcBnwBdm9rWZ/dvMrjOzA6syWBHJTvfdp8kVRSQa5bYIxRKbycDRwAbgv8AHwCZgO6FVaF+gPXAZcIWZ3QNc6u7bqi5sEckmCxbA3XfDMcfAgAFRRyMiucTcPfEOs0OA14DZwM3AXE9WOdTPA4YA/0doMTrB3bdnPOJqUlBQ4PPmzYs6DJGcUFQEffrAxo3w7rvQvHnUEYlItjGzN929IL484aOxWD+gx4H/5+5nuft/ykuCANy9yN0fB3oSEqGxGYhbRHJAXh489hisWgUXXwzl/7QREcmcZH2E6gNnuPuT6Z7Q3XcC5wGPViYwEcktffrA2LHw5z/Dm29W7Bz9+4dNRCRVCRMhd//G3RdV9KTuvtPdP6p4WCKSiy6/HF5/HQp2a7wWEakaKU+oaGZNzez7ZlZQqqylmY01szvM7KiqCVFEckXdunBU7CfJu+/Czp3RxiMi2S+dmaV/AswCrgIws0aEztS/BS4BZptZ/8yGJyK5aP586N1bs06LSNVLJxEaAox29yGx10OBrsDVwN6EBEkzgYhIpfXsCUOGhLmF3n036mhEJJulkwh1dPdJpV4XAsuA6919E3ALcEgmgxOR3GQGEyeGYfSFhfDNN1FHJCLZqkKLrppZK+D7wGPFw+rdfQdgGYxNRHJYixZhCY4FCzTrtIhUnXQSoaVmdmLs+ytix/6peKeZdQPWZzA2EclxJ54Iv/gFtGwZdSQikq3SWXT1FmCqmX0F7AdMdfcFAGZ2KXAx8GLGI6xmZjYYGNylS5eoQxER4A9/iDoCEclmKbcIufuzwEDCRIn/R+gjhJkZ0J2wFMf9VRBjtXL3ae4+qlmzZlGHIiKlTJsGl1wSdRQikm3SaRHC3V8FXo0rc2BEJoMSEYk3fz5MmBDmGTr33KijEZFskWytsTqxeYIqzMyaVOZ4EZHSrrwSjjwSfvIT+OyzqKMRkWyR7NFYQ+AtMzu6Iic1s3HArRWOSkQkTr168OijsGMHDB+uWadFJDOSrTW2GRhN6Bw93sxapHIyMzvczKYRJl+8NHNhiojAgQfC7bfDq6/C1KlRRyMi2SBpHyF3/4eZDSF0jr7MzN4DPgJWApuA7YRV6vcF2gM9gQ7A34B+7v51FccuIjnoRz+Cjh1hwICoIxGRbFBuZ2l3/2dsfqALgHOAM4AGCaquA2YCP3b3GZkOUkSkmBkcf3z4fuFCaNcOGlWqR6OI5LI9jhqLPSa7F7jXzOoC+YR5hBoARcAKd19ZlUGKiMRbsSKsSXbRRXCreiSKSAWlO3x+B7AwtomIRKZtWzj//LBC/cknwzHHRB2RiNRGFVprTESkJrjlFujSJSREX34ZdTQiUhspERKRWqtxY3jsMfj889AqNGcOzJ4N+fkwZUrU0YlIbaBESERqtb594bTTYO5c2LIllC1ZAqNGKRkSkT1TIiQitd7cubB9e9myoiIYMyaaeESk9lAiJCK1XrIlN5Yurd44RKT2yWgiVNn1yUREKqJjx/TKRUSKZbpFaFGGzyciskfjx0NeXtmyvLxQLiJSnqTzCJnZ1RU4n1acF5FqV1gYvo4cGTpMN2kC9923q1xEJBlz98Q7zBKt7eyAJSkP37jXzUxo0SooKPB58+ZFHYaIpKF/f/jgg7Ay/eefQ/36UUckIjWFmb3p7gXx5eXNLL0W+Hap1/sBE4DngVcJi68CtAGOAYYCozISrYhIBbVuDRs3wqJF0LVr1NGISE1XXiL0Z3dfUvzCzH4HXOnus+PqLQL+ZWb/BH4JvJ75MEVEUtO8ObzzDtTRmFgRSUHSHxXufklcUb8ESVDp+rOBvpkKTESkIsxCErR9O2zdGnU0IlLTpfM3U3Mza5dsp5m1B5pWPiQRkcpZujQsyvrEE1FHIiI1XTqJ0GvAC2Z2ipmVdEE0swZmNhiYHqsjIhKpDh1g773DOmQiIuVJJxH6OdAS+CvwjZl9YWZfAJuB52L7fp7xCEVE0mQWhs7//e+wYkXU0YhITZZyIuTuC4EewG3Af4G82PZf4FbgMHf/tCqCFBFJV2FhGEavx2MiUp60xlW4+xfu/ht37+bue8W2bu5+ubt/UVVBioikq2tXKCjQCvQiUr7yhs8nZWZtgHbu/qaZ1XH3RJMviohE6oYboF49cA+Py0RE4qXVImRmA8xsPrAceCVWfIyZvWNmgzIdXBTMbLCZTVy/fn3UoYhIJQ0YEGabVhIkIsmknAiZ2XeBFwgzTL8IbI/tmgtMAZ4ws/4Zjq/aufs0dx/VrFmzqEMRkQz48EO45prQX0hEJF46LUJXAfcDnd39FGALgLtvcPebgWHAbzMfoohIxc2bB7/7HfzrX1FHIiI1UTqJUB/gN+6+PdFOd38eOCgjUYmIZMjpp0NenjpNi0hi6SRCRqwVKOFOs3qE4fQiIjVGkyYwZAg8+aSW3BCR3aWTCC0CCsvZfzHwSeXCERHJvGHDYN06eOmlqCMRkZomneHzfwAeMbOBwAygeGmN9sAZwLHADzIfoohI5Rx/PBx4ICxfHnUkIlLTpJwIufufzKwjcB2hZcgIS2sYYQTZZe7+l6oIUkSkMurXh48/DqvSi4iUltaEiu5+o5k9DpwJdIkVfww84+5LMx2ciEim1KkTJlbcuBGaNo06GhGpKdKeWdrdlxDWGxMRqVWGDIGiIpg5M+pIRKSmSGdCxaZm9n0zKyhV1tLMxprZHWZ2VNWEKCKSGYcfDq+8or5CIrJLOk/MfwLMIkysiJk1Al6Lvb4EmJ0NM0uLSPYqLAyPxx5/POpIRKSmSCcRGgKMdvchsddDga6ERGjv2NcxmQ1PRCRzDjoI+vbV5Ioisks6iVBHd59U6nUhsAy43t03AbcAh2QyOBGRTCsshPnzYcGCqCMRkZog7c7SAGbWCvg+cIu7O4C77zDTGs8iUrOdey60aAEHHBB1JCJSE6TTIrTUzE6MfX9F7Ng/Fe80s27A+gzGJiKScS1bwtChYf0xEZF0EqFbgKlmtgb4BTDN3RcAmNmlwHTg1YxHKCKSYRs2wM03h5XpRSS3pTOz9LOx5TUGA6uBuwBij8O6A7OB+6siSBGRTKpbF373O1i4EAoK9lxfRLJXujNLv0pcq0+sj9CITAYlIlKVGjeG004LK9LfeSc0bBh1RCISlbQ7S5tZc2Ag0DlW9Cnwsruvy2RgIiLpmjUr9brDhoVh9C++GJIiEclNaSVCZnYFcDXQkLDYKoADW8xsrLvflOH4RESqxIAB0KoVPPaYEiGRXJZyImRmPweuB94GXgBWxna1AU4Brjezze5+Z8ajFBHJsHr1QqvQwoVhtmlN/iGSmyw2DdCeK5otAv7g7ncl2X8J8At3PzCD8UWmoKDA52lIiUhWUwIkkjvM7E133214RDrD55smS4JiJgBN045MRCQixUnQF19EG4eIRCedRGilmZU3BVljoMyazmbWsUJRiYhUk0cegTZt4LPPoo5ERKKQTiJ0F3Cvme0VvyNWdgdwXdyuOZWITUSkyh19NGzfrhXpRXJVOqPGjiAMm19uZv8BVsXKWwPfJkyyaGZ2cqljmmUkShGRKnLggXDEEWEo/WWXRR2NiFS3dBKh4aW+H5hg/75A17iy1Hpii4hEaNgw+OlP4d134bDDoo5GRKpTOo/GVrl7nXQ2QiuRiEiN9oMfhGU3pkyJOhIRqW7ptAjdU4HzV+QYEZFq1bIl/PnPcOSRUUciItUt5RYhd4/vCL2b2KKsaR0jIlITnHkmtG0bdRQiUt3SeTSWiskZPp+ISLV56qmwCKuI5I501xo7AvgV0B3YbRg90DITQYmIRGH6dPjrX2HUKGjUKOpoRKQ6pNwiZGbHAP8AjgbWAu2AJbFtJ5AP/CfzIYqIVI/CQli/Hl54IepIRKS6pPNo7GrgTqC9u/cD1rn7MbGtC3At8EoVxCgiUi2OPTbMMv3YY1FHIiLVJZ1EqCfwf75rldb4OYLGAycjIlJL1asH55wDzz8PX34ZdTQiUh3SSYQ2u/uWUq93mFnDUq93AlpbTERqtWHDoGtXWLo06khEpDqkkwitMbPvlXr9KTC61OuLgA0ZiUpEJCKHHx5mmO7ZM+pIRKQ6pJMIPQ88b2bXxF4/BNxuZgvM7D1gAvBEpgMUEalOZuHrpk2wQX/aiWS9dBKh24FzgZmx15OB2wijx1oSZpH+XQZjExGJxLp1odP0PZobXyTr2a6+z1JaQUGBz5s3L+owRCQiRx8dhtK/996uViIRqb3M7E13L4gvz+jM0mZ2ZSbPlylmdoCZPWtmenQnIikZNgzefz/0FxKR7JXpJTYuyfD5MuU7wEtRByEitcfZZ4fh9JpTSCS7JU2EzGxHuhvQOt0AzKyBmd1oZtvNLD/B/lPNbK6ZvWZmr5vZbs1ae+LujwNb9lhRRCSmRQs48UT4059gx46ooxGRqlLeWmNbgD+ncS4Dzk7n4rHE53HgY6Bugv19gD8Bfd39AzM7BfibmXV395WxOvNI/D5OcvcV6cQjIlLa2LGhf1CdTLedi0iNUV4itN7dR6RzMjMblOb1mwDnAe2B8xPsvxL4m7t/AODu081sFfD/gKtiZWm3EImIpKJ376gjEJGqVt7fOf0qcL60jnH3Be7+v3KqDADih27NBY5PNzARkYp4772wGv3mzVFHIiJVIWki5O4fJyo3s/pmdoSZnR573djM6pV3TEWYWXOgGfB53K6VwAFpnutkYDDQzcx+Vk69UWY2z8zmrVmzJt2QRSQLrVoFkybB9OlRRyIiVSGtJ99mdgGwDHidMKEiwBHACjMbmdHIoHHsa3wn5y1AXjoncvfn3f1Mdz/M3e8sp95Edy9w94KWLVumGa6IZKNjjoH994cpU6KORESqQsqJkJmdCjwILAbuA7bFdr0G/BS4zszOyGBsm2JfG8aVNwSKMngdEZGk6taFc8+FF14IM06LSHZJp0XoMuBKd/+Ou/8/YCuAu29z9yeBM4FLMxWYu68DvgLaxO1qAyzM1HVERPZk2DDYtg0OOiiMIMvPVwuRSLYob9RYvEOA/sl2uvsbZta+0hGVNROIHxVWADyT4euIiCT1wQchASpuEVqyJHSgBigsjC4uEam8dFqEyl1tx8waAXtVLpzd3AgMMrNDYtc4CdgfuDvD1xERSWrMGNi5s2xZUVEoF5HaLZ0WofeAXwM3Jdk/FngnnYubWQPgZWCfWNETZrbC3c8AcPc3zawQeMTMNhMmXRxUPJmiiEh1WLo0vXIRqT3SSYTGAy/GhqLPAPYys0sIkyGeRhjSntb8Pu6+lXIet8XqTAWmpnNeEZFM6tgxPA5LVC4itVvKj8bc/WXCLNCHElp/mgG3A78BWgBD3X1W5kOsXmY22Mwmrl+/PupQRKSGGD8e8uIm7WjYMJSLSO1m7p7eAWaNgUFAl1jRx8AMd9+U/Kjap6CgwOfNi5/UWkRy1ZQpoU/Q0qVhSP1ee8FHH0HbtlFHJiKpMLM3Ey3LlXYilCuUCIlIMh98AH37wre/Da+8EhZmFZGaLVkilE4fIRERAbp1g0cegZYtlQSJ1HZKhEREKuCMUvPor1kTkiIRqX3SWmtMRETKmjQpzDj9ccaWnBaR6qRESESkEgYNgnr14KyzwiSLIlK7KBESEamEjh3DiLIFC+AnPwGNPxGpXZQIxdE8QiKSrkGD4JprQgfqSZOijkZE0pFWImRmrc3sdjP7yMxWxcr6mtkdZtauakKsXu4+zd1HNWvWLOpQRKQWueoqGDJEo8hEapuUR42ZWUfg30BroAjYHtu1CjgImGtm33P3hRmPUkSkhqtTB559VomQSG2TTovQOOAToLu7NyEkQ7j7Enc/CZhAWHpDRCQnFSdBTz4J55yz+4r1IlLzpJMIDQB+6O4fJtl/C3BU5UMSEand1qyBP/8Zbrgh6khEZE/SSYQauPvnyXa6+zagceVDEhGp3S6+GM49F66+GmbOjDoaESlPOonQl2b27WQ7zew4YF3lQxIRqd3MYOJE6No1JETLlkUdkYgkk04i9BjwVzP7kZl1ADCzvc3sEDMbAzwFPFwVQYqI1DZNmsBf/gLffANPPx11NCKSTDprjV0PfAd4ACieMuyr2FcDphH6CYmICHDIIfD++2HSRRGpmVJuEYr1AToZOB94AfhvbJsGFLr7EHffUSVRViNNqCgimVScBL39NkyfHm0sIrI7c80Hn1BBQYHPmzcv6jBEJEsceyzMnQvz5oW+QyJSvczsTXcviC/P6BIbZqY+QiIiCTzyCDRqBGeeCZs2RR2NiBRLp48QAGbWBmgP5CXYfXKlIxIRyULt28Of/hTWJbvoopAYaRZqkeils8RGZ+BR4MhkVdjViVpEROIcfzyMHRvmFxoyBM46K+qIRCSdFqEHgDaEpTaWAVvj9hvwhwzFJSKSlcaMgf33D4mQiEQvnUSoD3Cwu69OVsHMTqt0RCIiWaxOHbjwwvD9mjXh9X77RRuTSC5Lp7P0J8AXe6hzeiViERHJGVu3wpFHwnnnaXFWkSilkwj9H3ClmZV3zPJKxiMikhMaNIBLL4UXX4Tx46OORiR3pfxozN1nmNn+wBIzewdYC8T/HdMsk8GJiGSz0aPh9dfhmmvgiCN2JUSzZkUalkhOSWfU2HDgQUKn6HZJqtX6UWNmNhgY3KVLl6hDEZEsZwb33RdmnT7jjPC4bOtWyM8PSVFhYdQRimS/dB6NXQ3cCXQGGrp7nfgNSNqRurZw92nuPqpZMzVuiUjVa9wYLrgANm8OSRDAkiUwahRMmRJpaCI5IZ1EqKm7/9Ldl8TWHUvkV5kISkQkl0yYADviVmosKgpD7UWkaqWTCP0nNqt0eWr9ozERkeq2dGl65SKSOekkQqOBW8zsRDNrkaTObRmISUQkpxSvUJ9quYhkTjqJ0BJgKDAdWGVmO+I3oHWVRCkiksXGj4e8uNUbzeCqq6KJRySXpDOz9Bbgz+XsN+DsyoUjIpJ7ikeHjRkTHoe1agWrVsHLL8OPfqTFWUWqUjqJ0Hp3H1FeBTMbVMl4RERyUmFh2eHyN90EDzwAa9dCy5bRxSWS7dJ5NNYvhTrdKhqIiIjsctll8NZbSoJEqlrKiZC7f5xCNU0ULyKSAWaw997wzTchKVq1KuqIRLJTOjNLpzJ+4Uzg4oqHIyIipS1cCHfdBfPmhT5D9dLp0CAie5S0RcjMXjWzd8ysbqxoMfDpHrZkw+pFRKQCuncPy3C8+qpGkYlUhfL+tjDKJkqbgFv3UF8zS4uIZNjw4fCvf8GNN4bFWYcMiToikexh7qlNBm1mn7v7/pWtU9OVWnT1x5988knU4YiIAKGv0He/C6tXw//+Bw0aRB2RSO1iZm+6e0F8eTqjxnqUc/K8PdWpLbToqojURI0awdNPw4wZSoJEMimdROh35ey7yczWAH0rGY+IiCSRnw9du4J7SIhSbNAXkXKkkwidVs6+3xH6B/2+UtGIiMgePf00DBwYJlwUkcpJJxFKyt3XAI8BzTNxPhERSe6MM0Ii9NOfhmH1IlJx5c5IYWYPlnrZzMz+SBgdlug8XYEVGYxNREQSqFsXpkyBww+Hs84KM1A315+hIhWyp6m5Lij1vQPJ1horAj5EkymKiFSLFi3CI7LvfhfOPx+mTdPirCIVUW4i5O4lj86yYWi8iEg26dsX7r8fWrdWEiRSUelM1n5llUUhIiIVMqJUO/3GjWF9MhFJXTqLrk6uwjhERKQSJk+Ggw6CpUujjkSkdsnIqDEREYnWUUdBURGcfTZs2RJ1NCK1hxIhEZEscPDBoVXoP/+BX2nVR5GUKRESEckSZ5wBl14K99wThteLyJ6l01laRERquBtuCPMKrV0bdSQitYMSIRGRLFKvHrz8cph0UUT2TI/G4pjZYDObuH79+qhDERGpkOIk6MUX4cILdy3O2r9/2ERkFyVCcdx9mruPatasWdShiIhUygcfwB//CLfdFnUkIjWXHo2JiGSpX/0K/vUvuPzyMAu1iOxOLUIiIlnKDB56CA44AE45Bd54A2bPhvx8jSoTKaZESEQkizVtCiNHwoYNsHVrKFuyBEaNUjIkAkqERESy3r337l5WVARjxlR/LCI1jRIhEZEsl2z9Ma1LJqJESEQk63XsmLi8Q4fqjUOkJlIiJCKS5caPh7y83ctbt4Zt26o/HpGaRImQiEiWKyyEiROhYcPwulOnUDZ3LvzgB7s6UYvkIs0jJCKSAwoLYdKk8P2sWeHrEUfA00+HVqEGDSILTSRSahESEclRP/0p/P3v0LgxbNwImzdHHZFI9VOLkIhIjihuCSqtbl3YuRNOPhkaNYLnnkvcn0gkW6lFSEQkx9WpEyZdnDkTBg+GTZuijkik+igREhERhg+HRx4JrUYnnQRffx11RCLVQ4mQiIgAMGxYWHbj9ddDC5FILlAfIRERKXHOOaGvUPfuUUciUj3UIiQiImWcdhocdBC4w223wVdfRR2RSNVRIiQiIgktWABXXAEDBsC6dVFHI1I1lAiJiEhCPXrAs8/Ce+/BccfB2rVRRySSeUqE4pjZYDObuH79+qhDERGJ3Mknw9Sp8NFHcOyxsHp11BGJZJYSoTjuPs3dRzVr1izqUEREaoRBg2DaNFi6FN55J+poRDJLo8ZERGSPBgyAxYthn33C661btT6ZZAe1CImISEqKk6Ann4SePWH58rL7+/cPm0htokRIRETS0r59SIL69QuPy0RqMyVCIiKSlqOOghkzwiiyfv3CIzOR2kqJkIiIpO073wmLtH71VUiGJk6EOXNg9mzIzw9LdYjUBkqERESkQgoK4O9/h7594Ze/hC1bQvmSJTBqlJIhqR2UCImISIUdfjjMnQtFRWXLi4pgzJhoYhJJhxIhERGplGQdptWRWmoDJUIiIlIpHTsmLm/ZsnrjEKkIJUIiIlIp48dDXl7ZMrOwHMfMmdHEJJIqzSwtIiKVUlgYvo4cGTpMd+oE114bZp8+9tiw75NPoEuXkCCJ1CRKhEREpNIKC2HSpPD9rFll961ZE0aW9e4N998PBx1U7eGJJKVHYyIiUqX22w9uvBHeegt69AiP0rZujToqkUCJkIiIVKk6dWD0aPjwQzj1VPjtb6FPH9i4MerIRPRoTEREMiT+kVi8/fcPC7ZOnx7q7r13KN++Heql8NuoeEHXPV1HJB1qERIRkWp1yilw663h+wULQifqZ54B92jjktykREhERCKzcyfssw+ceSacdhp89lnUEUmuUSIkIiKROeywsETHzTeHFe27dYN77ok6KsklSoRERCRS9evDb34D778PRx+tViGpXuosLSIiNULnzvDii7BjR3g9Y0ZY3f7qq+HZZ2HOnDBhY35+GIJfPJGjSGUoERIRkRrDbNcIsn/8A266CR58EDZsCEkQwJIlMGpU+F7JkFSWHo2JiEiN9LvfwSuvwJdf7kqCihUVwZgx0cQl2UWJkIiI1FjHHLPrUVm8pUurNxbJTkqERESkRuvYMb1ykXQoERIRkRpt/HjIyytbZgY/+EE08Uh2USIkIiI1WmEhTJwIDRuG1+3aQYcO8Ic/wOTJkYYmGdS//65lVKqTEqE4ZjbYzCauX78+6lBERCSmsBCOOAL69YNly+C998IvzREjQouRlueQilIiFMfdp7n7qGbNmkUdioiIJNG0KTz/PAwbFlazf+ihqCOS2kqJkIiI1EoNGsAjj8D999eO+YSievRTG0yZEibMnD07TJg5ZUr1XVuJkIiI1FpmYXLFhg1h3ToYPhzWrIk6KknHlCnh3zB+wszqSoaUCImISFZ45x148smwXtmiRVFHU/2qssVpyhRo1CgknplosXGHTZtC0nPFFWGCzNKqc8JMLbEhIiJZ4ZhjwtpkgwfDkUfCCy9Anz5RRxUUP/qpjWulJWuxgd3fwzvvwOrV8MUXsHZt2Lp1C1Md7NwJBQWhxW7tWvjmm/KvW10TZpqrq31CBQUFPm/evKjDEBGRNH30EZxwQvhl+/zzYaRZlIoTidKtHnl5YUqATCVDU6bAyJEhWenUqXKJ1ubN4bMrTmaGDYNVq3avV78+tGkTWuAefzyUtW4dEqHSRowI68UBnHNOeO/77QctWoTtt7+FlSt3P3+nTrB4ccXeQyJm9qa7F8SXq0VIRESyyre+BW+8ARdfDAcfHHU04RFPokc/l14K3buH/k1duoTEYvNm2L49dARv0CA8itqT8lpsTj8d1q+H/fcPr195BT74oGyi07Rp6HAO4dHa7Nmpva9t2+DYY+Hb395V9thj4RFacZKz7767FtEFeOKJ3c/TqFHiRHH8+NTiqCy1CCWhFiERkeywfTs89VRojUglscik5cuhffs911u2LEwUee21MHbsrvLihGj58pCw3HhjmESyYcNQ3rAhzJu3+6K0EN6re0hIijuQn3kmPPNM+H7ffcO+bt3guedC2R//GOoWt9jstx+cey6sWLH7+TPZYpPJFq1k1CIkIiI56bHHwuOZV1+Fe+4p20JRFbZtg+nTQ1Lx4ovJ67VsGVpitmwJCQfAoEGw996hbMsW2Lo1fG3UKOxv3x569iy7L1ESBCEJuuEGaNVqV9ndd8O990Lz5ok/h5Ejdy+7+eaqb7EpLIRJk8L3s2Zl7rypUItQEmoREhHJDu7h8dQNN8App4THM40bV821Nm4Mj+NWrgyPoy64ILSsXHVV1fURys8Pj8Pi1bYWm+IRb1WVCKlFSEREcpIZXH99aE255JLQr2X69NAiU1lffx0euy1cCOPGhdac0aPD6KgTTtjV6tK6ddUlEuPHV0+LTW0Z5ZYutQgloRYhEZHs89e/woUXwksvVXxovTv8+9/h0dcTT4RkqHt3eOut0G8nmaps8aiOFpvaTi1CIiKS84YMgeOOgyZNwusVK6Bt2/TOcffdoWUpLy/MjzNyZBhCXt0dsUuLso9NbaeZpUVEJKcUJ0EPPxz687zwQvK6O3aE1qOzz941sur000P/ns8/D4u9fve70SZBUjlqERIRkZw0aFBIhE49NYwqmzEjzGbcsSP84hdh7bLJk+Gzz0KH50GDwnHt2sGPf5z+9aq6pUYtQRWjPkJJqI+QiEj227gRjjoKFiwoW148B8+gQaFP0amnlt//R2o+9RESERGJs/fesGHD7uXuoeXnpZeqPyapXuojJCIiOe2zzxKXJ5pNWbKPEiEREclpHTumVy7ZRYmQiIjktPHjw1D40qpz0U+JlhIhERHJaYWFYTh8p06hk3SnTplb/kJqPnWWFhGRnJfNS0hI+dQiJCIiIjlLiZCIiIjkLCVCIiIikrOUCImIiEjOUiIkIiIiOUuJkIiIiOQsJUIiIiKSs5QIiYiISM5SIiQiIiI5S4mQiIiI5CwlQiIiIpKzlAiJiIhIzjJ3jzqGGsnM1gOf7KFaM2B9Bfe3ANZWILSo7ek918TrVOZc6R6bav1U6uXa/VVd91amr5WN91e23VuQe/dXTb239rS/Ku+vTu7ecrdSd9eWYAMmVrZOefuBeVG/x6r6XGradSpzrnSPTbW+7q+q/Tevzmtl4/2VbfdWpv/Nq/NaFT1XTb239rQ/ivtLj8aSm5aBOqmco7aprveUyetU5lzpHptqfd1fu6vO96P7q/w62XZvQe7dXzX13krnWtVCj8YiYmbz3L0g6jgkO+n+kqqie0uqUhT3l1qEojMx6gAkq+n+kqqie0uqUrXfX2oREhERkZylFiERERHJWUqEKsDMTjWzuWb2mpm9bmblPs80s6ZmNjl2zFtmdpOZ1Yurs7+Z/dXM3ojVuTTJuYaZ2ZdmdkEG35LUEFVxb8XqDTKzZWZ2bZLzfMvMXjGzf5jZm2Z2XobektQg6d5fsWO+bWYfmNnkJPv3+LNL91duMLMGZnajmW03s/wU6mfkd2Ol76/qHqZW2zegD/A10C32+hTgC6BNOcc8Azwa+74B8C/g+lL76wDzgOtir5sBC4FRperUB54G7gUcuCDqz0Jbzb+3YuV3xe6dVcC1Cc7RBFgKjIy9bg+sAwZF/Zloi/z+uhyYCbwPTE6wP5WfXbq/cmAD8oE3gIdjv6PyUzgmE78bK31/Rf7h1bYt9gvlL3FlHxT/QyWof2jspuhRquwHQBHQJPb6FGA7sE+pOpfF/nGL+3HlAQNj3ysRysKtKu6tWNnpsa+LkyRCPwXWAHVKld0D/CPqz0RbdPdXbP+pgAGzkiRCqfzs0v2VA1vs51EXoH8qiVAGfzdW+v7So7H0DSBkqKXNBY4vp/43wIK4+nsB3y1VZ6G7fxVXpwPQFcDdi9z95UpFLjVdVdxbuPuzKVz3LXffGXeeo8wsL4W4pXZI9/7C3ad67DdLOecs92cXur9ygrsvcPf/pXFIRn43koH7S4lQGsysOaFp7vO4XSuBA5IcdgCwKu6HycpS+4q/Jjpn6TqSxarw3kpFsvuvDqG5W2q5Ct5fqUjlZ5fuL0kkU78bK31/KRFKT+PY1y1x5VsIj66SHZOoPqWOSaWOZLequrdSvbbuv+xWkfsr1fPq55tURKbunUrfX0qE0rMp9rVhXHlDwnPNZMckqk+pY1KpI9mtqu6tVK+t+y+7VeT+SvW8+vkmFZGpe6fS95cSoTS4+zrgK6BN3K42hJ7siSwCWpmZxdWn1DGLkpyzdB3JYlV4b6Ui2f23k9DBWmq5Ct5fqUjlZ5fuL0kkU78bK31/KRFK30wgfu6Nglh5IjMInb+6x9XfDLxeqk4XM9snrs5n7v7fygYstUZV3FupmAEcbmalfx4UAP9yd/3Fnj3Svb9SkcrPLt1fkkimfjdW/v6KeshdbdsIc3FsBA6JvT6JMGdBm9jrcYRe8I1KHfMM8HDs+/rAP0k8V8LY2OumwCeUmishLgYNn8/CrSrurbjzL6b8eYRGxF63I8wvo3lesmiryP1V6thZlD+PUNKfXbq/cmsjyfD5qvrdmIn7K/IPrTZuhLk15gKvETLXb5fadyuhqW6vUmVNCZNMzQXeAm4G6sWdsy0wlTAh1dvAbxJc957YDyQHPop93zHqz0Nbjb+3ro7dK98QkqFZQN+4Ot8CXgX+ETvP+VF/FtpqxP01Kna/fEUYiTMLODXunKn87NL9leUbYULEWcD82O+oOcAze7i/MvW7sVL3lxZdFRERkZylPkIiIiKSs5QIiYiISM5SIiQiIiI5S4mQiIiI5CwlQiIiIpKzlAiJiIhIzlIiJCIiIjlLiZCIiIjkLCVCIiIRMDM3s8lRxyGS65QIiYiISM5SIiQiIiI5S4mQSBYxswFmttLMtprZYjPrbmYzzWxV6UcxZra/mf0xVnedmf3PzG4ws7wE5+xjZrPNbKOZLTWzB83sjNj51pvZ++XEMyJ2jR1mNit23NtmttrMVpjZODOrF3fMAWZ2u5l9WCq+F8ysdwXe5x7PleR8PWLv+Usze9/MTo7VKzSzBWb2lZk9a2YtE7zncj/b4s8kVv2HsXorzey+NM+T0meQLgt+ZGavm9kXZvaNmS0xs+lmVr8i5xSp0aJesVabNm2Z3wirQK8lrNrcHjDgMWAy0ApYQlh9vF2s/hHACsIKznVKnacrsAF4BdiP8MfT2bHjHbggxXgWA2uAvwEtY2VnAduBiXF1fxqLpUfs9T7A48BG4IBU32e654o736NAM6A+8DTwDTCcsBp7HeBgYB3wp7jj0/lsvTjOBHGkc55yP4MK3Dv3x2J7GvhJ7D2Pi3+v2rRlyxZ5ANq0acv8Fvvl6EDfUmU9gIHAA7F9B8cdc3Gs/IxSZU/FynrG1S0+xwUpxrM4lvR0iit/JnaeHqXKfgCMjqvXJHb8ram+z3TPFXe+Q0uV9Y+V/Tuu7iOxBMkSfC6pfLblJULpnKfczyDN+6YZsAO4P+p7WJu26tr0aEwke33j7v8pfuHu7wEzCS0xi93947j6c2NfTwAws7rAicBKd38nru6rFYhnkbsviSt7IfZ1cKk4n3T3+0tXcvevCa0h3ROcd7f36e4vV/BcAEXuvqDU6xWxr2/G1VsGNARaAphZHVL8bMtTwfMk/QzStI3QAtjHzPqaWSsza1aB84jUGvX2XEVEaqk1CcpaEv7qzyvVT6WYAZtidYrrNgY+SXCe+GNTsaqc83QuCcKsIeGRzLlAB3b1ZWwJLEpwjkTvs6LnAvgi7vXW2Ne1ScqL+1Wl89mWpyLnSfoZpMPdi8zsVODPwL9jxU8RWtZEspISIZHstbOcfR+4e6/qCiRNU4AhwHnAM+6+FcDMFiepX977TPdc5Z2vvOuUlqnPNp3zpBpbuczsTOAh4DlgOqEf1OJMnFukplIiJJJb1gBfAe0S7YyNptoUeySzhtD60CZB1URle9K6nPN8Grv+vsAZwEvu/kQFrlEik+dKUTqfbXWcJy2xz+thQr+ln2by3CI1mfoIieQQd99JeNTRwsy+V3qfme1F6EPUI1Z3B6EPTxsz6xl3qv4VuPwBZtYxruzk2Ndpsa9bCR1/PS62BiROpMqTyXPtUTqfbUwRsT9GY31xJphZkwqcJ1N6EB6FZjTBEqnplAiJ5J4xhMcdd5lZFwAza0EYMr6AMAy72FXA18BtZrafmdWJPT4pqMB1lwN3xq5V/BhmMDAp1pEbd98EvAgMNLPTYvXygLuARulcLJPnSkM6n+2HQNdY5+hBhH5Mmypwnj0ys3+a2fo9VPs4dv3rY3Mv/djMfmpmfzCzx9O5nkitEvWwNW3atGVuA3oSOiBvJQyDXgk8mqBea+A+QnLyOaFD9E1AswR1+wCzCXPvfArcSWjJceD8FONaTBjmfQyhE+4qwmiscUC9uLrNCXPZLAdWA28BvwKWxt7XyjTeZyrn6p7kfN8DLic8qnJCQvhW7LxvxV57bP/16X62wFHAO7Hj/wsMTuffKNXPIFZ3CXHD/5PUO5ow19M6whQDa4B/AMOjvre1aauqzdzLtBqLiOyRmf0QeAI4xd2fT6H+YsJw8P5VHJrEMbPvAHMI8ylNjDoekZpGj8ZEJKnYUhOTEuw6ktASMaeaQ5I0xJYAeQh4Fkj07yiS85QIiUh5GgMjzOyHsTWo6sb621wI3Ozu8XPuSM2ykdAn6mxX879IQno0JiJJmVlzQj+Zk4AWhMkDlwD3uvu9KRw/AriBMAHgDkLfk7Pd/R9VFrSISBqUCImIiEjO0qMxERERyVlKhERERCRnKRESERGRnKVESERERHKWEiERERHJWUqEREREJGcpERIREZGcpURIREREctb/BzqFADG3RiIKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (9, 6))\n",
    "\n",
    "plt.errorbar(regs, tm1, yerr = tm1_err, fmt = 'bo--', label = 'stabilized sinkhorn')\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.yticks(fontsize = 14)\n",
    "plt.xticks([1e-3, 1e-2, 1e-1, 1e0], fontsize = 14)\n",
    "plt.gca().get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())\n",
    "\n",
    "plt.xlabel(r\"reg parameter, $\\varepsilon$\", fontsize = 18)\n",
    "plt.ylabel(\"time elapsed (s)\", fontsize = 18)\n",
    "\n",
    "plt.legend(fontsize = 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measure the dependence of estimated time on the size of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.15 ms ± 142 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "18.7 ms ± 468 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "18.1 ms ± 198 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "5.72 ms ± 23.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "118 ms ± 17.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "169 ms ± 20.4 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "249 ms ± 46.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "ns = 2 ** np.arange(3, 10)\n",
    "\n",
    "tm1, tm1_err = [], []\n",
    "\n",
    "for n_samples in ns:\n",
    "    x_i = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "    x_f = empirical_distributions.generate_1d_normal_numpy(n_samples)\n",
    "\n",
    "    C = cdist(x_i, x_f)\n",
    "    C /= C.max()\n",
    "\n",
    "    a = np.ones(n_samples) / n_samples\n",
    "    b = np.ones(n_samples) / n_samples\n",
    "    \n",
    "    t = %timeit -o sinkhorn.sinkhorn_knopp(C, reg = 1e-2, a = a, b = b)\n",
    "    tm1.append(t.average)\n",
    "    tm1_err.append(t.stdev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9c54330e50>"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAF/CAYAAABKac3QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABKCklEQVR4nO3dd5hU5dnH8e8N0puhSG+KEEAEZcVooqAiKgq2WF4QjQ1sqLEb7II1dhCDRmxEY+8SG4hJBAVFRaKiIihIkybVhb3fP57ZZVl2l53d2TlTfp/rmmv2POfMOffMHph7n2rujoiIiEi2qhJ1ACIiIiJRUjIkIiIiWU3JkIiIiGQ1JUMiIiKS1ZQMiYiISFZTMiQiIiJZbYeoA0hVjRs39nbt2kUdhoiIiCTAjBkzlrl7k+L2KRkqQbt27Zg+fXrUYYiIiEgCmNm8kvapmUxERESympIhERERyWpKhkRERCSrKRkSERGRrKZkSERERLKakiERERHJahpaXwGrV69myZIl5ObmRh2KSFqqVq0aO+20E/Xr1486FBHJYkqGymn16tUsXryYli1bUqtWLcws6pBE0oq7s379ehYsWACghEhEIqNmsnJasmQJLVu2pHbt2kqERMrBzKhduzYtW7ZkyZIlUYcjIllMyVA55ebmUqtWrajDEEl7tWrVUlOziERKyVAFqEZIpOL070hEoqZkKIX06RMeIiIikjxKhoowswFmNm7VqlVRhyIiIpI1oqwQUDJUhLu/4u5DGzRoEHUoGWP69On06NGD6tWr86c//Snu1/fv359mzZpF2pwyfvx4unTpgpnxyCOPRBbHyy+/TI8ePTAzrrvuuhKPe+CBB1IiXhGRdKBkSMrsuuuuY/LkyXG/Licnh5kzZ9KiRYsyHb/nnnty8cUXF2y//vrrnHXWWdscd/fdd9OuXTvWrVsXd0zxOvXUU3n99dcTdr5jjjmGY489Nu7XDRw4kJkzZ273uLPOOiuh8YqIZDLNMyRldv311wPQp5LrMdu0acNOO+203eMaNmxImzZtqFq1aqXGUxlatGihjsMiIilCyVCKmDABpk6FjRuhXTsYNQoGD446qmi8+OKLZTru5JNP5uSTT67cYCrJ6NGjow5BRERi1EyWAiZMgKFDQyIEMG9e2J4wIXkx/O9//6N///706NGDHj160Lt3bx566CEApk2bRo8ePYDQFyX/mNmzZwMwbtw4fv/735OTk0O3bt04/PDD+fLLL0u81i233EKvXr1o0qQJffv25ZtvvinYt/fee9OwYUPatWtXaryjRo2iQ4cOmFlB0933339P3bp1C+Lr0aMHHTt2xMy2apJavnw5Z555Jm3btqVjx4706tWLN954Y5tr3HrrrbRs2ZLOnTvTv3//reIsTW5uLpdffjm77bYbe+yxBz169OCiiy5i6dKlABx99NHb9IG65pprCt7PM888w4knnkiXLl3Yddddee6550q93ldffcVuu+2GmbH77rtz++23b7V/7dq1DB06lO7du9OuXbtiE7FZs2ZxxBFH0K5dO9q3b0+/fv34+OOPC/YX7oN0//33M2zYMHr27EnVqlW58MILOeOMM2jTpg1mxjvvvMPAgQPZdddd6datG1OmTCnT5yYiEhl316OYR8+ePb00s2fPLra8d+9tH2PGhH1r1xa/v1Ejd9j2UaNG2P/UU+H18+cX//qXXw77v/yy1JBLtcsuu/iNN95YsD169Ghv27btVscAfu21127z2t/+9rc+ceLEgu077rjDW7Zs6atXr97quLZt23qTJk38wQcfdHf3DRs2eL9+/bx9+/a+cePGguNOOeWUba597bXXerhdt5g0aZIDPmnSJHd3nzt3rvfu3XurY0466SSvWbNmwe9rw4YNvscee/jee+/tq1atcnf3p59+2qtWrervvvtuwevGjh3rVatW9TfffNPd3ZcuXeoHHnigAz5+/PhtPoPCbrzxRu/WrZuvW7fO3d2/++47b9KkSUGc23s/+++/vy9btszd3a+88kqvU6eOL1++fKtjC/8uNm3a5CeccILffPPNWx0zd+5cB7xbt27+3Xffubv73/72Nzcz/7LQzTJnzhyvX7++X3zxxZ6Xl1cQX506dfyzzz7b5nydOnXyzz//3N3D7/qCCy5wd/fx48c74H/84x8L3vugQYO8devWnpubW+pnVtK/JxHJHvnfaZUFmO4lfOerZigF/Pxz8eX5NUWVbdmyZXz77bfsvPPOBWVnnHEGp512Wple/8ILL3DIIYcUbJ9zzjksWLCg2A68DRs25IwzzgCgRo0ajBw5krlz5/LYY49V8F1A48aNufTSSwu2n3nmGZ544gluueUWOnfuDMDjjz/OJ598wsiRIwvWwjruuOPIyckp6BOVl5fHyJEjOfjggzn44IMLzl3Wz2Pq1Kk0a9asYIby9u3bc+utt9KqVasyvf7oo4+mUaNGQOhovXbtWj766KNij83Ly+OUU06hS5cuXHHFFcUec+CBB9K+ffuC87k77733XsH+/FFpN954Y0Ft1YgRI6hTpw4jRozY5nwHHXQQu+22GwDDhg3b5rqDBw8ueO9HH300P/zwA999912Z3ruISBTUZyjBShtsVbt28fvbtQtNY0W1bbv18a1bl37+Tp3KFOI2GjVqRNeuXRk2bBgzZ85k0KBB9OjRg2uuuaZMr8/Ly2PQoEHMmjWLKlW25NfFfQF26dJlq+099tgDM+ODDz4oSJLKq27duhx++OEALFy4kLPOOou+ffty/vnnFxzz9ttvY2bsu+++W712t91247HHHiM3N5eFCxeyYMECTjnllG2OKYv999+fyy+/nMMPP5yhQ4fSr18/Tj311DK/j44dOxb83LBhQwAWL168zXHuzumnn87UqVN5/PHHy32+t99+m65du261vEy1atXYY489ePvtt3H3rZr08hNLgDp16lCnTp0yXa9wuYhIYVH3m1XNUAoYNSokSoXVrh3KkyG/382wYcN45JFH2GOPPejSpQsvvfTSdl+7YMEC9ttvP/Ly8pg6dSozZ84sGPq9sZiqraIrk++www7Uq1ePhQsXJuS9QEgSTj31VNydRx55ZKsv8mXLlhUkQ4X7Fk2aNImGDRuyYsUKFi1aBMCOO+641XnLOvfUZZddxqOPPsqSJUs46qijaNq0KZdcckmxn0dxahe6GfKTy82bN29z3OjRo/n555/59ttvGTNmTLnPt2zZMn7zm99s87qGDRuyfv36baYuqFu3bkLiFxGB1Og3q2QoBQweDOPGQY0aYbtt27CdzKy4cePG/PWvf2XhwoW8+OKLVK1alWOPPbbUjtAAr732GsuXL+eyyy7b6kuwJKtXr95qe9OmTfzyyy9lnoOoLEaPHs2bb77J2LFjadmy5Vb7GjduTJUqVfj4448LEreZM2fy7bffsmjRInbaaSeaN28OwIoVK7Z67cqVK8scw8knn8xHH33EF198wfHHH88dd9zBjTfeWOH3Vli/fv14+eWXGTJkCFdeeSXziqteLIPGjRuzfPnybcqXL19OrVq1yvR7FREpj59/hgsugKLTxa1bB8W00lcaJUMpYvBg+N3voHdv+P775CZCS5Ys4YILLgBCTc2RRx7Jo48+yubNmwtGjOXvC33QwuizmTNnFtR2FK59ya9ZKc4XX3yx1fbHH3+Mu7PPPvsk5L18+eWXXH755QwaNIgTTjihoPymm24C4OCDD2bTpk3bxPHJJ58wbNgwAFq1akWrVq22Gk1VXOwlufLKK/n++++B0Cz40EMP0a1bNz777LPyvq1idYq1i959993UqVOHoUOHlus8ffv25YsvvtiqBmjTpk3MnDmTvn37aj4kEUm4V1+F/feHnXYqud/s/PnJi0fJkLBu3ToeeOABpk6dWlD273//mzp16tCrV6+Csvbt2/Pjjz8CMHLkSF5++WUOPvhgqlWrxj333MPmzZtxd0aOHFnitX788UcefPBBIDSjXX311bRv3z4h8wXl5uZy0kkn0bhx422ajcaNGwfAkCFD6NmzJxdffDG//PILEGpAhg8fXpBcVKlShWuuuYa3336bt956C4Cff/6Z++67r0xxfPDBB9xzzz3k5eUBMG/ePH788UcOPPDACr/H4jRs2JAxY8bw5ptvlmvpjWuvvRYz4+qrry5Idm+66SZ++eUXRiWrrVZEMlZuLkyaBBddBPldSVetgtWr4S9/gWbNin9dmzbJizHyIeyp+ijv0PqKqOxhhSVZt26dX3PNNd6jRw/v0aOH77bbbt67d2+fMmXKVse9+OKLvvPOO3u3bt38gAMO8CVLlri7+wsvvOBdu3b1tm3bep8+fXz06NEOeNOmTX3IkCH+0Ucfeffu3b1atWo+ZMgQv/zyy33PPff0xo0b+0EHHeRz5swpuEavXr38N7/5jVerVs27d+/un3zyiR922GHetGlTB7x79+7+2muv+ciRI32XXXZxwHfZZRe//vrr/cknn3TA69ev723btt3mkW/FihU+bNgwb926te++++7es2dPH5M//0Eht912m7do0cI7derkffr08RdeeMEBb926tZ9++uklfp4vvvii9+vXz7t27erdu3f3rl27+s0331wwbP2oo47a5v3ceeedW72fhx9+2N966y3v3LlzwTUvv/xynzJlinfv3r3g8x00aJCvXr3au3fv7lWqVPE6dep49+7d/amnntrqtSNHjvQvvvhiq9cOGTKkIObPP//c+/fv723atPG2bdt63759ffr06QX7i56ve/fuvmnTpoL9l1xyibdu3doB79y5s0+cONEfffTRrd7TPffcU+JnpqH1Ipll7Vr3J55wP+EE9wYNwnQx1au7P/ts2B/779Ddw3G1a289tUzt2qE8kShlaL157C9B2VpOTo5Pnz69xP3/+9//thpVkwj5q1yUY/kvkbRWGf+eRCS5vvoK1qyBnj1D09dOO0HjxnDEEeFx8MFQ0viLCRPg9NNDJ+q2bStnNJmZzXD3nOL2aWh9ClESJCIi6SI3F/7zH3jllfCYMwcOOADefRcaNYLPPoPOnaFKGTrkDB4MsR4UkXwXKhkSERGRMlm7FvKnFjvySHjjDahePSRBF1wQaoDyde0aTYzloWRIRERESvT111tqf6ZNg0WLoEEDOP98OOOM0PxVr17UUVaMkiERERHZxuTJMGxYSIYAunULI8Jyc8P2oYdGFlrCKRkSERHJcitWwMSJofbnxBNh4EBo3hzat4fhw2HAgNCxOVMpGaoAL7Jmk4jETyNaRaKxaRPce29IgN5/HzZvDqO/8kc2d+oUEqRsoGSonKpVq8b69eu1VIFIBa1fv55q1apFHYZIxtu0KYz+WrAABg2CqlXh/vuhVi247LJQ+9OrVyjPNkqGymmnnXZiwYIFtGzZklq1aqmGSCRO7s769etZsGABTZs2jTockYy0cuWW5q833gjNYc2ahaawKlXg44+hyPrZWUnJUDnlr76+cOFCcvN7k4lIXKpVq0bTpk0L/j2JSMV98w20awc77ADXXw933x2avwYODLU//fptmfsnlf7pRTnXnmagLsLMBgADOnTocOacOXOiDkdERKRUmzbBf/8ban9efRW+/BKmTIH99gsTIS5dCnvvnZ3NX4WVNgO1kqESbG85DhERkai4gxn873/whz/A8uVQrVro/DxgABx/PKj1eWtajkNERCTNffPNlskP99knrN/VoQMccwwcckho/kqlZq90omRIRESkEiRq8e2RI8NCpl9+Gba7dt1S61Ot2pY1vaT8lAyJiIikiFWr4F//go8+gttvD2Vffw2tWsHZZ4e1v3beOdoYM5GSIRERkQj9+CM8+2zo/Pzee6FDdKNGcPnlYRTYo4+G/kFSeapEHYCIiEimmTABpk4NyU27dmE73+bN8O9/hwVPIRzz5z/DTz+Ftb/efx8WLw6JECgRSgbVDImIiCTQhAkwdChs3Bi2582DM88MydHKlfD662H01513hiRo4MDQOXqXXSINO6spGRIREUmgESNg3bqty9avh9GjoWFD6N8/DH8/5JCwr1698JDoKBkSERFJgDVrYNKkUBNUHLPQ/LWDvnlTjvoMiYiIlEPhOYsvuSTU+gwcWHIfnzZtlAilKiVDIiIiZbR8Ofzzn3DaaaFj9MqVobxbN7jwQnj7bRg/HmrX3vp1tWuHSRIlNSlHFRER2Y6pU0Nn5w8/hLw82HFHOPjgMC/QjjvCKadsffwOO8Dpp4dO1G3bhkRo8OAoIpeyUDIkIiJSyMKF8OabMHEinHgiHHUUNGgQmsWuugoOPRT22qv0Jq/Bg7fMDB3lauxSNkqGREQk6+XmhkRn4kT47LNQ1qwZ9O0bfu7cOdQOSWZSMiQiIlnnu+9C4rNuXej8XK1amAG6aVO49dYw7H333TXhYbZQMiQiIlnh/ffhmWdCEjRnTijba6+QDAF8+qlGe2UrjSYTEZGM4w6zZoVZnnNzQ9mLL8Lf/w4dO8J994UFUKdN2/KaRCdCkyerv1C6MC88UYIUyMnJ8enTp0cdhoiIlNHq1WHF94kTw/OCBaF82jTo1SsMi69dG2rWjDZOiYaZzXD3nOL2qUJQRETS0ubNMGNGWOF9l13gk0/g+OO3DHs/5JDwaNUqHN+wYaThSgpTMiQiImlj0aIttT9vvhlqey69FG67DfbdF/7zn1ALpL4/Eg/dLiIikrJ+/RV++CHU/OTlQZcusGJFGPV1xBFhzp+DDw7HVqsWEiKReCkZEhGRlDJ37pZ+P++8Ay1awFdfQZUqoQN0+/Zh2HsVDQGSBFEyJCIiZdKnT3hO9AipdeugVq0wp8/FF4cRYBDW/jrppFD74x72H310Yq8tAkqGREQkydxh9uwttT9TpoQ5fjp1gsMPD2t5HXoo7LqrJj2U5FAyJCIiSTN9eqjd+fHHsN21K5x7LtSoEbYPPDA8RJJJyZCIiCRcXl4Y9p5f+3PssWHV9513hr33hmuvDcPeW7eOOlIRJUMiIpJA7nDqqfDaa7BsWSjLydkyx0/DhvDss9HFJ1IcJUMiIrJdEyaEVds3bgwdm0eNguOOgw8+CLU/P/8M48aFPj7r1sFhh20Z9t6kSdTRi5ROy3GUQMtxiIgEEybA0KEhyclXtWqY2HDjxvC8//5hEsSqVaOLU6Q0pS3HoVkaRESkVCNGbJ0IQVgKo1o1eP750Bz2zjtKhCR9qZlMRERKNX9+8eVr12reH8kMqhkSEZFStWkTX7lIulEyJCIipbrhhm0nP6xdO3SiFskEWZMMmdnOZvaCmT0VdSwiIumkZs0wZL5atbDdtm0YOTZ4cLRxiSRKNvUZ2huYCBwQdSAiIulk9OiwOGrr1qGGKNFrk4lELSVqhsysupndYmabzKxdMfsHmtlHZjbFzP5jZsUOjSuNuz8JbExEvCIi2WLDBmjcGIYP1zphkrkirxmKJT9PAl8D2wzMNLOewD+AXu4+28yOAP5lZl3dfVHsmOkU/176u/vCSgteRCTD1awZhs8DvPRStLGIVJbIkyGgLjAEaAWcXMz+K4F/uftsAHd/1cwWA+cCV8fK4q4pEhGR0q1aBUuWhNXjRTJZ5M1k7j7L3b8p5ZC+QNGpoD8CDq68qERE5KGHoGNH+O67qCMRqVyRJ0OlMbOGQAPgpyK7FgE7x3muw4EBQBczOz8xEYqIZKbNm+H+++EPfwgrzUPoOK3O05KJUqGZrDR1Ys9FOz5vBGrHcyJ3fw14rbRjzGwoMBSgjWYTE5Es9sYboUbo5pujjkSk8qV0zRCwNvZco0h5DaDISjkV5+7j3D3H3XOaaJllEclio0dDixZabkOyQ0onQ+6+HFgJNCuyqxnwbdIDEhHJAkuXwnvvwVlnbZloUSSTpXozGcDbQNHRYjnA8xHEIiKS8Zo0gR9+gB3S4RtCJAFSumYo5hbgEDPrDGBm/YHmwJhIoxIRyUB5eeG5cWPYccdIQxFJmsjzfjOrDrwJ7BgresrMFrr7MQDuPsPMBgOPmdl6wsSMh+RPuCgiIonzt7/B44/D668rGZLsEXky5O6/An22c8zLwMvJiMfMBgADOnTokIzLiYikDPfQcbpWLWjQIOpoRJInHZrJksrdX3H3oQ30P4GIZJlJk2D2bDjvPK1DJtlFyZCIiAChVqhRIzjhhKgjEUkuJUMiIsK8eWEh1jPPDM1kItkk8j5DIiISvcaNYcwY6N8/6khEkk/JkIiIUKdOmGRRJBupmUxEJMu98gqMHQubNkUdiUg04k6GzKyqmTU2sxZmtmMlxBQpMxtgZuNWrVoVdSgiIpXOHa6/Hu67D6pWjToakWiUKRkys4PM7O9m9jWwHlgM/AD8bGZrzGyamd1oZrtUZrDJoKH1IpJNpk2DGTM0nF6yW6l9hmLJzSPA74HVwFfAbMJq8puA6sBvgFbAZcAVZnY/cIm751Ze2CIikgijR0O9ejBkSNSRiESnxGQothbYFOA94HfAR+7upRxfGzgS+Avwhpkd6u5qgRYRSVGLF8PTT8PZZ4eESCRbFdtMZmZVgSeBc939j+7+YWmJEIC7r3P3J4HuhGa06xMerYiIJMySJbDXXnDOOVFHIhKtkmqGqgHHuPt38Z7Q3fPMbAjQsUKRiYhIperWDf7zn6ijEIlesTVD7r6hPIlQodfnufuX5Q9LREQq0+zZsHRp1FGIpIYyT7poZvWBHsA6d58eK2sCnAfsCPzT3f9bCTGKiEiCnXUWLFsGX3yhUWQi8cwzdDYwGbgawMxqEjpYXwUMB94zsz6JDS/5NM+QiGS6Tz+F99+H005TIiQC8SVDRwLD3P3I2PYgoBNwDVCPkCSNSGx4yad5hkQk040ZExZjPe20qCMRSQ3xJENt3P3BQtuDgR+Bm9x9LXA70DmRwYmISGItXw5PPAGDB0PDhlFHI5IayrU2mZntBOwPPJE/5N7dNwOqcBURSWHvvAMbNoQZp0UkiCcZmm9mh8V+viL22n/k7zSzLoA62oiIpLDjjoN586B796gjEUkdZR5NRmgGe9nMVgKNgJfdfRaAmV0CnAO8kfAIRUQkITZvDouxtm4ddSQiqaXMyZC7v2Bm/YABwBLgPgAzM6ArYdmOv1VGkCIiUnFHHhkSobFjo45EJLXEUzOEu08CJhUpc+DURAYlIiKJ9fXX8NprcMMNUUciknpKWpusSmweoXIzs7oVeb2IiCTO/fdDtWpw5plRRyKSekrqQF0D+NjMfl+ek5rZSOCv5Y4qQpp0UUQyzZo1MH586DzdrFnU0YiknpLWJlsPDCN0mB5lZo3LcjIz29PMXiFM0HhJ4sJMHk26KCKZ5oknYPVqGD486khEUlOJfYbc/X0zOxJ4HLjMzD4HvgQWAWuBTYTV7X8DtAK6A62BfwG93X1NJccuIiJlcOSRkJcHe+8ddSQiqanUDtTu/u/Y/EF/Ak4EjgGqF3PocuBt4Ex3fyvRQYqISPk1bw7nnBN1FCKpa7ujyWJNZmOBsWZWFWhHmGeoOrAOWOjuiyozSBERKZ8bboCcHOjfP+pIRFJXXMtxuPtmd//W3T9093+7+8dKhEREUtP8+XD99TBlStSRiKS2cq1NJiIiqe+BB8Lz2WdHG4dIqlMyJCKSgTZsgAcfhIEDoW3bqKMRSW1KhkREMtDTT8OyZRpOL1IWSoZERDJQzZowYAAccEDUkYikPiVDIiIZ6Pjj4eWXwSzqSERSX0KToYquZ5YKtByHiKS7yZNh/fqooxBJH4muGfouwedLOi3HISLpbPFi6NcPrrkm6khE0keJky6aWXn+KWmlehGRCD34IOTmwhlnRB2JSPoobQbq64opc6C4FmhPSDQiIlJuubkwdmyoGerUKepoRNJHacnQMmCvQtuNgNHAa8AkwoKtAM2AA4BBwNBKiFFERMrgxRdh4UL429+ijkQkvZSWDP3T3eflb5jZDcCV7v5ekeO+A/5rZv8G/gz8J/FhiojI9rz9NrRvD4cdFnUkIunF3MvWwmVm37t7u+0cM9/d2yQisKjl5OT49OnTow5DRKTM3GHJEmjaNOpIRFKPmc1w95zi9sUzmqyhmbUs5SKtgPrxBiciIhX3669hTiElQiLxiycZmgK8bmZHmFm1/EIzq25mA4BXY8eIiEgSrVgBLVvC449HHYlIeiqtz1BRFwDvAy8BmNnKWPmOsedFwNGJCkxERMpm/PiwDtnuu0cdiUh6KnPNkLt/C3QD7gS+AmrHHl8BfwV2d/e5lRGkiIgUb/NmGDMG9tsPunePOhqR9BRPzRDu/jNwaewhIiIRmzgRvvsObr456khE0le5luMws2Zm1jP2sxZ7FRGJyOjR0KIFHK1OCiLlFlciY2Z9zWwmsAB4N1Z8gJl9amaHJDq4KGihVhFJJ3fcAQ89BNWqbf9YESlemZMhM/sD8DphJuo3gE2xXR8BE4CnzKxPguNLOi3UKiLppEsXTbIoUlHx1AxdDfwNaO/uRwAbAdx9tbvfBpwEXJX4EEVEpKg1a+CUU+CLL6KORCT9xZMM9QQudfdNxe1099eAXRMSlYiIlOrxx+Gxx2D16qgjEUl/8SRDRqw2qNidZjsQhtqLiEglcg8dp3v2hN/9LupoRNJfPMnQd8DgUvafA8ypWDgiIrI9kyfD7Nlw3nlhCQ4RqZh45hm6C3jMzPoBbwH5y3C0Ao4BDgSOT3yIIiJS2OjR0KgRnHhi1JGIZIYyJ0Pu/g8zawPcSKghMuDF2PMm4DJ3f64yghQRkcA9jCDr1Qtq1ow6GpHMYO4e3wvM2gLHAh1iRV8Dz7v7/ATHFqmcnByfPn161GGIiIhIApjZDHfPKW5fXMtxALj7PML6ZCIikkQbNsC778Khh0IVzf0vkjDxTLpY38z2N7OcQmVNzOx6M7vHzPatnBBFRATgn/+Eww+H99+POhKRzBLP3xZnA5MJky9iZjWBKbHt4cB7mTADtYhIKnKH++6Dzp1h//2jjkYks8STDB0JDHP3I2Pbg4BOhGSoXux5RGLDExERgA8/hBkzNJxepDLEkwy1cfcHC20PBn4EbnL3tcDtQOdEBiciIsHo0VC/Ppx8ctSRiGSecnXBM7OdgP2BJzw2HM3dNxOG2YuISALl5sK0afCnP0HdulFHI5J54hlNNt/MDnP3N4ArCInUP/J3mlkXYFWC4xMRyXrVqoUZp9etizoSkcwUTzJ0O/Cyma0EGgEvu/ssADO7hLAcxxsJj1BEJItt2gR5eVC9emgmE5HEK3Mzmbu/APQDHgf+QmydMjMzoCvwHvC3SogxqcxsgJmNW7VKlVwiEr0XXoA2bWCOVn4UqTRxz0CdLTQDtYikgt694YcfQjJUtWrU0Yikr4TOQG1mDQk1RO1jRXOBN919eflDFBGRoj77DKZMgdtvVyIkUpniSobM7ArgGqAGW0aOObDRzK5391sTHJ+ISNYaPRpq1YLTTos6EpHMVuZkyMwuAG4CPgFeBxbFdjUDjgBuMrP17n5vwqMUEckyK1fCE0/A4MHQsGHU0Yhktnhqhi4ALnD3+4rZd7WZDQcuBJQMiYhUUIMG8Prr0LJl1JGIZL54Jl2sX0IilG80oIGfIiIJYAZ9+sCuu0YdiUjmiycZWmRmtUvZXwdYULjAzNqUKyoRkSz21ltwwQWgGT5EkiOeZOg+YKyZ1Sq6I1Z2D3BjkV1TKxCbiEhWuvNOeOYZqF3an58ikjDx9Bn6HWFI/QIz+xBYHCtvCuwFLCHMwXh4odc0SEiUIiJZ4uuvYeJEuP76sAyHiFS+eJKhUwr93K+Y/b8BOhUp04yOIiJxuP/+kAQNHRp1JCLZI55kaLG7N4/n5Gb2U5zxiIhkrTVrYPx4OO44aNYs6mhEskc8ydD95Th/eV4jIpKV1q6FP/5RtUIiyZbQtcnMrJ+7v5mwE0ZIa5OJiIhkjtLWJotnNFlZPJLg84mIZIVZs0B/f4lEI961yX4HXAR0BbYZYg80SURQIiLZ5pprwqKsCxZAjRpRRyOSXcpcM2RmBwDvA78HlgEtgXmxRx7QDvgw8SGKiGS2+fPhpZfgzDOVCIlEIZ6aoWsI645d4u5uZj+5+wH5O83sGkCzYoiIxOmBB8LzWWdFG4dItoqnz1B34C++pcd10Z7Xo4DDERGRMtuwAR58EAYOhLZto45GJDvFkwytd/eNhbY3m1nhCt08QGuRiYjEYdYsyM2F4cOjjkQke8WTDC01s/0Kbc8FhhXaPgtYnZCoRESyRE5O6DR9wAHbP1ZEKkc8ydBrwGtmdm1sezxwt5nNMrPPgdHAU4kOUEQkU/3yC7hDnTpgFnU0Itkrng7UdwP/BVbGth8hDLE/HdhImG36hsSFJiKS2c4+G775Bj74QMmQSJTKnAy5+1JC7VD+tgOXxB4iIhKHxYvh6afDCDIlQiLRSugM1GZ2ZSLPFwUzG2Bm41atWhV1KCKSwR58MHScPvfcqCMRkUQvx5H24yHc/RV3H9qgQYOoQ8l6EyZAu3ZQpUp4njAh6ohEEiM3N8wt1K8fdOoUdTQiUmIzmZltTmYgIoVNmBBW7l63LmzPm7dlJe/Bg6OLSyQRXnkljCAbOzbqSEQESu8ztBH4ZxznMuC4ioUjEowYsSURyrduHZx3HixdGpYsaNwYjovdcdOmwerVULNm2FejBtSvD+3bh/2//AJVq4byqlWT+16SZcKE8LnNnw9t2sCoUUocU9URR8Bzz0H//lFHIiIAtmVC6SI7wnIbzeM6WTlek6pycnJ8upaQjkyVKmHIcWk6doSvvgo/9+4dFrksbI894OOPw885OTBjRvh5hx1CUtSnD7z6aig7/HBYuHBLIlWjBvzhD2HxTIBLLw3JWH6yVbMmdO8ORx4Z9j/xRIg5f1+NGiEh6dgx7P/6663PnX9clQQ1VBetSQOoXRvGjVNCJCICYGYz3D2nuH2l1Qz1Lse1yvMakW20bh1qOIor//RT2LgR8vK2lI8dC8uXh6UNNm4Mj7p1t+y/8MLQLJG/b8OGLbVGEPokVa26Zf8vv0DhPvT/+ldIlvLPv2kTnHTSlmTozDPDvsLOPhvuvz8cW1y/kEsugdtvDzVaO++8dZJUo0aoBTv9dFi2LDwX3X/ccSEJXLYsHFtcTdqIEUqGUs3VV8NvfgMXXRR1JCKSr8RkyN2/Lq7czKoBPYHm7v6CmdUBNrr7ppJeIxKvffbZNhmqXRtuvjl8kRTVpUvp5zvppNL3jxlT+v7PPtt6e/Pm8Mj35ZdbEqX856ZNwz6zUHNTeN/GjbDXXmF/lSpwwgnb7s9P5jZuDH2mCidyGzfC7ruHZGjBAli5svi4i0soJTorVsAddyhBFUk1JTaTFXuw2Z+AW4HGwBp3b2BmBwFPAle6+98rJcoIqJksOgsWhJqUTp3g55/VB2Z78vJCzdYPP2y7r21b+P77ZEckJbnzTrj4YvjkE+jRI+poRLJLac1kZe6xYGYDgYeB74EHgNzYrinAecCNZnZMxUIVCf1zNm2CZ54JX+R5eeFZiVDxqlQJNWa1a29dXrt2SCAlNWzeHGog//AHJUIiqSae7puXEWp/9nb3c4FfAdw9192fBo5Fs1FLBf36a+hDc+mloR+NlM3gwaGzdNu2oVmubduwvXkzXHfd9jujS+WbOBG++y707xKR1BLP2mSdgT4l7XT3D8ysVYUjkqxWvXoY4bVZs1zFbfDgbWvPhg0LSdHChaEz9w7x/IuXhGrSJPx+jj466khEpKh4/mssdfUcM6sJ1KpYOJLNXn8dfvvbUCOUqXMBJdsDD4Qv4VGjwlpYTz65bXOaJEevXmEKBhFJPfE0k30OXFzK/uuBTysWjmSrZcvCiK/hab+gS2oxg5EjQ1+VV16Bvn3DSDRJrldegblzo45CREoST83QKOANMzsceAuoZWbDgVbAUcDOwMEJj1CywogRYW6f22+POpLMdM45Yaj/7NlhjiJJnjVrQqJ/xBFaX08kVZU5GXL3N81sCDAa+EOs+G5C89kKYJC7T050gJL5pk8PK3hfeOH25wuS8jv22PAA+OijkBTtvnu0MWWDJ54IgwLUcVokdcU1zxBAbJLFQ4AOsaKvgbfcfW2CY4uU5hlKjrw8+P3vQxPCV19BgwZRR5T53GHPPcPIppdeCsuSSOVwh27dQuI5fXpothSRaJR3OY5ixZKe5ysclQhhKP1ee4VmHCVCyWEWkqBDD4VDDgk1F8dpieVKMXkyfPEFjB+vREgklWmgrUSqZk24996oo8g+bdrAv/8NAweGpUAWLVLn9crw+efQvHn4jEUkdSVozWyR+N17b/hClmg0bAhvvRUSorff3nrhW0mM888PzZG1NOmISEpTMiSR+OKLsGr3449HHUl2q1ULnn0WnnoqLOuxbBnk5m7/dbJ9P/8cnmvWjDYOEdk+JUOSdO7hL+b69bV2VirYYYeQFP36a5iHaODAMBxcym/DBujcGa6+OupIRKQslAxJ0j33HLz7bpgMsHHjqKORfNWrh+Hfb74JBx4IS5dGHVH6evrp8PlppJ5Ieoh7aH220ND6yrF2bfiLuWFDmDFDy26kopdfDh1+W7WCf/1LC+bGyz0svbF2bWgO1igykdRQ2tD6uGqGzKypmd1tZl+a2eJYWS8zu8fMWiYiWMls1avDn/8Mo0crEUpVAwfCO+/A8uVw5plRR5N+PvwwzCl03nlKhETSRZmH1ptZG2Aa0BRYB2yK7VoM7Ap8ZGb7ufu3CY9SMka1aiEZktS2775hpF/dulFHkn7+9jeoVw+GDIk6EhEpq3hqhkYCc4Cu7l6XkBDh7vPcvT9hmY7rEx+iZIqhQ0NfCkkPnTtD69aweTMMHgz/+EfUEaWHu+4KC7PWqxd1JCJSVvFMutgX6OnuP5Ww/3bgq4qHJJno9dfD+mMdOmz/WEkt69fDwoUhIfrpJ7j44qgjSm0NGkDv3lFHISLxiKdmqHopiRDungvUqXhIkmk2boQLLoCOHcNirJJe6taFiRPDkh2XXBKSIU3QuK3cXDjssPBZiUh6iScZWmFme5W008wOApZXPCTJNHfdBd98E2acrl496mikPGrUCBMzDh8Od94J554bdUSp56WXQiK0adP2jxWR1BJPM9kTwEtmdhXwFoCZ1QNaAccAFwO3JTxCSWtLloT5hI48MiwKKumrShW4557Qj2jffaOOJvXcdx+0bx9qh0QkvcSTDN0E7A08BORPTrQy9mzAK4R+QyIFmjQJS2507x51JJIIZnDppVu2//536N8/LEaazT77DKZMgdtv15QRIumozM1ksT5BhwMnA68TOkt/RUiCBrv7ke6+uVKilLSUlxe+PI8+WhP3ZaKFC0NfsH33ha+/jjqaaI0ZE9YgO+20qCMRkfKIa9JFD55w9wHu3iX2OMrdn6ysACU9bdoUviTHjo06EqksLVrApElhpuXf/x6mTYs6ougccABce22YWV1E0k9C1yYzs0cTeT5JXw88EL4cmzaNOhKpTHvtBf/5T5hT58ADwxQK2ejEE+GKK6KOQkTKK54+QwCYWTNCp+naxew+vMIRSdpbujSs1t23b2gik8y2667w3//C4YfDggVRR5NceXlhxukTTlCtkEg6i2c5jvbA48A+JR3Clo7VKcPMdgP+AnxMWDbkQ3f/e7RRZba//AXWrAlD6bU2U3Zo1gymTg3LrUDoQ7Trrpn/+584Ec45Bxo1guOPjzoaESmveGqGHgKaEZbl+BH4tch+A+5KUFyJ1AR40N0nmVk1YLGZPe/uK6IOLBP98AOMHx8mV+zcOepoJJnyE6Evv4Q994RTTw0JcSaPrrrvvjCSTjWgIuktnmSoJ9DR3ZeUdICZHRVvAGZWHbgBuATo4O7fF9k/ELgaWA9UBS5w9+llPb+7TypSlMuWRWYlwVq3Dn2Fdt016kgkKh07hkkZ//pXWLQIJkwII60yzZw5oWbo+uu3JIIikp7iSYbmAD9v55i4/j4ys3bAk8DXhESn6P6ewD+AXu4+28yOAP5lZl3dfVHsmOkU/z76u/vCImVnAze5+y/xxClls3Il7Lgj9OwZdSQSpSpVwnw7LVrARReFyTZfeincG5lkzJiQBA0dGnUkIlJR8Ywm+wtwpZmV9pp4u0/WBYYA40vYfyXwL3efDeDurwKLgYLFANw9x917FPPYKhEys2OBeu5+T5wxShmsXAm//W34EhQB+POf4ckn4YMP4P77o44m8X78MazX1qxZ1JGISEWVuWbI3d8ys+bAPDP7FFgGFF2usUE8F3f3WQBm1qqEQ/oCtxYp+wg4mNB0ViZmNgho7u6jzKw7sN7ds3yauMS67rqw9MZBB0UdiaSSE0+ETp1g993Ddl5eqDnKBM8+GxZnFZH0F89oslOAhwkdpVuWcFjCRpOZWUNCcvVTkV2LgDKv/mNmBwBjgU/MbADQGDiP0DQnCTBrFoweDcOGhY6zIoXtsUd4/vFHGDAgNC+l89pm7qEvVPPm6iskkini+RvtGuBeoD1Qw92rFH0AJXauLoc6seeNRco3UvwcR8Vy90nu3sDd+8Qeu7n75OKONbOhZjbdzKYvXbq0fFFnGXc4/3xo0CAsyCpSkk2bYN26UHv40ktRR1N+kyeHgQLvvht1JCKSKPEkQ/Xd/c/uPi+2TllxLkpEUDFrY881ipTXANYl8DoF3H1crA9STpMmTSrjEhlnzhz48MOQCDVqFHU0ksratYN//zs0mR1zDIwbF3VE5TN6dOgMvk9JM66JSNqJZzTZh2bWLH8UVwkS1kzm7svNbCVhbqPCmgHfJuo6UjEdO4YJ9rTshpRFkyahRuX440OzavXq8Kc/RR1V2c2fDy++CJdeCrVqRR2NiCRKPDVDw4DbzewwM2tcwjF3JiCmwt4GcoqU5cTKJWKzZoVmshYtMntiPUmsOnVCQnHddXDUUREHE6cHHgjPZ58dbRwikljxJEPzgEHAq4RZnDcXfQCJrh+4BTjEzDoDmFl/oDkwJsHXkTh9802YT+jmm6OORNJRtWphlfcdd4T168NadusqpfE7cTZvDrOrDxwIbdtGHY2IJFI8zWQbgX+Wst+A4+K5eGz26TeBHWNFT5nZQnc/BsDdZ5jZYOAxM8ufgfqQ7TTVSRL8+c+hiePUU6OORNLd5MkwahS88w688krq9j2rWjX0j1u/PupIRCTRzL1s3XzM7Cd3b17RY1JdbPj9gA4dOpw5Z86cqMNJSa+9BkccAbfdFvpOiFTUc8/B4MHQvn1Y4kI1LyKSaGY2w92Ldr0B4msm612GY7rEcb6U5O6vuPvQBg3imj8ya2zcGBZh7dQJLrgg6mgkUxx7LLz5Jvz0U5iD6PPPo45oa9OmhWVF5s6NOhIRqQxlTobKOGPzqArEImng229hw4awGnn16lFHI5lk//3h/ffDiLNUu7dGjw7LijQuaeiIiKS1eGagblOGw44Fzil/OJLqunQJcwtl4irkEr1u3eDjj8OSHe4wYwbkFFupnTyLF8PTT8NZZ0G9etHGIiKVo8SaITObZGafmln+oOnvgbnbeejvpgz2wguhmUyJkFSm/LXLHn4YevWC++6LNp4HH4Rff4Vz9GeeSMYqrWbI2DpZWgv8dTvHJ3IGakkhkyeHWYPvuAMu0m9ZkmDQIHj11bDcy8KFcNNNYJbcGHJzw9xC/fqFfnIikplKTIbcvU+RojXufn1pJzOzYYkISlLLpk3hC6ldO002J8lTq1ZYGf7cc+GWW0JC9NBDyV0cddOmMFBACxCLZLZ45hnqVtIOM6vt7utKOyZdFBpaH3UoKWPs2DC65/nntQSBJFfVquH+a9kSbrgh9NtJ5ppgtWpp+giRbBDP0PobStl3q5ktBXpVMJ7IaWj91pYsCbMDH3xw+i2dIJnBLNyDs2dvSYR+/bXyr/vll/DEE8m5lohEK55k6KhS9t1A6C90R4WikZSzYgV07RqG0ie7v4ZIYbvuGp5feCGsfP/dd5V7vbvugqFDYc2ayr2OiEQvnmSoRO6+FHgCaJiI80nq6NQJ/vMf+O1vo45EJGjWDJYuDbVEH39cOddYsSLUCg0eDA31v5pIxiu1z5CZPVxos4GZ/Z0waqy483QCFiYwNolQXl7otDp0qCaak9Syzz4hQT/kEOjdO/RlO/jgxF5j/PiwcOy55yb2vCKSmkpdm8zM8gptOsUnQgDrgP8B57v7B4kLLzo5OTk+ffr0qMOIzMMPw+mnb/nrWCTVLFwIhx0W+hLNnBmacxNh82bo2BFatAgzYotIZihtbbJSa4bcvaAZLRMWYZWyWbkSrrgCfv/7MNeLSCpq0QKmTIEJE8LM6Iny009Qpw6cd17izikiqS2ePkNXVloUklKuvRZ+/jmsx6RO05LKGjQIM0ObwaxZMGJEaOKtiFat4NNP4bjjEhOjiKS+eBZqfaQS40gZZjbAzMatWrUq6lAi8fnnMGYMDBsGPXpEHY1I2T3/fJil+qSTwrIx5bFsGaxeHZKrKgkZXiIi6UD/3IvI9nmGdtwRhgyBG2+MOhKR+Fx9dej0/+STcPjhIamJ16hRsPPOofO0iGQPJUOyldatw0iaRo2ijkQkPmZw+eXwyCNhLb3evcMQ/LJasyYMHDjkEKhdu7KiFJFUpGRIgPBFcNJJ8NVXUUciUjGnnBIWeG3dGurXL/vrnngi1Cap47RI9lEyJEDoazFhQug4LZLuDj0UXn4ZatQI9/T2ZslwDwMG9twTfve75MQoIqlDyZAwZw7ccUfoK7TvvlFHI5JYw4eHJrPXXiv5mBkz4IsvwrEaQSmSfZQMCRdeGP6CvvXWqCMRSby77grLyRx5ZOgPV5ycHPjkEzjhhOTGJiKpQclQlnv99fC49lporik1JQM1bRo6VB94IJx2WmgSLm7i/R49oFatZEcnIqlAyVCW22+/MBx5+PCoIxGpPPXqhU7VgwfDQw9tPez+hhvg1FMrPlmjiKSvUtcmy0ZmNgAY0KFDhzPnzJkTdTiVyl39IyS75OWF4fZNm8Jjj8FVV8EPP4QaoQcf1Dp8IpmstLXJVDNURLZMuvjDD7DXXqHjqEi2qFIlJEITJoQmsx9+COXr18PQoaFcRLKPkqEsdcklYfRM48ZRRyKSfCNGhNXpC1u3LpSLSPZRMpSFJk2Cp5+GK6+Etm2jjkYk+ebPj69cRDKbkqEsk5sL558P7drBpZdGHY1INNq0ia9cRDKbkqEs849/wKxZcPfdGkYs2WvUqG3XH6tdO5SLSPbZIeoAJLlOOimsTD9wYNSRiEQnf9TYiBGhaaxNm5AIaTSZSHbS0PoS5OTk+PTtLWiUZtavV22QiIhkJw2tF6ZNC52lp02LOhIREZHUomQoC+TlwXnnwQ47QJcuUUcjIiKSWtRnKAs8/DBMnx4mlKtXL+poREREUotqhoowswFmNm7VqlVRh5IQK1aE+YT22w/+7/+ijkZERCT1KBkqItOW43j2WVi+HO67T+uQiYiIFEfJUIY788wwr1D37lFHIiIikpqUDGUod5g3L/zcuXO0sYiIiKQyJUMZ6qmnoGNH+PDDqCMRERFJbUqGMtCaNWFV+t12g549o45GREQktWlofQYaNQoWLgydp6tWjToaERGR1KaaoQzz9ddwxx1wyimwzz5RRyMiIpL6lAxlmPffDxMr3nJL1JGIiIikByVDGeb002HuXGjWLOpIRERE0oOSoQyxYUOoFQKoXz/aWERERNKJkqEMcccdsP/+8PnnUUciIiKSXpQMZYD588MIsj/+Ebp1izoaERGR9KJkqIh0XKj1kkvC8x13RBuHiIhIOlIyVES6LdT67rvwzDNhZfo2baKORkREJP0oGUpzS5bAnnvCpZdGHYmIiEh6UjKU5k48EaZPh5o1o45EREQkPSkZSlOLF8Ojj0JeHphFHY2IiEj6UjKUpq64As48E+bNizoSERGR9KZkKA1NnQqPPAIXXQTt20cdjYiISHpTMpRmNm+G886DFi3gqquijkZERCT97RB1ABKfhx+GGTPgH/+AunWjjkZERCT9qWYozbRpAyefHEaRiYiISMWpZijNHHJIeIiIiEhiqGYoTXz6KVxzDaxbF3UkIiIimUXJUBpwh+HDYexY2LAh6mhEREQyi5rJ0sCTT8L778O4cdCwYdTRiIiIZBbVDKW4X34J64717AmnnRZ1NCIiIplHNUMpbuRIWLgQnnsOqlaNOhoREZHMo2Qoxf3f/0GTJvC730UdiYiISGZSMlSEmQ0ABnTo0CHqUADo0SM8REREpHKoz1AR7v6Kuw9t0KBBpHG88goMGQKrVkUahoiISMZTzVAK2rABLrgAatWC2rWjjkZERCSzKRlKQbffDnPnwjvvQLVqUUcjIiKS2dRMlmLmzYObb4bjjoMDD4w6GhERkcynZCjFXHVVeP7rX6ONQ0REJFuomSzF3HZbqBVq0ybqSERERLKDaoZSxObNYQ2y5s1h4MCooxEREckeSoZSxL33wv77h+U3REREJHmUDKWARYvg2muhfn2oVy/qaERERLKLkqEUcMUVsHEj3H131JGIiIhkHyVDEfvgA3j0Ubj4Yth116ijERERyT5KhiJ2223QsiX85S9RRyIiIpKdNLQ+YhMmwDffQN26UUciIiKSnZQMReSXX8JSG7Vrw+67Rx2NiIhI9lIzWUSuuAK6dw+LsoqIiEh0VDOUZH36wJo18MkncO65ULNm1BGJiIhkNyVDSTRhAkydGobRV6kCu+0WdUQiIiKiZrIkmTABhg4NiRBAXh78+c+hXERERKKjZChJRoyAdeu2Llu3LpSLiIhIdJQMJcn8+fGVi4iISHIoGUqSNm3iKxcREZHkUDKUJKNGhTmFCqtdO5SLiIhIdJQMJcngwTBuHNSoEbbbtg3bgwdHG5eIiEi2M3ePOoaUYmYDgAEdOnQ4c86cOVGHIyIiIglgZjPcPae4faoZKsLdX3H3oQ0aNIg6FBEREUkCJUMiIiKS1ZQMiYiISFZTMiQiIiJZTcmQiIiIZDUlQyIiIpLVlAyJiIhIVlMyJCIiIllNyZCIiIhkNSVDIiIiktWUDImIiEhWUzIkIiIiWU3JkIiIiGQ1JUMiIiKS1czdo44hJZnZUmBeKYc0AFaV8/SNgWVxvqYi15PSpeNnW1rM5bm/EnXtVL1GIs5X3nOU53WpfE9W9v2VDKn8+ZYkypiT+W8+UfdXcTG3dfcmxR7t7nqU4wGMq8Brpyfzenpk3mdbWszlub9S7fNK9DUScb7ynqM8r0vle7Ky768kvYeU/XxTMeZk/ptP1P0Vb8xqJiu/VzL8etkkHT/bKGNOxrUTfY1EnK+85yjP69Lxnkwn6fj56t98JZ5PzWQRMLPp7p4TdRySmXR/SWXS/SWVKar7SzVD0RgXdQCS0XR/SWXS/SWVKZL7SzVDIiIiktVUMyQiIiJZTclQkphZDTO7y8xmmtl7ZjbNzI6OOi5JP2a2v5k9b2aTzex9M/vczC4otP8IM3vdzN4xs6lm9oaZ7R5lzJLazKy6md1iZpvMrF2hcjOzk2L30jtm9qGZPV34mNhxbc3sOTP7yMymmNlbZtYt2e9DUo+Z/cnMvoz9f1X40aDQMcXef7F9ZboHK2qHRJ5MSnUVcCTQ3d1/MbM9gKlm1svdP404Nkkvg4CZ7n4DgJl1Bz42s2/d/VXgEeB8d/9HbP8twDtmtpu7L44qaElNsS+VJ4GvgapFdtch3E8HuPv7ZlYVeAz4l5n1cPf1seOeAH4C9nb3vFhy/rqZdXD3jcl4H5LSbnH3R4rbsZ37D8p+D1aIaoaSpwfwkbv/AuDunxAmhDowyqAkLd0L3JW/EUumVwIdYkVT8hOhmDsIE5n1S1aAklbqAkOA8cXs2ww85e7vA7j7ZuBuoCNQeMRPD2CSu+fFtt8CWgFdKidkySCl3X9Q9nuwQpQMJc9zwH5m1grAzA4BmgD6S13i4u6z85NqM6tiZmcCG4FnYvuPKfKS/L+caiQvSkkX7j7L3b8pYd96dz+pSHFx99NzwNFmViu2PRjII/1nqpZKVtr9F9tf1nuwQtRMliTu/oiZ1QZmmdlPhKz2GeDpaCOTdGVmVwHnAcuB/u6+oIRD9wE2AC8nKzbJaPsAi4AphcpOJwyJXmhmq4GmwJXu/kME8UnqOcLMTgaqE+6dW939owqcr7h7sEJUM5QkZnYG8Bcgx907A3sAUwl/PYnEzd1HAs2BUcB7ZrZv0WPMzICrgavcfUmSQ5QMY2Y1gcuA89z910K7HgHaAa1jzycB/0tyeJKaFgPfEP5g+wPwAvCBmf2uPCcr5R6sEM0zlASxL6SfgXvc/fpC5ZOBt2NfaiLlZmavAA3cff8i5TcDrdx9SDSRSbowsz7AJKC9u39fzH4DJgDfuftVhcp7AJ8Afdz9vVhZTUKNZT93/3dlxy7pxczeB35296MKlfWhlPsvdkyx92AiqGYoOZoAvwG+L1I+Fzg26dFIWjOz6sUUzwa6FjnuQkIH1lOTEJZkvruAFcV8CXWMPX+fX+DuGwg1Avr/TYrzLVsGfMSjpHuwwpQMJccyQgfX5kXKmwPrkh+OpLkZsb+QCmsBFPQZijXL9geOd/dNZrazmfVNZpCSOczsRsIfdOfFtnuaWc/Y7vz7rnmh46sS/gjU/29ZzsxujvWXLawlMD/O85R2D1aYkqEkiA03fRQ4zcwaApjZnsBBqAO1xK8esf8QIPynAPwR+Hts+0RgBKEvUTczywEOBv6Q/FAl3ZnZpcAA4H6gZ+x+GgDkT6r4IfAFcFEsCQI4lzDS5/kkhyupZx9CB3sAzKw3YUqZ+8t6gjLcgxWmPkNJEsuMrwP6Ev5aqkdIkO5y/RIkDmY2CDgDqEWYg6MWIREa6+5uZrkUP1L0ene/LmmBSlqINbu+CewIdAemAQvd/Rgzawn8WMJLT82fSM/M2gK3A+2BTbH9N7r765UYuqQBMzsUGE74zqtC+L/pr+7+bGx/ifdfbH+Z7sEKx6nvYREREclmaiYTERGRrKZkSERERLKakiERERHJakqGREREJKspGRIREZGspmRIREREspqSIRFJGjObaGbLzczN7E9Rx1MeZlbHzMaY2XwzW2xmc8xsaARx3GBmC82sWbKvLZJplAyJSNK4+6HAMVHHUUHXE9Z7O8DdmwLjgX0jiKMh0IAw03O5mdn3sUWjRbJWcbPUiohIyQ4CPnf3b2PbdwBF115KhuHAZe6u9b9EKkg1QyIi8fkNsD5/w903uvuKZAfhgRIhkQRQMiSSIcysr5ktMrNfY00f3cxsspktM7NvzezPRY5/vLj+O7F+Pati5X2KlOcff5qZXWdmc81spZk9FutL09jMnjezJbG+NCeXEnJVM7s2dtxqM5thZv2KeV8NzOxOM/shdv0fYn12GhU65prYe3cze8TMjo6db3VZ+ieZWX0z+2vsc1tsZvPM7C4za1DomFPNbBHQGtg3dr1FZta3hHN2iu1fH4uhm5m9Fevns8rMJpjZTsW8rqeZvRqLY4mZTY+tR1f4mPdL+N19Ueh3d6iZ3Rt7Lz+b2TNm1rjQsX1LeD+vFTrmaDP7wMwWxOL+wMyuMLMdS/s8RdKOu+uhhx4Z9AAmA8uAZ4AmhD96rgQcOLrIsX1i5X8qUv6nWHmfEo7/BDgSMGBvIBcYA9wLtACqxn7OAzqVcI7/AecB1YA6wITYeXoXOrYW8DHwJfDbWNlvY6+dDdQpcm6Plf8tds66wFdF31+R1+Rf4ytg11jZrrHtmUDtIsd/D0yO4/fxSCyuD4HusbKuhMUnPwNqFPlsNsQ+y1qxz/dkwoK8I+L83U0F+sbKdgNWAC8VE1+x7wc4gLDo6hGx7SrA0OLuCz30SPeHaoZEMlMjYJS7L3X3POBOwhfbUQk6/7fu/pIH04D/AMOAN919obtvBu4mfJkfW8I5Frj7aHfPdfe1wAWEL9o7Cx1zMbAHMNzdvwSIPV8MdI5ds6imwIXuvtbd1wAXEhKRkuRfY4S7z4ldYw4wgrCK9mWlfxRldq+7fxo7/xfAbUA3QoKBmVUBxgFrgYvcfX3s830MeAO43sx2ieN6U9397dj1ZgGvAodZWCW8LAYCq9z91dg58tx9XOw8ap6TjKJkSCQzrXf3mfkb7r4RWEqotUmEosnFQkJtUOHyH2PPrUo4x7uFN9x9Wez1e5pZfpwnAL8SarsK+yj2fGgx5/3I3Qv36XnD3WeXEEP+NQBeK1L+CiE5O76U18bj3SLbr8eeB8Se9yDUSL0V+30V9jLh8y0psSzOf4ts/0CohWtSxtcvARrGmgub5he6+wB3Ly25FEk7SoZEMtOyYsp+JXwZVsb5fy1a7u75ZSWNtFpcTNmi2HP72HMHwqjXHwr1aVkEfE6oQWlUzDmWbCf2ojoAKwsnUFCQQC6P7U+Eou+3uPcK8FMxr11Y5JiyKOl3VNZ74D5CQnghsMDM3jWz082sVhwxiKQFDa0XyUx5FXz99v5QKvb8sSa5RPvF3eOZWLAyYkhHFfocYk2MA81sd2AQcBLwEHCpmfV29+KSWZG0pJohkey2KfZc9A+jbUY5VYKmxZTlJz1zY89fAw3MrE7RA83st2bWPQFxzAF2LFrjYWY1CBMbzknANWDb91v0veZfp7imzBZFjql0ZlbFzMzdP3P3K4C2wDVAJ0LHd5GMoWRIJLstiD23LlK+TxKu3afwRmzYdy9ghrvnNwv9M/Z8TJFjqwDPAsUOa4/T07Hnw4uUDyB0AH+axOhTZDv/ei/HnmcSkp2+sUSssIGEEWXPJyiWwtYSS4bNbAczG21mbYCHgUvyD4p1ir8rtrljJcQhEhklQyJZzN3nAl8Ag82sVaw24P+ANkm4/O5mdo6ZVTOz2sA9hOTj4kLH3AlMB0aZWU8AM6tL6M9SldBsU1F3EKYKGGlmu8au0QEYCXwK3J6AawAMiTU5YWZdgUsJfZ8ehIImxjMJUwLcaWY1LTgZOAy41rfMep1I/wN2NrOahCR4KFsmlTw3FitmtgOhRiiPxCWIIqkh6rH9euihR2IehGHgiwgdZTfHfu5LmC9mUazs19jPXQu9bhfgTWAl8A1wFXAaYSTVcuC12HGPx7YdWAW8REhIFhG+PD3286DYY1GsbH3s54bAxELnOIuQ7HwHrAZmAP2KeV/1CMPQvyd0Qv4OeABoXuiYs4u53tlxfHb1CUlR/jXmE2pBGhQ65tRiPsfXynDuR2JxtSXUAv0Y+6z/AexUzPE5hJFtSwgjAKcDg4sc836R38VnsfLJse38392DsfLPgDWx8qXADYXO9VvCyLOlhCa7M2LlPYCxhDmefiJ04n4HOCjqe10PPRL9MHdHREQqh5k9Apzi7hZ1LCJSPDWTiYiISFZTMiQiIiJZTcmQiEglyF+oldgM17EJIy+POCwRKYb6DImIiEhWU82QiIiIZDUlQyIiIpLVlAyJiIhIVlMyJCIiIllNyZCIiIhkNSVDIiIiktX+H+MaIFX8xAKmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (9, 6))\n",
    "\n",
    "plt.errorbar(ns, tm1, yerr = tm1_err, fmt = 'bo--', label = 'stabilized sinkhorn')\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.yticks(fontsize = 14)\n",
    "plt.xticks([8, 32, 128, 512], fontsize = 14)\n",
    "plt.gca().get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())\n",
    "\n",
    "plt.xlabel(r\"number of points\", fontsize = 18)\n",
    "plt.ylabel(\"time elapsed (s)\", fontsize = 18)\n",
    "\n",
    "plt.legend(fontsize = 16)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}