Figure_1_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "417f4911",
   "metadata": {},
   "source": [
    "# Figures 1 and 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea90ba86",
   "metadata": {},
   "source": [
    "Here, using XGI, we reproduce figures 1 and 2 of the paper  \n",
    "\"Do higher-order interactions promote synchronization?\"  \n",
    "Yuanzhao Zhang *, Maxime Lucas *, Federico Battiston  \n",
    "[arXiv:2203.03060](https://doi.org/10.48550/arXiv.2203.03060)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8da1197d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sb\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "\n",
    "import xgi\n",
    "from utils import (\n",
    "    compute_eigenvalues_multi,\n",
    "    random_flag_complex_d2,\n",
    ")\n",
    "\n",
    "sb.set_theme(style=\"ticks\", context=\"paper\")\n",
    "\n",
    "fig_dir = \"figures/\"\n",
    "results_dir = \"results/\"\n",
    "\n",
    "Path(fig_dir).mkdir(parents=True, exist_ok=True)\n",
    "Path(results_dir).mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b70265e6",
   "metadata": {},
   "source": [
    "## Fig. 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dee920f0",
   "metadata": {},
   "source": [
    "### Generate random structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f568cc20",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100  # number of nodes\n",
    "ps = [\n",
    "    0.1,\n",
    "    0.1,\n",
    "]  # ps[i] is the wiring probability of any i+2 nodes (ps[0] is for edges, e.g.)\n",
    "\n",
    "n_repetitions = 10  # number of realisations of random structures\n",
    "\n",
    "# generate random hypergraphs\n",
    "HGs = [xgi.random_hypergraph(N, ps) for i in range(n_repetitions)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8a6dcbc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate random simplicial complex\n",
    "MSCs = [random_flag_complex_d2(N, p=0.5) for i in range(n_repetitions)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "326410db",
   "metadata": {},
   "source": [
    "### Compute Lyapunov exponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "84b71678",
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = np.arange(0, 1.01, 0.1)\n",
    "n_alpha = len(alphas)\n",
    "\n",
    "lyaps_HG = np.zeros((n_alpha, N, n_repetitions))\n",
    "\n",
    "# compute Lyapunov exponents for all alpha values\n",
    "for j, HG in enumerate(HGs):  # for all hypergraphs\n",
    "\n",
    "    for i, alpha in enumerate(alphas):\n",
    "        lyap_multi = compute_eigenvalues_multi(\n",
    "            HG, orders=[1, 2], weights=[1 - alpha, alpha]\n",
    "        )\n",
    "\n",
    "        lyap_multi = np.sort(lyap_multi)[::-1]\n",
    "        lyaps_HG[i, :, j] = lyap_multi\n",
    "\n",
    "lyaps_MSC = np.zeros((n_alpha, N, n_repetitions))\n",
    "\n",
    "# compute Lyapunov exponents for all alpha values\n",
    "for j, MSC in enumerate(MSCs):  # for all simplicial complexes\n",
    "\n",
    "    for i, alpha in enumerate(alphas):\n",
    "        lyap_multi = compute_eigenvalues_multi(\n",
    "            MSC, orders=[1, 2], weights=[1 - alpha, alpha]\n",
    "        )\n",
    "\n",
    "        lyap_multi = np.sort(lyap_multi)[::-1]\n",
    "        lyaps_MSC[i, :, j] = lyap_multi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4d4c0e17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# average and std over the random realisations\n",
    "# consider the second largest exponents only\n",
    "means_HG = np.mean(lyaps_HG[:, 1, :], axis=1)\n",
    "std_HG = np.std(lyaps_HG[:, 1, :], axis=1)\n",
    "\n",
    "means_MSC = np.mean(lyaps_MSC[:, 1, :], axis=1)\n",
    "std_MSC = np.std(lyaps_MSC[:, 1, :], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd87fef7",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "34c89003",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAADwCAYAAACwqvfHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA330lEQVR4nO3dd1yV9fv48dcZ7KVscOIAxQJXloqZWJFaao4caZ92hvTTUrNhpabliNIc1SdN/fbJkalpmpYDVDLNreBK3CLbxTjAOdy/P5ATxBDk4Dng9Xw8fCA39zn3dTycy/e8bpWiKApCCCFMQm3uAIQQojaRpCqEECYkSVUIIUxIkqoQQpiQJFUhhDAhSapCCGFCklSFEMKEJKkKIYQJSVItIjs7m7i4OLKzs80dihCihpKkWsSZM2fo168fZ86cMXcoQogaSpKqEEKYkCRVIYQwIUmqQghhQpJUhRDChCSpCiGECUlSFUIIE5KkKoQQJqQ1dwDliY2NZdKkSWRmZuLg4MD06dNp0qRJifOuXLnCxIkTuXLlCnq9nq5duzJ27Fg0Gk21xJWZmclnn31W7NjYsWNxcHColusJIWoOi02qubm5REREMGXKFEJCQti6dSsRERFs2LABlUpV7Nz333+fhx56iFdffRWdTseLL77IihUrGDp0qJmiF+Yk/+kJc7LY7v/Ro0fRaDSEhIQA0L17d7Kzszly5EiJc3v37s0zzzwDgK2tLf7+/ly+fPmuxiuEEGABLdXt27fz+uuvlzgeHh6Oj49PsWNeXl4kJCQQHBxc7Hjfvn2Nf4+Li+OXX35h0aJFZV4zOTmZlJSUEsfj4+MrGb2wVAZFzdW8OuTkW2OjzkWXo0caquJuMHtS7dq1K8eOHStxfP369ezZs6fEcbW67Mb1tm3bePfdd/nwww8JCgoq87wVK1Ywd+7cOwtYWLyT569x5GYgBjSoUFBQER4Zw8SXO9LSz9Xc4ZVr/vz5XLhwgWnTppn8uVevXs26detYvHixyZ/bEgQEBLB9+3a8vb3NGofZk2pZfH19SUpKKnYsOTm5ROu10Ny5c1m6dClz5syhQ4cO5T73oEGDCA0NLXE8Pj6ecePG3XnQwuyyc/RM/+EQBjSACgWV8fikBX+y6MMw7Gws9tde1AIW+9sVFBRETk4Ou3btolOnTkRHR6NSqWjVqlWJcyMjI9m6dSurVq0qM+kW5enpiaenZ3WELarRst9PsvNQ+WPlGdl5ZOn0QPHJTEWBTJ2e16ZtxdHOqszHd2ldjyGPB9w2lj179jB58mRcXFy4evUqy5cvJyoqiqVLl6LX60lLS2PAgAG88cYb7NmzhxkzZuDv78+JEyfIzs5mwoQJhISEkJmZyYcffkhsbCyurq7UqVMHFxcXAJKSkpg8eTLnz59HpVLRtWtXRo0ahZWVFaGhoTz11FNERUVx9epVXnjhBVJSUti7dy83btzg008/pV27diXiTk9P5/XXX+fy5cvk5uYyefJkOnToQK9evXjzzTd59NFHAfjiiy/Q6XQEBASwdu1atFotSUlJODo68sknn9CkSRNyc3OZNWsWf/75J/n5+TRq1IgPPvgADw8Phg8fjouLC2fPnqVnz54MHTqUCRMmcPr0aZycnGjevDmKojBt2rQS5w4ZMoSPPvqI5ORk0tLScHFxITIyksaNGzN8+HD8/Pw4fvw4aWlphIaG8s4776DVFqSyb775hsOHD5Oens5TTz3FmDFjbvtemprFTlRptVq++uorZs+eTa9evZg7dy5z5swxLpN65ZVX2Lp1K2lpaSxYsICsrCxGjBhBnz596NOnD7NnzzbzKxDmYDDkV+nnlXH69Gk+/fRTNm7ciLW1NUuXLmXevHmsXr2aZcuWMX/+fNLT04GCsf4hQ4awZs0ahg4dyqxZs4CCHpbBYGDjxo18++23XLp0yfj8Y8eOJTAwkPXr17Ny5Uri4uL49ttvjT+/du0a69atIzIykunTp9OuXTt++uknnn76aRYsWFBqzBcvXmTMmDGsW7eOZ555xhjH0KFDWbFiBQB5eXmsWrWKIUOGAHDo0CHGjx/P+vXreeyxx4y9uQULFpCbm8uqVatYu3Ytbdu25f333zdey87Ojg0bNjBy5EimTp2Kl5cXmzZt4rvvviM2NrZYXEXP3bhxI/7+/qxYsYItW7YQEBDA999/X+w1/PDDD2zYsIHDhw8b4wZwd3dn9erV/PjjjyxatIgLFy5U6j01BYttqQIEBgYW+wcrqugv1/Hjx+9WSMKMhjwecNtW5OY95/nyx0Nl/vyFJwN5tEMjk8Tj4eFBo0YFz2VnZ8e3335LdHQ058+fJz4+HkVRjAXPPTw8jOP8gYGBxnHNmJgYxowZg1qtxtHRkaeeeoozZ86QlZXFvn37+Oqrr4CCVS1Dhw7lv//9L+Hh4QA88cQTADRs2BAomJ8o/H7Hjh2lxtyqVSuaNWtmjGP58uVAwWTvrFmzuHz5MocPHyYgIIDGjRtz4MABOnbsiL+/PwCDBw9mxowZXL16lW3btnH9+nX27t0LQH5+frEC7+3btzf+PSoqiuXLl6NSqXB2dqZfv36cPHmy1HOfffZZDhw4wOLFi7l48SL79++nbdu2xp8PGTIEa2trAPr06cPOnTt59tlngYKVQFDQG/Xw8CAlJcX473O3WHRSFaKyQlrXY8G6WLJ0eRQdAlCpwN5GS+fgeia7lp2dnfHviYmJDBo0iP79+xMcHMzAgQP57bffUBQFKEiKZSk8B8DKysp4rOhxKEhaBoOhxLllfV+awm4ygEqlMl7DwcGBPn36sGrVKg4ePMiwYcNKfUzh+RqNhvz8fN5++20ee+wxAHQ6HTdv3jSeW/Q1a7Va8vP/6SX8e6150XMjIyPZt28fzzzzDG3btsXBwYHk5GTjz4tu6lEUpdjkddF/g6Kv726y2O6/EHfCzkbL+GGt0WAAFFTkAwp2Nlo+erljtU1SxcbG4uDgQEREBF27dmXHjh0oilIskZTmkUceYdWqVej1erKzs9m4cSNQkOTatWvHkiVLgIKEtXz5cuO67erw7LPPsmbNGi5evEi3bt2Mx3fv3m1c971s2TLatWuHs7MzXbp0YenSpeh0OgA+++wz3nnnnVKfu1u3bqxcuRKArKwsfvnllxKJtdCOHTsYNmwYTz/9NA0aNCAqKqrYfyZr167FYDCQlZXF2rVri8VqCaSlKmqdgIZ1CHI6Vmyd6sQxz+Hm6lJt1wwJCWHNmjWEhYXh4OBAs2bN8PPz49y5c9jY2JT5uJEjRzJ58mR69OiBq6trsa7qZ599xscff8xTTz1Fbm4uXbp0ISIiotpeg5+fH76+vjz88MPFWn9eXl588MEHJCUl4e3tbdytFh4ezsyZM+nfvz+KotCwYcMyl4K9++67fPTRRzz55JM4Ojri4eFRZuv9jTfeYMaMGSxYsACVSkXr1q2LDRXo9XoGDhxIRkYGPXv2ZMCAASb8V6g6lWKO9rGFiouLo1+/fqxevbrUVQaFZBukZZP3584kJiYyYMAA1q1bh6trwXpeU61t/eGHH/D39+eBBx4gPz+f8PBwunXrxqBBgyr1PMOHD2fAgAH06dOnSvFUJ2mpilrHwcGBjz76yNxh1CiF67zfeustY0I1JX9/fz799FPy8/PJy8ujY8eO9O/f3+TXsQTSUi2iMi3V6TM/v6vdSyFEzSATVXegcBvkOV0DEnM9OadrQHhkDMfPpps7NCGEmUlSraSS2yDVgMq4DTI7R2/uEIUQZiRJtZJiDl0udxvkH4el5KAQ9zKZqKqkxPQsNGoVhvySQ9EatYrEtCwzRCWKktl/YU6SVCvJ29W+1IQKYMhXqOtc9prEipCEIETNJt3/SgppXQ97Wy1QemJdFXWaE+dlwsrcDKp8rrpkkeR+k6suWej0OeYOqULmz59f5q4kU5g4cSKhoaFMnTqV0NBQ9u3bV23XuhsuXbpEYGCgucMoRpJqJZW1DdLeVkvfrk25fjOH8XNjWLH5ZJktWlG9/k4/y6lmKSR43yDVLZME7xuM2fIxJ1Plzg7Lly9n0aJFxapJCdOS7v8dKG8b5KMPNGTm//bxv00nOHgqhbeGtsWzrr25Q75n6PJ0zPprIflqpdhcYrY+h093zOPrpz7B1qrs4iYVVRPrqQ4ePBhFUQgPDy/RGv75559LjR3gu+++Y/ny5Tg4OPDQQw+xZMmSUu/WcfToUaZMmUJGRgZqtZqRI0fyxBNPEB8fz5QpU0hNTQUKKmK99NJLQEGlrBEjRrB582Zu3LjBqFGjOHDgALGxsRgMBmbPnk2TJk0IDQ0lNDSUgwcPcvXqVfr378/IkSNLxLBmzRq+//578vPzcXBw4O233yY4OJj33nuP9PR0vv76axISEhgwYABffPEFDz74YJV/F/5Nkuod0qjycbf+p5tve6tQRyMfZyJHd2Xx+jjWx5zl/0VG88bA1nQO9jVXqLXGT3Eb2HVhf7nnZOZmkaXX/XtxBgoKWXnZjPp1Ig7WZf8n16lhOwa06lWheE6fPs3vv/9Oo0aNyM7ONtZT9fDwICEhge7duxtL0sXFxfHRRx8RFBTE//3f/zFr1ixCQkKK1VPNyspiyJAhxqQ6duxYHnroIebNm4dOp+P111/n22+/NZb+K6yn+tdffzF8+HDmzZvH+PHj+eqrr1iwYEGJpLp8+XICAgJYuHBhsVuOlBf7iRMn+N///sdPP/2Eq6srs2bNKlbcpJBeryc8PJwJEyYQFhbGuXPnePHFF+ncuTMjRoxg5MiR9O3bl+vXrzNs2DDc3d3p06cPBoMBW1tbfvnlF9asWcN7773HTz/9RKtWrZgwYQJLly5lwoQJAGRkZLBy5UquXbvGgAEDipUxBNi3bx/Lly/n+++/x8HBgbi4OF599VW2bNnCBx98wMCBA1m2bBlr167lpZdeqpaECtL9rxY2VhpeezqID196EI1axbT/28uXKw6ikzWs1U6fX/IDX5mfV0Zp9VR37drFl19+yfTp08utp1pYvDomJoa+ffsWq6cKGOup/uc//wH+qacaFRVlvH559VSvXbtW4ddRXuzbtm0jLCzMuHX1+eefL/U5Tp48SW5uLmFhYQA0btyYbdu2kZyczI0bN4w353RxcaF///6lvo4GDRrg7u5u3M3479cxfPhw1Go1rq6uhIWFsXPnzmIxbNu2jUuXLjF06FD69OnDe++9h1qt5vz589jZ2TFr1iw+/fRT6tSpY2wpVwdpqVajBwK9mTO2G7OWHWDzXxc4djaNsc+2p1mDOuYOrUYa0KrXbVuR287s4uu935f582HBT9OtSSeTxFMT66mWprzYtVptsTjKuvGmVqstUcovPj6+1NKH/34dhQWnb/cayqujWvi8PXr0MLZsARISEvDy8gLg3LlzODs7c+bMGTIyMnB0dCzzWlUhLdVq5upsy8RXOvJS71YkpWcxbs4OVkedJl8msapFpwZtsdfallicoUKFvZUdHRu0Lf2BVVST66mWF3u3bt34/fffuXr1KoDxTgH/5ufnh729Pdu2bQMKbnkyePBg3NzccHZ25ueffwbg+vXrrFmz5o5ex+rVq4GC+2z99ttvJW7e2aVLFzZt2sSVK1cA+P3333n66afJyckhISGBDz/8kPnz59O5c+diidfUpKV6F6jVKvp2bcb9Td2Z+b/9LFofx8GTyYwe0gY3F7vbP4GoMFsrW0Z3eJlpO+cVTFbdYmdlw7sPjzTJJFVpanI91fJif/jhh3nxxRcZNmwYNjY2BAYGFmuhF7K2tmbevHlMnTqVWbNmoVarmTlzJq6ursyfP5+pU6eycOFC8vLyePLJJ3nmmWcqHWdycjL9+vUjMzOTl156iY4dOxa7p1fnzp2JiIjg1VdfNcb09ddfY21tzVtvvcXQoUMJCgqiefPmPP300yxbtsx4Hy5TkipVRdyNeqq6HD0L1sXy2+7zONlbM2pQax687587wMri/6rLzMxkeuQMbjjryLUyYJ2n4YMXxuPmYvqSdrXd0aNH+euvv4xjkBs3bmTJkiVltlirS2hoKDNmzCh2LytLJS3Vu8zWRkvEwNa0DfBkzo+HmLLoL3p2asyLve/Dxkpz+ycQFaJR1NS9/s8sv622ajvd7lV+fn4sXLiQXr16odFocHFxYerUqeYOy6JJUjWTTkG++Desy+dLD/DrrnMcjU9j3LB2eLrIW1JVUqTadBwdHY23sTanwrHamkAmqszIvY4dH4/oxHM9W5KQksGY2TvYuPsC+nw1qbmuXNZ5k5rrKkuxhKhBJKmamUatYmB3f2a80QV3FzuW/HqKQxn3SQFsIWooSaoWwr9hXaaN7IxGo4Jbha+lALYQNY8kVQuy/0QyBkPJxRhSAFuImkNmRe5AdU2ESAFsIWo+aalakNsVwE6+mlVi66IQwrJIUrUgZRXAVgEqFUTtv8TURX9xPaNmFFwW4l4k3X8LUlgAe/KC3RjQoEJBQYWdrRXjn3uArXsvsOPgZf6OjOLNIW1p7e9Z6WvIji0hqpckVQtTXgHsNv4etGvhxderj/DBN3/y9CPNGN6jJVZa6XAIYSkkqVqgsgpgq1QqQts3INDPlc/+t5810ac5/HcK44a1o76nk7nCFUIUIU2cGsjbzYFpESEMetSfcwnXGf3Fdn7bfU4msYSwAJJUayitRs2wHi35JDwEJ3tr5q48zKdL9nIjM9fcoQlxT5OkWsO1auLGnDGPEBLsy59Hr/D/IqM4cjrF3GEJcc+SpFoLONpb8/bw9owa1IbM7DwmfL2LxevjyNOXX3VeCGF6klRrCZVKxaMdGjL7rUdoVr8Oq6JO8/bcnSSkZJg7NCHuKRadVGNjYxk4cCA9e/Zk4MCBnDlzptzzFUVh5MiRvP/++3cpQsvj6+HIjDe6MLB7c+IvXWPU59Fs3nNeJrGEuEssNqnm5uYSERHBqFGj+PXXXxkxYgQRERHlJoevv/6a/fvLvy/8vUCrUfNcz0CmjuiMg50VX/54iOnf7yMjSyaxhKhuFrtO9ejRo2g0GuNdF7t3786UKVM4cuQIwcHBJc7fuXMne/bsYfDgwaSklD9Rk5ycXOo58fHxpgneQtzfzJ05Y7sxd+Uh/jicwMnzVxnZLxCDoi62uUCXo0c2VAlhGmZPqtu3b+f1118vcTw8PBwfH59ix7y8vEhISCiRVC9evMiMGTNYvHgxS5cuve01V6xYwdy5c6sWeA3hZG/NO889wOa/LvDfn48y6bv9qLgPBZVxG2x4ZAwTX+5ISz+5MZ4QVWX2pNq1a1eOHTtW4vj69evZs2dPieNqdfERC51Ox5tvvsmkSZNwc3Or0DUHDRpU4p7hUNBSHTduXAUjrx7VUVZQpVLx+IONaFLPmTGzdpCvFPwbKqgAjEWwF30Yhp2N2X8lhKjRLPYT5OvrS1JSUrFjycnJJVqve/fuJTU1lUmTJgGQmpqKXq9Hp9MRGRlZ6nN7enri6Vn5YiQ13dnLNyitsmDRItiPdmh09wMTohax2KQaFBRETk4Ou3btolOnTkRHR6NSqWjVqlWx87p06UJ0dLTx+zlz5pCYmCi30S2FFMEWovpZ7Oy/Vqvlq6++Yvbs2fTq1Yu5c+cyZ84cNBoNAK+88gpbt241c5Q1y+2KYJf1MyFExakUWcBoFBcXR79+/Vi9enWJFnFtkJ2j5/nJv5Gly4Nb46lFqYBnHvVn8OMBaDUW+/+tEBZNPjn3kMIi2BoMFExT5QMK9rZaxjzbjib1XVix5RTvzI3hSmqmucMVokaSpHqPKSyC3dj2It7WyTS2vcj8MSE80rY+M994mP7dmnHq4lVGfR7F1r0XZCeWEJVksRNVovqUVQTbSqvm+Sdb0SbAk8+XHmDW8oPsP5FMeP8gHO2tK/TccrsWca+TlqooIbi5B3PGdqPj/T7sPHSZNyKjiY1PNXdYQtQIklRFqZwdrHn3Pw8QMTCYm1m5vP/VH3y/8Th6g5QTFKI8klRFmVQqFWEPNWbWm13xq+fCj1tOMX7uThJSpZygEGWRpCpuq76nk3ES6++L1xj9eTRb/pJJLCFKI0lVVEjhJNbHr3XCzsaK2SsOMkPKCQpRgiRVUSlFJ7FiDifwRmQ0R2USSwgjSaqi0v6ZxGptnMT6v1+PySSWEMg61XuOqUoLFkxiNaJVE1cif9jPyq1/c+hUihTBFvc8aamKKqnv6cSMW5NYpy9dY9y83Ry+eR/ndA1IzPXknK4B4ZExHD+bfvsnE6IWkKQqqqxwEmvCiw+i1+eTjxpQodz6WlgEOztHb+5Qhah2klSFyVy/mUNpi6yKFsEWoraTpCpMprAIdmmkCLa4V0hSFSZzuyLY1lby6yZqP/ktFyYT0roe9rZaKHUQAFZsOcVvu8/JTixRq0lSFSZTXhHs156+H2cHG+auPMwni//iekaOucMVolrIOlVhUoVFsIuuU5045jncXF14pG195q86ws5Dlzl1IYrRg9vSJuDeu6utqN2kpSpMrrAIdj3bRNyt041FsB3trRk3rB1vDmlLdo6BD//7JwvXxZKnN5g5YiFMR5KquKtUKhWh7Rvw5ZhHaNGoLj9vj+etWTs4n3jD3KEJYRLS/Rdm4e3mwLSRIfy45RTLN5/krS+288JTrejV2Q+VqvRlWSC3axGWr0ot1SNHjjBw4EAGDx7M9u3bjcdHjhxZ5cBE7afRqBkS1oLpEV2o62zLN2uOMnnhHq7e1Jk7NCHuWJVaqtOmTeOTTz5Bq9Xy8ccfk5KSwoABA7hxQ7pyouJaNHblyzGP8M2ao2zbd5E3Poti1KA2PBDobe7QRC1zN3o6VWqpajQamjdvjp+fH9988w0xMTEsXry43O6bEKWxt7XizSFteXtYe/T6fCYv3MPXq4+QkyeTWKJmqVJS1ev1ZGZmAmBlZcXnn39ObGwsR48eNUlw4t7TpU09vhzbjVZN3Njwx1ne/GI7ZxOumzssISqsSt3/CRMmkJ2dbWw6q9VqZs6cycaNG00SnKh5TFGv1bOuPVNf78zqqL/5YdMJ3pq1g//0aknvLk1NFKUQ1adKSbVVq1YljqlUKnr27FmVpxUCjVrFwO7+tPb34LP/7Wfhujj2H0/m1T4tpAi2sGgmWVJ148YNTpw4wYkTJ3juuedM8ZRCANC8QV1mv/UIC9bF8tvu84z5Mp2c3Fbko0aFgoKK8MgYJr7ckZZ+ruYOV1RRbVgyV6mkqigK586dMybQkydPcvLkSRITE9FoNPj5+UlSFSZna6MlYmBr7m/qxmc/HAA0QEF1AcBYBHvRh2HY2cjSa1E+gyqfG846cq0MWOdp0OlzcMB0SbvCv4GzZs1iyZIl6HQ6FEXBycmJPn360L17d1q2bIm/vz/W1tYmC0yIf8vNK/3GgkWLYD/aodFdjkrUJH+nn+VUsxTy1f9UShuz5WPe7RpBgLtpxuwrPPu/aNEixo8fz7Zt25g3bx62trYkJiYSGhrKfffdJwlVVDspgi2qQpenY9ZfCwsSqgrjn2x9Dp/umIcuzzSbTircUu3bty+DBw8GwMfHhwceeIBPPvmEXr168eGHH8rklKh2tyuCLV3/6leTxzx3nv+LLL0O/vX/soJCVl42f148QLcmnap8nQr/Fk6aNKnY987OzkybNo3t27fz0UcfsWnTJiZOnIirq0wWiOoR0roeC9bFkqXLo8QnA1i++SR1nW3o1q6BbECpwao65pmjzyXhZhKXrl/h0o1//ly5mVzmYzQqDUmZqaYIv+qz/127dmX9+vVMmzaNXr168eeff5oiLiFKKCyCPXnBbgxojLP/9rZWDH28Bauj/+aLZQfZeyyJ8AHBONnLkFRNU5kxz+w8HZdvJP6TOG8l0ZTMdJQid59Qq9R4O3rgV6cBZ69dLPW6BsWAl4O7SV6DSfpLjo6OTJkyhV69epni6YQoU3lFsLu1b8DclYeIOZzA8XPpvDm4LcH+HuYOWVRQiTHPW7L1OUzdPoch9/chKTOVyzeucOl6ImnZV4s9XqvW4uPkyUMN2lLf2Zv6Lj7Ud/bB29EDK40VujwdI9a9S1Ze8SEAFSrsrGzp2KCtSV6HSQehOnbsaMqnIzY2lkmTJpGZmYmDgwPTp0+nSZMmJc7Lz89n/vz5bNmyBZ1OR6dOnZgwYQJqtZSLrY0Ki2AXKiyC7exgzbv/eYCtey/w35+PMuGbXfTt2pTnerbESqsxV7h3VU0d81QUhS1nYsoc89Tpc1h08EcArDVW1HPypqVn84Lk6exDfRcfvBzc0ajLfp9trWwZ3eFlpu2cV6wlbGdlw7sPj8TWytYkr8ViR/Zzc3OJiIhgypQphISEsHXrViIiItiwYUOJ8bLvv/+emJgYli1bhkaj4bnnnmP16tUMGDDATNELc1GpVDzaoRGtmrjz+dL9/Lw9nkOnUhj7bDsa+TibO7xaoSpjnvn5+aRkpRXptidy+dafrLzsMh+nQkXHBm0ZGtQXdwdX1Ko7azA1d22M/2mPYvF/8MJ43FxMNxdksUn16NGjaDQaQkJCAOjevTtTpkzhyJEjBAcHFzt3zZo1jBo1Cjs7OwBmz54trdR7nI/7rSLYW/9m+eaTvDlrO8/3CuTJkCaoy1iWJW6vomOeeoOeKxnJJZJnws0k8gx5xZ7T2caRRnXqA3A85e9Sr6ugEOwdiKdj1cc9NYqautftjd/bam2q/JxFmT2pbt++nddff73E8fDwcHx8fIod8/LyIiEhoURSPXv2LGfPnmXJkiWkpaXx8MMPM3r06DKvmZycTEpKSonj8fHxd/YihEXSaNQMeTyAtgEeRC49wLdrY9l7PInRg9vg5mJn7vBqnPLGPD+O/pLHm3UhMaNgzDMpI5V8pfhmDXd7VwI9mlOvsMvu7E09Z2+cbByNz383xjyrm9mTateuXTl27FiJ4+vXr2fPnj0ljpfWAtXr9Rw6dIivv/4ag8FAREQECxYsKDVZA6xYsYK5c+dWPXhRIwQ0ci2oH7A2lt/3nOeNz6KIGNiaTkG+Jc6t7jHJmjjmqSgK13U3WH9qa5ljnrmGXNaf3GqcaW/ne3+x5Onr5HXbMcu7NeZZ3cyeVMvi6+tLUlJSsWPJycklWq9Q0IJ98sknsbUt+EcPCwtj27ZtZT73oEGDCA0NLXE8Pj6ecePGVTFyUZ3utLSgnY2WN55pTfuWXsz58RCfLtnLow805JW+92Fva1UNkZrPnY55GvINJGemcfnGFS7fSOLyzVtd9huJZJYz3gkFy5a6N+nMC22eQau587RyN8Y8q5vFJtWgoCBycnLYtWsXnTp1Ijo6GpVKVWq5wbCwMH7++WdCQ0NRFIWoqCiCgoLKfG5PT088PeV+8/eijvf7ENCoLrOXH2TL3gvEnkllzNB2tGhccz605anImKcuT0fCzaRbifNWAr2RyJWMZAz5xe+04GzjSMM69ann5EVWXja7Lu4v9br5Sj7NXBtXKaEWqu4xz+pmsUlVq9Xy1VdfMWnSJKZOnYqdnR1z5sxBoylYMvHKK68wePBgunfvzujRo4mMjKR3794YDAYefPBBXn31VTO/AmGpXJ1tmfjKQ2z44yyLfolj/NydPPNoAIMe8zd3aFVS1phnll7HpKhZBLg1JTEzmbSs4us7VSoVng7uBHu1pN6tcc56t7rsheOdhc9/6EpcjR7zNEUR9dux2KQKEBgYyIoVK0r92bfffmv8u42NDe+9997dCkvUAiqViidDmnB/M3cif9jP8s0nOXgymed7NC5RBDv92k2TjnmasvRcniGPxIwUEm4msfNc6XvbAfT5ek6knqahiy8t3Jv+kzydvPF28sRac/shkNoy5lndLDqpClHdGnk7EznqYX7YdILV0af5cME19Ib7UFAZt8GOnbePia90MkkR7DspPacoCjdybhq77Ak3k0i4tTwpKTMVRSm9yExRGpWa3i0eY0hQnyrFXxvGPKubJFVxz7PSanj+yVa0aurGxwv2oNyqiGksgp1rMEkR7PKWJH26Yx5zn5zCdd0N43rOhFsJ9PLNRDJzi5c11Kq1eDt68EC9YHydvKjn5E3CzSTWHN9U6rUNSj7ejqbZslvTxzyrmyRVIW65diOH0tp8piqCvevigXJLz720ZmyxQiBwa6LIxRdfp4IxznrOXvg6eeHh4FZiS6YuT8dvf0fLmKeZSVIVNUp1rvMsLIJtIA+NayIqmyyUHHsM6d5osKpUEexcfa5xrLPwz9GkE+U+xtvJgw71WhsniXydvHC0qfjrkjFPyyBJVYhbvF3tUezTsfXfDxo9KCpQKVg1PEHuqXYlxi7zlXzSsq5y5WayMXFeudVtT826WqLVaactP6n1bRFW5SLJMuZpfpJUhbil/X3u2Jzej6LRo1IBqoKkqGj0WAfsY9XhmxzP+RNXDwNJGclcyUgm91/72K01Vvg4efGQmx8+Tp7GFqePkycalfqubMOszjHP2tA9r26SVIWgoLu+9Ww0aPUlViSpVIDGgHWTY5zOBS6Bq50rLT2aGxOmr5MXvs5euNrVKbeCknTPaz9JqsKk7sbe9jtd56nPN5CSmcaVm8lcuZlU8DUjmSs3CxbE/7u7XpQKFQ/Ua439tUB+255MstqKJ3q2pHfrppWqeiXd89pPkqqoUW63zjNfySc9+1qRxJlS8DUjmeSMVAz/qpxko7XBx9GD5m5+5OhzOHAlttTrKii0872PbiGdCL0/nc+XHmDhujj+iiuoeuXpal/q40ojS5JqN0mqosYobxvmxG1f4OPkSVJmaol6nVq1Fi9Hd9r43o+Powc+Tl74Onni7eRJXVsXY9FzXZ6O19a9S/ZtxjxbNHLly7ce4bv1cWzcdY43IqN47en7LeKGgzLmaX6SVO8xltw9L6qw3FxiRgpXbiaTmJHC4cRjZW7DNCgGbuhuEujR3DjG6e3oia+TJ+72rhUqWm5rZUt48DC+KEzcxuMlxzxtbbSE9w/mwVbefLniIF8sO8ju2ERGDgjGxVFanvcySarCpCqzDbMwcV7JSCbxZkrB14wUEm8lUZ0+p8LX1ag0PNqsC4Pv712l+O/zaVGpMc92LbyYMzaUr1YdNt5w8I1nWtMh0LtKcYiaS5KqMJnb3Q3z2aCnScu+SuLNFBIzSk+cKlS42deluVtjvBw98XH0xNvJAx9HT46n/M23+5eVem1T3mK4smOezg7WvD28PQ/ed5mvVx/h44V7ePzBRrzUu1WJWq3SPa/9JKneg0xZJQkKbuaWmn21YItkOXfDXHhgOVCQON3t69LczQ9vRw98nDzxdvTA28kTTwf3MismudvX5YfDayxyG6ZKpeKRtvW5r4kbs5cf5Pc95zlyOoU3h7Ql0M/NeF5NrPwvKkeSqoWp7g/dnVRJgoJ1nMmZaSRmpJCUkWL8mpSRSnJWWonixv+mQsVDDdoy8L5eeDm4Y1WBUnP/VhO2YbrXsWPSqx3Z8MdZFq+P4915MfTr1pyhYS2w0srNKO8FklTvIberkvTZExO4rrtJYkYySRmptxJnKkkZKaRnXyvxfFZqLV6OHrT2DsTL0YMbupvEXNhb6rUVFFp7B1LfueTtcCqjutd5mqJ7rlareKpLE1r7e/D5sgP8tO1v9p9IYszQdrg7a0rUa9Xl6JGGau0hSdUCmbp7Xijmwr5yqySF//J+icc4WNnh7ehJC/emeDl64OXogbejO16OHtS1cym2e0iXp+NAwtEavQ3TlBp4OTHzjS78uOUUK7acYvQX2+nW1pcjNwMxoDHWaw2PjGHiyx1NUq9VmJ8k1TtQnV30O+2eQ8Fs+s3cTJIzUknOTCU5M62ge56ZajxWnkZ16tGxQTu8HN3xcvDA29FDqiRVkVajZmhYC9q39OKzH/axee8lQAOo/qnXmqM3Sb1WYRnkHbQgt+uef/3UJ6hVapIz00jOTL2VMNNIykwl5dbfs/W6Es9rq7XB08GdRnXqc+7apTKv37N5qMVXSaqps+f+DevSt2szvlp1hH93FUxVr1VYBkmqd6g6uui3K2L82rp3S02aapUad/u6NHNrhIeDO14O7ng6uuHl4IGngxtONo6oVCp0eboaXyWpJku7riuo15pfssaARq2qVL1WYbkkqd6BqnTRcw15pGalk5KZRkrmra9Z6aRmpnH+2uVyH2untaGt7314Orjj6eCGl6M7ng7uuNnXLVEFvjTSPTcvb1f7UhMqgCFfwb2O3V2OSFQHSaqVdLsu+qweH3EzN/OfpJlV8DU1M43krHSu626U+rxONo44WjuU2hItNOj+3tI9r8FCWtdjwbpYsnR5lLbX9uft8TTwcqJVE7eSDxY1hiTVSrpdF/3Vde+U+ri6ti54ObjTytMfD3tXPBzc8HAo+Opu74qt1ka657WcnY2W8cNaM3nB7mKz//a2VoQE+xK1/xLvzIsh7KFGPN8rEEd7a3OHLO6AJNVKSs5MQaNSlyghV6ihiy/t6wXjeStZeji44WZfV+6rLgAIaFiHIKdjxdapThzzHG6uLvTr1pz5Px3mt93n2ROXyCt97qNL63pmr3wlKkeSaiV5OniUmVABevl3r1IXXbrntZ9GlY+7dbrxe9tby6jqeTgyZUQntu27yMJ1ccz833627rvI6/2C8HaT3QE1heybq6RODdpir7Xl30XiVaiwt7IzSRe9sHvulepE3ev20j2/h6hUKro/0JCvxocS2r4BB04kM3JmFKu2/Y3eUPZ/5sJySEu1kqSLLqqioj0FF0cb3hzSltD2DZj/02EWbzhG9IFLRAwMJqCR7LyyZCrl3/fdvYfFxcXRr18/Vq9eTatWrco8LzMzk+mRM+Q+Q+KuyM0z8OPWU6za9jeGfIWenfwY3qMlDnaVL0ojqp+0VO+QzKCLu8XaSsOwJ1rycOt6zPvpMBv+OMufR6/w6tP30+l+H5nIsjAypipEDdHQ25lPw0OIGNianDwD05bsZcp3f5F8VXZiWRJpqd4BmUEX5qJWqwh7qBEdWnmxcG0c2w9e4sjpFIb1aMmTnf3Q6bKlCLaZSUtViBqorpMtY4e1Y9IrHXFxtGHB2ljGfLmDMwk3MChqUnNduazzJjXXFV2O3tzh3lNkoqqIik5UCWFJdLl6lv9+kjXb48nPV1CRj4Kq2I4tqdd690hLVYgaztZay/NPtmL6yM6oVaCgpqBea8HXwnqt2dJivSskqQpRS1xMyqC0IlhF67WK6idJVYhaIjE9C4269OVVKuDUhWt3NZ57lUXP/sfGxjJp0iQyMzNxcHBg+vTpNGnSpMR5mZmZvP/++5w+fZr8/Hx69erFyJEjzRCxEOZTXr1WBdj05zkM+QpDwwJwc5HardXFYluqubm5REREMGrUKH799VdGjBhBREQEpc2rfffddwD88ssvrFy5klWrVrF79+67HbIQZhXSuh72tlr+XZhCpQJbGw0BjV35fc95Xv1kC4vXx5GRlWueQGs5i02qR48eRaPREBISAkD37t3Jzs7myJEjJc41GAxkZWWh1+vJycnBYDBgayt78MW9pbBeqwYDULAKABTsbLRMfqUTMyJC+Ojlh/D1cGRV1Gle+WQLq6P+JifPYO7QaxWzd/+3b9/O66+/XuJ4eHg4Pj7F7xHv5eVFQkICwcHBxY6/8sorDB8+nC5dupCdnU3v3r1p3bp1mddMTk4mJSWlxPH4+Pg7exFCWIi2Leux7JM+/HH4MolpWXi72dM5uJ7xLq3tW3rRJsCT7Qcu8cOm4yxaf4xfdp5haFgLQts3QKOx2HZWjWH2pNq1a1eOHTtW4vj69evZs2dPieNqdck3fdKkSbRv357x48dz8+ZNXnvtNRYvXszzzz9f6jVXrFjB3Llzqxy7EJbIzkZb7l1ZNWoVoe0b0KW1L7/uOseKzaf48sdDrNkez3M9W/JgK2+pJ1AFZk+qZfH19SUpKanYseTk5BKtV4DNmzezcuVKNBoNderUoW/fvkRHR5eZVAcNGkRoaGiJ4/Hx8YwbN84k8Qth6ay0Gvo83JRHH2jImujT/LwjnqmL/qJlY1f+0ytQ7pV1hyw2qQYFBZGTk8OuXbvo1KkT0dHRqFSqUnc6tWrVivXr1zN69Ghyc3OJioqiTZs2ZT63p6cnnp6e1Rm+EDWGg50Vw3q0pGdnP5ZvPsnvu8/zzrwYOgR681zPljTycTZ3iDWKRW9TPXbsGJMmTSIjIwM7OzsmT55MYGAgUDCOOnjwYLp3705CQgKTJ0/m3LlzaLVaOnfuzNixY7Gyqly9SdmmKgQkpGTw/cbjxBxOQKWC0PYNGBrWAs+6BaUus3P0xBy6TGJ6Fj7/GrMVFp5U7zZJqkL84++LV1m8/hhHTqdipVXTq7MfQc08iPxhH5k6PRq1CkO+goOtlo+ktoCRJNUiJKkKUZyiKBw8lcKSDcc4c/l64VGK3kNdpQJ7Gy2LPgyTFisWvE5VCGF+KpWKtgGefDG6K088VLiioPjKAKktUJwkVSHEbanVKpwdbcqsLaBRq0hMkzsQgCRVIUQFlVdbwJCvsOPgZaL3XyT3Ht+hJUlVCFEhZdYWALQaFSnXsohceoDnJ//Ool/iSEjNMEuc5iajykKICrGz0TLx5Y5MWvBnsdl/+1uz//U8Hdm27wIbd51jdfRpVkefpo2/Bz06NaZDoPc9swVWZv+LkNl/IW4vO0dfZm0BgPx8haOnU/n1z7Psjk0kP1/B1dmWsIca8fiDjXCvU7vLDkpLVQhRKberLaBWqwj29yDY34O069ls/usCv/15jmW/n2TFllN0CPSiR0c/Wvt7oC5j4qsmk6QqhKg2bi52DH4sgIGhzdl3PImNf55jT1wiu2MT8Xaz54mHGvNoh4a4ONoYH1PTd2xJ978I6f4LUf0S0zL5bfd5Nv91nusZuWg1ajoH+dKjU2N0umw+XrgHA5oaezdYSapFSFIV4u7J0xv48+gVNv55jtj4NKBgd1ZBSqq5O7bujek4IYTFsdJqeLhNfT4ND2HeuG4EN3enoIlXs3dsWX7aF0LUeg29nQlo5EpsfFqZGwy++yWOg6dSaFqvDs0auNCkXh0c7SpXie5ukKQqhLAI5e3YgoJVBTsOXmbHwX9arD5uDjSp70LTei40q1+HpvXr4OxgXe51qnsiTJKqEMIihLSux4J1sWTp8ihtTPXb9x4jP1/hzOXrxF++Rvyl65y+dI1dRxL443CC8XzPunY0rV+HpvVcCr7Wd6GuU8GNQA8cu8zkhbuLTYR9uzbWpBNhMlFVhExUCWFex8+ml9ixdbt6rdk5es4mFCTY+EvXOXP5OheSbpJfpNXr6myLn68zR06nkqc3UJ0TYdJSFUJYjJZ+riz6MKzcHVv/ZmejJdDPjUC/f+6ppcvVc/7KDeIvX+f0xWvEX77OwZPJ5N9mIqy8TQ0VJUlVCGFRbrdjqyJsrbUENHIloNE/rdslG+JYHR1frAVbyJSlC2VJlRDinuDr7lhqQoWC0oXebvYmuY4kVSHEPaHM0oUqcLDV0jm4nkmuI0lVCHFPKCxd6GBbsLa18C4G9jYFE2GmWlYlY6pCiHvGnUyEVZYkVSHEPcUUE2Hlke6/EEKYkCRVIYQwIUmqQghhQpJUhRDChGSiqoicnBwA4uPjzRyJEMISNWnSBDu78m9cKEm1iEuXLgEwbtw4M0cihLBEFSm2JFWqikhPTycmJob69etjY2Nz2/Pj4+MZN24cM2fOpGnTpnchQnG3yHtbu93p+yst1UpydXWld+/elX5c06ZNpVRgLSXvbe1WHe+vTFQJIYQJSVIVQggTkqQqhBAmJElVCCFMSJJqFXh4eBAREYGHh4e5QxEmJu9t7Vad768sqRJCCBOSlqoQQpiQJFUhhDAhSaq3sXPnTvr06cMTTzzBCy+8QHJycolzrl69yogRI+jZsydhYWFs2bLFDJGKO1GR93fv3r20bduWPn36GP+cPn3aDNGKOzVr1izef//9Un9m8s+vIsqUlpamdOjQQTlx4oSiKIqyZMkS5cUXXyxx3qhRo5TIyEhFURTl4sWLSqdOnZSEhIS7GquovIq+v/Pnz1emTZt2t8MTJnDx4kUlPDxcCQoKUt57771SzzH151daquWIiYkhICCAgIAAAAYPHsyePXtISUkxnqPX64mKimLw4MEA1K9fn5CQEH755RezxCwqriLvL8D+/fuJjY2lf//+DBgwgE2bNpkjXHEHVqxYQadOnXjhhRdK/Xl1fH5l7385EhMT8fHxMX5vbW1N3bp1uXLlinEpxtWrV9HpdHh7exvP8/b25sqVK3c9XlE5FXl/AZycnHjqqaeM3f7hw4fj5eVFmzZtzBG2qIQxY8YAMGfOnFJ/Xh2fX0mq5VAUBZVKVeK4Wq0udg5Q4ryi5wjLVJH3F+CLL74w/r1Zs2b07NmTLVu2SFKtBarj8yuf/HL4+vqSlJRk/D43N5erV6/i6+trPObm5oaNjU2xCY6kpKRiLSBhmSry/mZmZjJ//nz0er3xmKIoaLXSHqkNquPzK0m1HJ07d+bYsWOcOnUKgJUrVxIcHIyrq6vxHI1GQ/fu3Vm6dCkAly9fZufOnTz66KNmiVlUXEXeX3t7e1atWsXPP/8MwMWLF9m0aRM9evQwR8jCxKrl83vHU1z3iD/++EPp27ev0qNHD2XYsGHKpUuXFEVRlN69eytHjhxRFKVgFnnkyJFKr169lLCwMGXdunXmDFlUQkXe3+PHjytDhw5VnnzySaVHjx7Khg0bzBmyuANffvllsdn/6vz8yjZVIYQwIen+CyGECUlSFUIIE5KkKoQQJiRJVQghTEiSqhBCmJAkVSGEMCFJqkIIYUKSVIUQwoQkqQohhAlJUhXiXzZv3ky/fv1o27Yt/fv3JzExkVdffZUNGzaYOzRRA8g2VSGKWLt2LdOmTeObb76hRYsWvPzyy1hZWZGRkcGKFSvMHZ6oAaSlKsQtBoOBmTNnMm7cOIKCgrC2tqZ169bExMTwzjvvmDs8UUNIUhXilkOHDpGenl6srF9ubi6PP/64FKQWFSZJVYhbEhMTcXFxwc7ODoBr166xZs0amjdvbubIRE0iSVWIW3x8fEhPT2fv3r1kZWXxzjvvYGVlRUJCgrlDEzWIJFUhbmnTpg0vvfQSo0ePJiQkhEaNGjF37lyio6OZPn26ucMTNYTM/gshhAlJS1UIIUxIkqoQQpiQJFUhhDAhSapCCGFCklSFEMKEJKkKIYQJSVIVQggTkqQqhBAmJElVCCFMSJKqEEKYkCRVIYQwIUmqQghhQv8f/8LWMxkXRPAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 340x235 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.4, 2.35))\n",
    "\n",
    "ax.errorbar(\n",
    "    alphas,\n",
    "    means_HG,\n",
    "    yerr=std_HG,\n",
    "    fmt=\"-o\",\n",
    "    color=\"C0\",\n",
    "    ecolor=\"gray\",\n",
    "    elinewidth=3,\n",
    "    capsize=0,\n",
    "    label=\"random hypergraph\",\n",
    ")\n",
    "\n",
    "ax.errorbar(\n",
    "    alphas,\n",
    "    means_MSC,\n",
    "    yerr=std_MSC,\n",
    "    fmt=\"-o\",\n",
    "    color=\"C2\",\n",
    "    ecolor=\"gray\",\n",
    "    elinewidth=3,\n",
    "    capsize=0,\n",
    "    label=\"random flag complex\",\n",
    ")\n",
    "\n",
    "ax.set_ylabel(r\"$\\lambda_2$\")\n",
    "ax.set_xlabel(r\"$\\alpha$\")\n",
    "\n",
    "ax.set_xticks([0, 0.5, 1])\n",
    "\n",
    "sb.despine()\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "fig_name = f\"figure_1_lambda2_HG_SC_N_{N}_ps_{ps}_nrep_{n_repetitions}\"\n",
    "plt.savefig(f\"{fig_dir}{fig_name}.pdf\", dpi=250, bbox_inches=\"tight\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6532c771",
   "metadata": {},
   "source": [
    "## Fig. 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "952fe950",
   "metadata": {},
   "source": [
    "### Generate random structures and compute exponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6a10eb47",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 50  # number of nodes\n",
    "\n",
    "p_1s = [0.2, 0.4, 0.6, 0.8]  # wiring probability of 1-hyperedges\n",
    "p_2 = 0.05  # wiring probability of 2-hyperedges\n",
    "alphas = np.arange(0, 1.01, 0.05)\n",
    "\n",
    "HGs_2 = []\n",
    "lyaps_HG_2 = np.zeros((len(p_1s), len(alphas), N))\n",
    "\n",
    "for j, p_1 in enumerate(p_1s):\n",
    "\n",
    "    ps = [\n",
    "        p_1,\n",
    "        p_2,\n",
    "    ]  # ps[i] is the wiring probability of any i+2 nodes (ps[0] is for edges, e.g.)\n",
    "\n",
    "    # generate hyperedges\n",
    "    HG = xgi.random_hypergraph(N, ps, seed=0)\n",
    "    HGs_2.append(HG)\n",
    "\n",
    "    # compute exponents\n",
    "    for i, alpha in enumerate(alphas):\n",
    "        lyap_multi = compute_eigenvalues_multi(\n",
    "            HG, orders=[1, 2], weights=[1 - alpha, alpha]\n",
    "        )\n",
    "\n",
    "        lyap_multi = np.sort(lyap_multi)[::-1]\n",
    "        lyaps_HG_2[j, i, :] = lyap_multi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "22f40a59",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bound_multi(H, alpha, rescale_per_node=True):\n",
    "    \"\"\"Returns the lower bound N/(N-1) k_min\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    HG : xgi.HyperGraph\n",
    "        Hypergraph\n",
    "    orders : list of int\n",
    "        Orders of interactions to consider.\n",
    "    weights: list of float\n",
    "        Weight of each order, i.e coupling strenghts gamma_i in [1]_.\n",
    "    rescale_per_node: bool, (default=False)\n",
    "        Whether to rescale each Laplacian of order d by d (per node).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    float, bound of the Lyapunov exponent\n",
    "    \"\"\"\n",
    "    L_multi = xgi.multiorder_laplacian(\n",
    "        H, orders=[1, 2], weights=[1 - alpha, alpha], rescale_per_node=rescale_per_node\n",
    "    )\n",
    "    K_multi = np.diag(L_multi)\n",
    "    N = H.num_nodes\n",
    "    return -(N / (N - 1)) * np.min(K_multi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c1c1df22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute theoretical bound\n",
    "bound = np.array([bound_multi(HGs_2[-1], alpha) for alpha in alphas])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0446955",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b04ab8c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAEKCAYAAAAlye1PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABtf0lEQVR4nO3dd3hUVfrA8e+dSe+9EiC0hBpCJ4SOdAREiiKCioAIu64uKywIgoiUVVRY+WFZBUVApUkvkd5bgFACBEJID+m9TO7vj5CBkJkhfVLO53l4kty5c+cMmcw755z3vEeSZVlGEARBEIQqp9B3AwRBEAShrhJBWBAEQRD0RARhQRAEQdATEYQFQRAEQU9EEBYEQRAEPRFBWBAEQRD0RARhQRAEQdATEYQfy8zM5Pr162RmZuq7KYIgCEIdIYLwY/fu3eOll17i3r17+m6KIAiCUEeIICwIgiAIeiKCsCAIgiDoiQjCgiAIgqAnIggLgiAIgp4Y6LsBNV1Mci6nglOJT83F3tIQPy9LnK0N9d0soY7LzMzk3r17NGrUCFNTU303RxAELUQQLoeTwan8fCwOCZABCdh/JYnXezji52Wp59YJdVlhtv/WrVtp2bKlvpsjVJDwqDj2/HWWqLhEXB1tGdynM/VcHfXdLKEcRBAuo5jkXH4+FocsFwRgePJ1/bE4mriY4CR6xIIgVJA9h8+xYs1mkCSQZZAkNu44zL/eGcug3p303TyhjMSccBmdCk5F0nKbREEvWRAEoSKER8WxYs1m8mWZ/Pz8Il+Xr9lMeFScvpsolJEIwmUUn5qr7vk+S358uyAIQkXY89fZgh6wJpJUcLtQI4kgXEb2loY6e8L2lmIoWhCEihEVl4gsa/nYL8tExSVWbYOECiOCcBn5eVlq7wnL0E0kZgmCUEGc7W20B2FJwtXRtmobJFQYEYTLyNnakNd7OCJJoJAKRooKR4skCcLjs/XbQEEQagVZlol+lKDrBAb36Vx1DRIqlMiOLgc/L0uauJhw8ql1wp5Oxmw4/ohvD8Uyxk9Fn1bW+m6mIAg1zNNLkZKS07gUdIeG9ZwJi4hFkqSCXvHjLOl/vTNWLFOqwUQQLicna0NGdrIrcqyevRFf74lm86l4EtPyGNnZDoW2pAqhVgkKCmLhwoWkp6djbm7OsmXLaNSoUbHz0tPTmTt3Lnfv3iU/P58hQ4bw7rvvAvDw4UP+/e9/k5CQgEKhYNGiRfj6+lb1UxH05OmlSLIsq4ehRwzwp6NPM7FOuJYRQbgSOFga8q/hbnyzP5oDV5NJTM9jYi8nDJUiENdmOTk5zJgxg8WLF+Pv709AQAAzZsxg9+7dSM98CPvf//4HwM6dO8nIyGDYsGG0b9+eLl268P777/Pyyy8zduxYgoKCmDZtGgcOHMDMzEwfT0uoQk8vReKZOeCv/7eVn7+azZTxQ/XUOqEyiDnhSmJhouS9Ia60bWjG+ZB0Vu2NIjMnX9/NEirRtWvXUCqV+Pv7A9C3b18yMzO5evVqsXNVKhUZGRnk5eWRnZ2NSqXCxMSEmJgYbt26xUsvvQRAq1ataNiwIUeOHNH4mLGxsVy/fr3Yv5CQkEp7nkLlEUuR6h7RE65ERgYKpvZzZvOpeI7cSGHFn5HMHOSCrbn4b6/Jjh49yjvvvFPs+PTp03F1dS1yzNnZmcjISHx8fIocf/vtt5kwYQLdu3cnMzOTF198kbZt2xIYGIiDgwOGhobFrqHJ5s2bWb16dQU8K6E6EEuR6h4RDSqZQiExrps9thYGbDuXwLLtEfxtkCtKpSQ2fqihevbsyY0bN4od37VrF2fPFu+pKBTFB5wWLlxIhw4d+PDDD0lNTWXq1Kn89NNP+Pj4FBu61nYNgLFjx9KnT59ix0NCQpg1a1ZJno5QjTg72IqlSHWMCMJVQJIkBra1wcZcyfqjcXy2LYJclVyQ3IjY+KG2cHNzIyYmpsix2NjYYr1jgIMHD/L777+jVCqxsbFhxIgRHDlyhEGDBvHo0SPy8vIwMDBQX6Nfv34aH9PJyQknJ6eKfzJClXh2Q4awyFjtJ4ulSLWSCMJVqEtTS/LzZdYdfQQ8ybsQGz/UDm3atCE7O5tTp07h5+fHkSNHkCRJ4y5GLVu2ZNeuXbz33nvk5ORw+PBhfH19cXZ2xtvbm23btjF69Ghu3LjBnTt38PPzK3f7xA481cuzGzLIFKwJtrOxJDE5DYVYilQniCBcxWKS89SboDyrcOOHZ5c8CTWDgYEBa9asYeHChXz66aeYmpqyatUqlEolUDAPPG7cOPr27cvy5ctZtGgRAwcOxMDAgG7dujF58mQAPv/8cz766CPWr18PwBdffIG1dfnWm4sdeKoXXVnQiclp/GfeFC5duyM+MNUBIghXMV0bO4iNH2q+Fi1asHnzZo23fffdd+rv3dzc+L//+z+N53l4ePDTTz9VWJs0vuE//rp8zWZae3uKN/gqps6C1vBpXJIkLl27I5Yi1RFiiVIVExs/CFVNLHupfqLiEjUPh4HIgq5jRBCuYro2fsiXwbehKMggVCzxhl/9uDraov3TuMiCrktEEK5imjZ+UDz1x/jjkTgeiSFpoQIVvOFr7wmbmxrz7YZdLPzyZ77dsEtsEF8F+nVvT36+9g9GIgu67hBzwnqgaeOHbl6WXAvL4PfT8SzbXlDUo76Dsb6bKtQCg/t0ZuOOwxpvy8/PZ1fA2YK1ySJhq1Kps9NjE7gb+qT4ilKhEFnQdZgIwnqiaeOHvq2tsTVX8sPhOP7zZyRTXnCmlYcYnhbKp56rI/96ZyzL12wusgNPfn5BGdWnNwkQCVuV4+ns9ML/d4Aprw4mPSNLZEHXYSIIVzPtGllgZWbAN/uj+e++aF7r7kA3byt9N0uo4Qb17kRrb88i64STUzPYHXBGc47C44QtkaFbfrqWI32/cS8/fzVbBN46TMwJV0NNXEz413A3bC0MWH/sETsv6qgnKwglVM/VkSnjh7LgvQlMGT+UjKxsjSUyAZGwVYFEdrqgi+gJV1MuNkZ8ONyN1fui2XUxkcS0PMZ3d+BRap6oOS1UCHXCloYPeDIQHZvAwi9/FsOk5SQ2ZRB0EUG4GrM2M+CDoW58GxDDyeBU7sdmEZWUi4SoOS2Un66ELVmWuRUSxq27YSJZq5yMjQzEpgyCViIIV3MmRgreHeDC9wGxXLqfDjypNS1qTgvloSlhK/+pYKFeQiOStUrl6RrdFmamHD8bpP1ksRypzhNBuAZQKiQcrQzUPeBniZrTQlk9m7AVHZvAzbthmntuIlnruZ7Ogn4667xdqyYEXg8pkp0uliMJIIJwjZGQloe2KCxqTgvlUZiwBbDwy5+5dTdMc8a0mL/USVcWdOD1EFZU0KYMsixrT6irBDV9962UlBQ++ugjjh07hrm5OZMnT2bSpElazw8PD2fRokUEBgaiVCrp3r078+fPx8LColLaJ4JwDVFYc1pbT1jUnBYqgq5kLSTU1bVq6htyZdK1KQMVsCmDLMv89ddfpKSkMGLEiCoJxLVh961FixaRk5PD8ePHiYiIYNKkSXh6etKzZ0+N58+fPx97e3uOHTtGdnY2M2fO5KuvvmLu3LmV0j6xRKmGeF7Nab9mlfMpTahbBvfprLXOdH6+zK6AM2z88whHTgWy8c8jTPj7UvYePlfFrayeKrNGt0qlYvv27Zw4cYL4+HhycnLKfK2Serpnn5+fX+Tr8jWbK7y86f79+xkxYgSzZ8+mZ8+edO3ald9//71c18zIyGDfvn384x//wMLCAi8vL8aMGcOWLVu03ic8PJyhQ4diYmKCtbU1AwYM4Pbt2+Vqhy6iJ1xDFNacXn8srkiPuPBv/vD1FMZ0tUehqLphKqH2eX51LZDlxxWfRMJWES4OtpWSBZ2VlcVvv/3G/fv38fLyYtSoURgaln3k66ff93P4VOBzz0tMTi2SqPe0fFlm+tyvsLXWviqjt19bJo0eUOJ2BQcHc/fuXSZOnMiSJUvYs2cPc+bMYeDAgVhaPnmcqVOncvHiRY3XcHV1ZefOneqfQ0NDkWWZZs2aqY95e3tz4MABre2YOHEiO3fupEOHDmRnZ7Nv3z569epV4udRWiII1yCaak77epqx5UwCh6+nkJiex1t9nDAyEAMcQtlpqq6VkpbJ7r/OaN50oA4nbD1dDzr0YbTW0aqyZkGnpKSwYcMGYmNj6dChA4MGDUKhqJq/77w8VbluL63g4GDGjBnDyJEjARg6dChz584lPDyc5s2bq89bu3Ztia+ZkZFRbC7XysqK9PR0rffp1KkTW7ZsoUOHDuTn5+Pv78+ECRNK+WxKTgThGkZTzem/DXZl/dE4zt1N44tdUbw7wAVLU6WeWijUBk8na0FBwpbWCCPLhIRF1bm54iKZ0Pn56v+egp3RKmZTBqVSSV5eHv369cPPz69C5oEnjR5Qoh7qtxt2sfHPI0VqXRdSKBSMGNCtQj94BQcHFwl2qampZGVlYWdX9lUfZmZmxQJuamoq5ubmGs9XqVRMnjyZUaNGsXHjRnJzc/n000+ZNWsWX331VZnboYsIwrWAoVLijd6O2FkYsC8wiWU7Ipg5yFVU0hIqjK6ErXxZ5sylm5wLDK6xyTulpSsTGhkG9+1MemZWmT+QpKWlYWFhgbm5OdOmTSvX8HNZ6SrmUtHrm9PT0wkPD8fR8cn/0759+2jRogXOzs5Fzp08ebLW4Wg3Nzd2796t/rlhw4YA3Llzh6ZNmwJw8+ZN9ffPSk5OJjo6mtdeew1jY2OMjY155ZVXmDhxYnmenk5i3LKWUEgSIzvZ8ap/QWnL5TsiCInJ0nezhFpCV8JWoapI3qkudNWDlhQKrCxM1TW6SxuAAwMD+eqrr7h79y6AXgIwPMkPUEgSSoUChSShePy1otc33759G4VCwc6dO1GpVJw/f54vv/ySWbNmFTv3+++/5/Llyxr/PR2AoaAnPGDAAFauXElaWhq3b9/mjz/+YNSoURrbYWdnh4eHB7/++is5OTlkZGTw22+/4eXlVWHP9VmiJ1zL9Gxhha25ku8CYlm5K4q3+jjh66l56EUQSup5CVuayMCnq37Fxcmu1g1RF9SD1vLcy5gJLcsyR48e5ejRo1hbW2NtbV3OVpafpvyAyvg9BgcHM2jQIMLCwujQoQPu7u4sXrwYPz+/cl97wYIFzJs3j+7du2Nubs7bb79dZHnS5MmT6dChA9OmTQNg9erVfPbZZ6xfvx5JkvDx8WH58uXlboc2IgjXQm0amPPBsILNH9YejGGMnz0tPczExg/A1atX+eSTT1AqlbzzzjvqP8Z3332X//73v3puXfWm6Q35XlgUZy/f0phFK8syN+48qJX1p3Nz87QPDJQhE1qlUrFr1y4CAwNxdXXllVdeKZIRrE/P5gdUhuDgYLy8vJgyZUqFX9vKyoqvv/5a6+3ff/99kZ+9vb1Zt25dhbdDGxGEa6mGjsbMHu7G13uj2XwqHohHIYmNH5YuXcqSJUswMDDgk08+IS4ujpdffpmUlJQKuX5QUBALFy4kPT0dc3Nzli1bRqNGjYqdl56ezty5c7l79y75+fkMGTKEd999F4Dt27fz2Wef4eLioj7/559/xspK//tKP/uG/O2GXZwtnAvW4ul505q6nOnpqlEqlYqT5yuuHrQsy2zatIm7d+/SpEkTRo8ejZGRUQW0uuYIDg7G399f383QCxGEazEHK0Pe6O3I0u2RQEFRD6jbGz8olUp1UsbatWuZNWsWaWlpFZJ1mpOTw4wZM1i8eDH+/v4EBAQwY8YMdu/eXez6//vf/wDYuXMnGRkZDBs2jPbt29OlSxcuXLjAzJkzee2118rdpsqmM3lHE0li05+HsbIwqzGZ1NrqQXfv1JqT54PKXQ9akiQaN26MpaUlQ4YMQamseysbbt++TZMmTfTdDL0QQbiWCwzNQCE9CcBPq4sbP+Tl5al7qYaGhnzxxRf861//4tq1a+W+9rVr11AqlepP9H379mXx4sVcvXoVHx+fIueqVCoyMjLIy8sjOzsblUqFiYkJABcvXiQqKoqtW7diYmLCe++9R6dOmodwY2NjiYsrnvwUEhJS7udTEs/bielZcn4+Ow+dKVjrWgMyqXVlQZ88H1SuetBxcXFYWFhgampKly5dqrwmdHVy4cIFfTdBb0QQruXiU3O1L++k7m38MG/ePDIzM9XrBBUKBStWrGDv3r0lvsbRo0d55513ih2fPn06rq6uRY45OzsTGRlZLAi//fbbTJgwge7du5OZmcmLL75I27ZtycnJwc3Njbfeegs/Pz/Onz/P9OnT2bp1Kx4eHsUec/PmzaxevbrEba8MmnZiuhUSprGwR+ERdUJXNR+mrqx60Pfv32fz5s24u7vz2muvIUlSnQ3AdZ0IwrWcro0fCm+vS1q2bFnsmCRJDB48uMTX6NmzJzdu3Ch2fNeuXZw9e7bYcU0VjhYuXEiHDh348MMPSU1NZerUqfz0009MmjSJH374QX1ex44dad++PcePH+fVV18tdp2xY8fSp0+fYsdDQkI0Lu+oLE/PFYdHxTHh70tLdf/qmkldkAVdsfWgr127xvbt2zE0NKRbt24i+NZxIgjXcn5eluy/kqTxNlkGa7O6uVQ8JSWFW7ducevWLV5//fUKuaabmxsxMTFFjsXGxhbrHQMcPHiQ33//HaVSiY2NDSNGjODIkSMMGjSI3bt38+abb6rPlWVZ61pRJycnnJycKqT9FUXXciZtncrqkkn97LZ9aemZFVYPWpZlTp48SUBAAJaWlowfP75YIQqh7hFBuJbTtPGDRMEbobGhxG+nElDlQ7/W1rXyE7ksy4SGhqoDbnBwMMHBwURHR6NUKvH09KywINymTRuys7M5deoUfn5+HDlyBEmSNPa+W7Zsya5du3jvvffIycnh8OHD+Pr6YmZmxjfffEPz5s3p2rUr165dIzAwkE8//bRC2lhVtNefPqt9fS3FM6kd7ay5FFT+PXhL4tlt+2TQHoCh1FnQZ8+eJSAgACcnJ8aPH18tst0F/ZNkna+yuuP69eu89NJLbN26VeObZk0Xm5xbZOOHbl6WKBWwal80UYm59GphxRg/e5S1aBemL7/8knXr1pGVlYUsy1haWjJ8+HC8vLxo3rw5zZo1q/ClIDdu3GDhwoWkpaVhamrKokWLaNGiBVAwDzxu3Dj69u1LZGQkixYtIjQ0FAMDA7p168Y///lPDA0NOXv2LMuXLyc7OxsDAwM+/PBDunbtWqp2VMfXc+Ewta7EracV9qKfTuIqzD6u6B7y89omSRIKDVnQpWlHRkYGAQEBvPDCC+okPEEQQfix6vimVRUyslX838EYgiOzaF3fjLf7OmFsWDuGqH18fJgzZ456Dvfjjz/Gx8eHhQsXYm9vr+/mVarq+nree/hciTOptVFIEivmTqnQHvLzNisY2qcTlmVYVpWWlsaDBw+q1e9AqF7EcHQdZ2as5G+DXPn5WBxn7qTxn52RzBjogrVZzX9pjBgxgnHjxgEF+4x27NiRJUuWMGTIEObPn1+qZCyhYpQmk1obGfhg8dpiy5wmvzKY9IzMYoHy2XleTcdDQiN0JmClZWbzwdQxpXqujx49YsOGDSQnJ+Pi4lLrP/gJZSN6wo9V155DVZFlmZ0XE9l9KQl7CwNmDHLBzbZ2Vu05evQoCxYsoE2bNnz88cfl2iqtuqopr+fSDlE/j0JdFq5gyHhgr47sO3L+yTIjLcd1Pb5CoeCVF3uVailSWFgYGzduJDs7m6FDh9KuXbsKeHZCbVQ7xh2FcpMkiRc72PF6T0cS0/NYviOS4MhMYpJz2XYuge8DYth2LoGY5Jq/rrhnz57s2rULKysrhgwZou/m1GnaduoBypQomJ8vF9nJac/hc0V+1nZcp1ImYF2/fp3169eTn5/Pq6++KgLwc5w9e5YXXnihyh83Ojq6UndHKqmaP+YoVKhuXpbYmSv5v4MxrNwdBYWdB2pXzWkLCwsWL15cJ4NwTk5OtapNrCmTul2rpsz69Fut69srkyRJSAXflLoMZVJSElu3bsXMzIxXX31V4/I0QXiaCMJCMc3rmfFWHyf+u79gzWthR6E21pwubdZxTZeens4333xD69at6dGjB2ZmZvpuEqB5px7ta40l3UuHykEC6rs50bihW5kSvgrXfHt4eGBjY1MpbawI33zzjcbj48aNw87OjoSEBDZt2qTxnOnTpwNw9+5dDhw4UOx2W1tbXnnllVK1Jysri/fff5+QkBBMTEyYN28ezZs3Z8WKFRw/fhylUkmjRo2YN28ejo6OzJ49m/r166vbMnfuXFxcXJg5cyZ9+vRh+PDhnDt3jsjISPr27cu8efMA+O233/jhhx8wNzenVatWpWpjZRFBWNAoJCZb1JyuhTIyMrCxseHs2bNcuXIFf39/OnfujIFB9Xsr0EcPWVIo8O/YslTzv7m5uQQEBKg/1LRu3bqSWld7xcXFMWbMGLp06cKhQ4eYOXMmL7/8Mvfu3WP79u0YGRmxcuVKZs2axU8//fTc66WkpLBhwwbi4uJ44YUXGDduHCqVis8//5wdO3bg4uKic3vDqlT9/vKEakHUnK6dHB0dmTx5MkFBQQQEBHDo0CHOnz/PwIED8fb21nfziilJD/m5RTVKo5TzvxkZGWzcuJHw8HDMzMzo0aNHxbSjkhX2ILWxs7N77jlNmjSpsJ2PmjRpQpcuXQDo168f8+bNY926dSxbtkw9dfLmm2/SuXNnMjMzn3u9wlKujo6OODg4kJCQwPXr1+ncubN6i9BXX321WuwhLoKwoJGoOV17SZJE69atad68OefOneP48eNkZWXpu1klpqmHbG5uyve/7im2reDA3h3Zd/h8iY+XZv43ISGBDRs2kJCQQLdu3ejevXslP/Pa69n66rIsF9vSMT8/H1mWUalUxaYkcnOLdgqMjY3V32ubvqguoz/VoxVCtfO8mtP17atPYo9QNgYGBvj5+eHr66t+03r06BGHDh2ib9++ODrqfwMFbTT1kHt2bqNxPfD4EX1LdbwkwsPD2bhxI5mZmQwePJiOHTtWxtOsM+7cucOVK1fw8fFh9+7d2NjYMGDAADZu3Ii/vz9GRkb89NNP+Pr6YmFhga2tLTdv3gQgMTGRCxcuaNxl7Gndu3dnzZo1PHz4EA8PD7Zu3VoVT+25RBAWNNJVc1qhgB+PxIEE7RtZ6LmlQnmZmpqqv79z5w7BwcHcvn0bX19fevfujYVFzfgdawrMZTleEseOHSM3N5dx48bRrFmzMl1DeKJRo0b88MMPhIaGYmFhwapVq2jYsCH/+c9/GDlyJCqVCk9PT1auXAnAa6+9xvvvv8+AAQNwd3enc+fnTyE0adKEefPm8fbbb2NmZlZtPjiJYh2P1ZTiBlVNU83pjJx8/rsvmpRMFSM72THAp3Zu/lCTlef1HB4ezsGDBwkLC1Nvt9e1a9dqtaxJ37KyskhMTBRLkIRyq9Y94aCgIBYuXEh6ejrm5uYsW7aMRo0aFTsvMjKSAQMGFLntww8/xM/PryqbWys5WRtqzIKeM9KdVXuj2HYugbiUXF71d6hVmz/UZfXq1WPSpEkEBwdz6NAhjhw5QoMGDWjYsKG+m6Y3sixz8OBBrKys6NKlCyYmJiIACxWi2gbhnJwcZsyYweLFi/H39ycgIIAZM2awe/fuYr2uCxcu0Lt372qTcl4X2FkY8K/h7nx7KIYTt1KJT81j6gvOmBqJImy1gSRJeHt707RpU+7evasOwA8fPiQrK4smTZrUmdGPvLw8tm/fzvXr16lfvz6dOnUqlkgkCGVVbYPwtWvXUCqV+Pv7A9C3b18WL17M1atX8fHxKXLuxYsXiYmJYdy4cWRlZTF69GjGjx+v8bqxsbHExcUVOx4SElLxT6KWMzVSMGOgC5tOPuLYzVSW7YhgxkAXHETmdK2hVCqLlPY7dOgQYWFheHp68sILL9T63mBmZiabNm0iLCyMFi1aMGLECBGAhQql9yB89OhR3nnnnWLHp0+fXuwP3NnZmcjIyGJB2NjYmAEDBvDGG28QHR3NpEmTsLW11bhLzubNm1m9enXFPok6TKmQeNXfAUcrQ7aeTWDp9kjeHeCMmbGSU0/NJft5WeJcCyps1XWjRo3i8OHDBAYG8u2339KmTRv69OmDtbW1vptW4RITE9mwYQPx8fF06dKF/v3715nev1B19B6EC/d6fdauXbs4e/ZsseOaPoX++9//Vn/v6urK2LFjOXDggMYgPHbsWPVC7qeFhIQwa9as0jZfoGDosr+PDQ5Whvzvr1hW7Igkv5bWnK7rrKysGD58OF26dOHgwYNcvXqVW7du8d577xXJsq4NEhISSEpKYsCAAepCEoJQ0fQehLVxc3MjJiamyLHY2NhivWNZlvn2228ZPXq0eks6WZa1LsR2cnLCycmpchpdx7XzNEfu7ci3h2KB2l1zuq5zdnbmtddeIyQkhJiYGHUAjouLw87OrlihhZokLy8PAwMDGjduzMyZM2tlL1+oPqrt5EabNm3Izs7m1KlTABw5cgRJkoott5AkiRMnTvDzzz8DBZ9ef/vtN4YOLdv6P6F8wh7loG3ErrDmtFB7NG7cWL0KITc3l59//pn//ve/3Lhxo9I2WahMFy9eZM2aNaSlpQGIACxUumobhA0MDFizZg1fffUVQ4YMYfXq1axatUr9Cfvtt98mICAAgGXLlhEUFMTQoUN55ZVXGD9+PL169dJj6+suXTWlRc3p2k2SJDp27Eh6ejq///47//vf/3j48KG+m1UisiwTEBDArl27kGWZnJwcfTdJqCOq7XA0QIsWLdi8ebPG27777jv1925ubkV+FvRH1JyuuwwMDOjevTvt2rXj6NGjXLhwgf/973+0bduW4cOH67t5WqlUKnbs2MG1a9dwd3fnlVdewdzcXN/NEuqIatsTFmomPy9L7bsvyeBiU7uDcFBQEKNHj2bw4MGMHj2ae/fuaTwvMjKSqVOnMnjwYIYPH87GjRvVtz18+JAJEyYwZMgQhg0bxuXLl6uq+RXC3NycwYMHM336dLy9vbG3t1ffplKp9Niy4rKzs9mwYQPXrl3Dy8uLiRMnigAsVCkRhIUKVVhzWpJAIfHkK2ColFh/NI6jN1L03cxKUVhg5u9//zt79uxh2rRpzJgxQ+Pc6Icffoi3tzd79uxhw4YNbN68mWPHjgHw/vvvM3ToUHbv3s1nn33GzJkzycjIqOqnU24ODg6MHTuWbt26AZCWlsZXX33FiRMniu16oy8GBgbqYfQxY8ZgaFi7PyQK1U+1Ho4WaiY/L0uauJgUqzmtypf57/5ofj3xiKjEHEZ3ta9VpS5LU2CmsCQrgIWFBV27dmXv3r14eXlx69YtXnrpJQBatWpFw4YNOXLkiMYldzWh+Ezh2trk5GQMDAwICAjg/Pnz9OnThzZt2uhl7W1ycjJWVlYolUpeeeUVlEqlWAMs6IUIwkKl0FZzevYId9YejOHw9RRiknOZ0q/mlbqsiAIzbdq0Yfv27bz33nskJSVx/PhxnJ2diYqKwsHBoUiPrPAampS2+MzXX3/N3r171Y+/fv16jTvQvPXWWzg4OLBs2bISX/t53N3dmT59OhcuXODYsWNs376dM2fOMHDgQBo0aFBhj/M8ISEh/Pbbb/To0YNu3bpVm31lhbpJvPqEKmVhouTvg1359cQjTgansmx7BO8OdMHRquYMA1ZEgZlly5bx2Wef8eKLL1KvXj26dOlCeHg4sixr7JFpK5VYmuIzkZGRfPvtt8yZM4dr165pfX5QsFXctGnTGDduHL6+vjrPLQ0DAwO6dOmCj48PJ06c4OzZs6SkVN30RGBgIDt37sTIyIh69epV2eMKgjYiCAtVzkApMaGHA262hvxxJoHPtkUwrb8zzVxrdsWlkhaYgYKt8D799FP1Xr2ffPIJDRs2xM3NjUePHqkLRhReo1+/fhofszTFZ9atW4csywwaNEgdhM+dO8f8+fOJj4+nR48ezJ8/HxsbG7p37461tTU//vhjhQbhQqamprzwwgt07twZS8uCKmqxsbGcOXOGXr16YWVlVaGPJ8syR48e5ejRo1hbWzN+/HgcHR0r9DEEoSxq1jigUGtIkkS/Nja8O9CFvHyZL3dHcfJWwRD1tnMJfB8Qw7ZzCcQkV48EnpIoaYEZgNWrV7N+/XoAHjx4wJ49exg+fDjOzs54e3uzbds2AG7cuMGdO3cqZFvOgwcP0qpVK3VlOYAff/yR3r1707x5c3bv3q2epy7ssR45coSsrKxyP7Y2VlZW6p7/tWvXuHz5MqtWreKvv/4iOzu7wh5n//79HD16FFdXV9566y0RgIVqQ/SEBb1qXd+MD4e789990aw/9ggoyKauiTWnCwvMLFy4kE8//RRTU9NiBWbGjRtH3759mTVrFrNnz2bYsGHIssyCBQto3rw5AJ9//jkfffSROkh/8cUX5a7cFBMTQ0RERLFg/sEHHzB+/HhycnLo1asXBw8eRKVSqXdP2r9/P0FBQXTo0KFcj18Sffr0wcPDg0OHDnH8+HEuXbpEz549ad++fbl3LvLy8iI5OZmRI0diZGRUQS0WhPITQVjQO3c7I97q68TyHQXJR/k1uOZ0SQvMODs78+OPP2o8z8PDg59++qlC2xUbW1DP29bWtsjx+vXrA2BkZISHhweBgYEkJSVhb2+v7jEX3reySZJEs2bNaNKkCZcvX+bw4cPs2bMHV1fXMs3fpqamkp2djYODA56ennh6elZCqwWhfMRwtFAtXH2QgbbVSqLmdPkVrlXOy8srcvzmzZtAQd3n8PBwjI2NiwxXg/aksMqiUCho3749f/vb3xg1apQ6AD98+JCIiIgSXSM2Npbvv/+eX375pdqsSRYETURPWKgW4lNztVfaQtScLq/C5LBn1xSvWbOGmJgY7t27x6NHjxg9erR6jjYhIQEoSDjTByMjI1q1agUUfIjYs2cP0dHRtGrVij59+hTr1Re6f/8+mzdvJjc3lxdffFEU4BCqNRGEhWpBZ81pGewsxEu1PBwdHalfv36xMpoLFizgu+++IyYmhhEjRjB79mz1bcHBwZiYmODt7V3VzS1GkiRGjhzJoUOHCAoK4ubNm3Ts2BF3d3d1oL516xbnz58nNDQUQ0NDxo8fT6NGjfTcckHQTbyzCdWCn5cl+68kabxNBsITcsjOzcfYUMyglNWgQYP47rvvSExMZOnSpSxduhSAESNGFDs3Pz+fCxcu8MILL1SbRCYnJydeffVV7t+/z8GDBzlz5gwmJibqIJyenk54eDjOzs7qTHNBqO5EEBaqhcKa0+uPxal7xIVf69sbcf1hJst2RPJOf+caVdijOnnllVf48ccfCQgI4OWXX9Z57oULF0hISGDChAlV1LqS8/T05O233+b+/ftF6nK3a9eO9u3b67FlglB6IggL1Ya2mtOOVgbsDUziz/OJLNkWwdt9nWhRz0zfza1xXF1defPNN/nll1+eG4TXrVvHsGHDipXbrC4kSSo21CxqPws1kQjCQrWireb0YF9b6tsb8/1fsXy9N5qRHe3o72Mt3ni1KCx08exGDv3796d///5cv35d5/2nT58O8NzzhOqlUaNGmJrW7MpzdY0IwuUUFZ/O0cBI4pIzcbQ2pWdbN1ztxX6klaFVfTPmjHRnzYFotp5LIOxRNq/3dBTzxBqEh4cDFKsfLdRuW7du1VihTai+RBAuhyOBEXy36zoSEjIyEhI7T99nyrCW9PRx13fzaiVna0NmD3fnpyOxXLiXTnRSLu/0d0Ylw6mnhrH9vCxxrgHFPSqLv78/K1asoF69ehgbG6uPF27ssGLFCho3bqzHFlaO2v78QPdzFNngNY8IwmUUFZ/Od7uuI8sgP15YU/j1253X8fKwxcVOzFtWBhMjBVNfcGbv5ST+vJDIwj/Cyc2TkWpoucvKYGdnx4svvqj19saNG9fqHlNtf35QN55jXSDG8croaGAkEprnIyXgyOWSVfYRykaSJAa3s2V8dwdy8go+/uTLIMtPvq4/FkdsDdoAQhCEukcE4TKKS85U93yfJT++Xah8j1Lz0JabJcpdCoJQ3YkgXEaO1qZae8KyXLBnrlD5dJWzFOUuBUGo7kQQLqOebd209oQBjl+NYsWmS9yLSqnCVtU9heUuNZILbhcEQaiuRGJWGbnamzNlWEu+3fl0dnRB72t0zyaExaZy5kYMl+88on0zR0b1bExDFyv1/cXSporxvHKXCWl55KpkDMXIBFBQQ3rGjBm1dlP72v78oG48x7pEkp+u+1aHXb9+nZdeeqnU6+yiEzI4cjlCHUx7+bqrs6LDYlLZciyE87cK9mPt5O3ESz0acy8qpdjSJhlZLG0qo1PBqcXLXcrgZG1ATHIeDR2NmdLPSfSKBUGodkQQfqysQbgkQqNS2HIshIu343SeJ0nw+XR/sbSpDGKTc4uVu7S3NODPC4nsC0zC3FjBm32caOUh/m8FQag+RBB+rDKDcKF7kcms2nqVmETNmdMKCYZ29WRc36aV8vh11ZUH6fx0OI7MnHwGt7NhaDtbFAoxPC0Igv6VKzHr6tWrjB49mnHjxnH06FH18XfffbfcDauNGrlZ08jNWuuSGrG0qXL4NDBn7kvueDgYsftSEl/vjSY1U6XvZgmCIJQvMWvp0qUsWbIEAwMDPvnkE+Li4nj55ZdJSREZwdoULm3SlFkty2BqpNRDq2o/BytD/vWiG7+djufYzVQWbw1nSl9nLEyVotylIAh6U64grFQqadq0YOh07dq1zJo1i7S0NLGzjQ4927qx8/R9rbf/dTmCXFU+I/wbiWzpCmZooGB8d0caO5vwy/FHrPgzEpmCaYDaVu7y+PHj/Oc//yE7OxtXV1eWLVuGk5NTkXMSExOZM2cOYWFhqFQqZs2aRb9+/fTU4tIpyfM7f/48U6dOxcPDQ33s888/p0mTJlXd3DL78ssviYuL49NPPy12W03+/QlPlGs4Oi8vj/T0dAAMDQ354osvCAoK4tq1axXSuNqocGmTJIFCkh5/LUjKerGbJy097Th+NYpZa07xfzuCiE7IKHL/qPh0NgXcYdXWq2wKuENUfLqenknN1aWZJVP6OanHImpbucuEhAT++c9/snz5cvbt20fv3r2ZM2dOsfMWLlxIs2bN2LNnDz/88AMLFiwgKipKDy0unZI+vwsXLjB27Fh27Nih/ldTAnB4eDjvvvsuP/74o9ZzaurvTyiqXD3hefPmkZmZibl5QY9NoVCwYsUK9u7dWyGNq616+rjj5WGrdWnTzQcJbD12j2NXIzlxLQr/Nq6M9G/EzbBEsWtTBQmJyUYhFQTeZxWWu9S0r3FNcOLECby8vPDy8gJg3LhxLF++nLi4OPXa0ry8PA4fPqz+W61Xrx7+/v7s3LmTKVOm6K3tJVGS5wdw8eJFsrOzGTVqFJIkMXnyZAYOHKivZpfK5s2b8fPzo2nTpsTFFV9VUZN/f0JR5QrCmrKIJUli8ODB5blsneBiZ6Y1C7p5AzvmTrDj5oMEthwN4diVSI5fiVT33MSuTeUXn5qrtd5ZTS93GR0djaurq/pnIyMjbG1tiYqKUgepxMREsrKycHFxUZ/n4uJSI3pSJXl+AJaWlgwbNozhw4dz9+5dJkyYgLOzM76+vvpodql88MEHAKxatUrj7TX59ycUVSFlK1NSUjh37hzr16+viMsJjzVvYMe81zsyb0IHbK1MtJ4ndm0qPV3lLmUZrMxqbjE5WZY15mUoFIoi5wDFznv6nOqqJM8PYOXKlQwfPhyAJk2aMHjwYA4dOlQlbaxsNfn3JxRVqt+YLMvcv3+fvXv3snLlSqZNm0bv3r3p3Lkzb775Jr///ntltbNOa9HQDi8PG+1BA7G0qbT8vCx1VP6GS/fSCI3NqrL2VCQ3NzdiYmLUP+fk5JCYmIibm5v6mL29PcbGxsTGxqqPxcTEFOlhVlcleX7p6el888035OXlqY/JsoyBQc39cPW0mvz7E4oqcRD+8ssvadeuHYMHD+Yf//gHv/76K/Xq1WP69On8/vvvXLp0iZ07d1ZmW+s0R2tTrVnnsgymxrXjzaWqOFsb8noPxyKJcYVfOzcxJzVTxfI/Izl0LZmaVs+mW7du3Lhxg9u3bwPw+++/4+Pjg53dkzlupVJJ3759+fXXXwGIiIjg+PHjNSK7tiTPz8zMjC1btrB9+3YAHj58yL59+xg0aJA+mlzhavLvTyiqxBWzfHx8mDNnDj179uTGjRt8/PHH+Pj4sHDhQuzt7Su7nZWuKipmlUdUfDr/XHMSbb8tCfBv4yqWNpWSpnKXTtaGPHyUzbcBscQm5+LTwIyJPR0xN6k5a7hPnTrFihUryM7Oxt7enqVLl+Lu7s7w4cNZvHgxrVu3JiEhgfnz5xMaGkpeXh7vvvsuw4YN03fTS6Qkz+/WrVt88sknpKSkoFKpmDFjRo3LV1m1ahXR0dHqJUq15fcnPFHiILxgwQIWLlyo/jklJYUlS5Zw5MgR5s+fX+Ne3M+q7kEY4OiVCI27Ng3r2pCQyBSuhyYgSdCtlSsjuxcNxmLXptLLyslnw4lHnLubhp2FAZP7OtHYWfvcvCAIQmmVu3b00aNHWbBgAW3atOHjjz8uMiRUk9SEIAy6d226FZbIlmMhXL9fEIz9Wroyonsj7oQniV2bykiWZU4Gp7LpZDyqfJkRHe14wceauJQ8UWlLEIRyq5ANHNLS0li6dCkBAQGcPn26ItpV5WpKEC6J4IeJbD12j2v34nWeJ3ZtKrmIhBy+OxRDVFIu7raGRCbmIj1VaUuGWlFpSxCEqlUh+ewWFhYsXryYL774oiIuJ5STl4ctc8a35+NJnXC0EUubKoK7nRFzRrrj29CMiMSCNca1rdKWIAhVr0IXlXXt2rUiLyeUUzMPG5q424hdmyqIsaECZxsjrUvFCittCYIglJRY2V3LFe7apIksF2xqIJRcfGou2qJwTa+0JQhC1RPvwLVcz7ZuGrdNLHTsSiQrNl3ibkRS1TWqBtNVaQu54HZBEISSEhUearnCXZs0LW0a1aMxD2PTOHczhst3HtG6kT0v9WiEl4et+v5iaVNRfl6W7L+SpPE2GYhOyiErJx8TI/H5VhCE56uQ7OjaoDZlR2uia2lTeFwa20/c4/T1aGQZWjS05aXujYlNzOC73TfE0qZnnApOZf2xOPWHGYmCoX03O0MiEnJxtDLgrT5OeDqJNcWCIOgmgvBjtT0Il0Tko3R2nLzHyWvR5Ot4WYilTZorbTlaGXD4egpbziaQny8zrIMtA31sUCi0DmALFejgwYOsWbOG0NBQPD09+e9//8v8+fMZPnw4Q4YM0XfzBEEjEYQfE0H4iZiEDFb+HkhYbJrG2xUSDO3qqXUrxrouIiGH7wNiiEzMpamrCW/2dsLOQsz8VKYdO3awdOlS1q5di7e3N5MnT8bQ0JC0tDQ2b96s7+YJglZi4kooxtnODHdHC7G0qYwK1xT3amnFnagsPvkjnIv3NH+gEcpPpVKxYsUKZs2aRZs2bTAyMqJt27acOHGC2bNn67t5gqCT+HguaFS4tElTZrUsQ3aOClV+Pkqxf6lGRgYKXunmQCsPM9YdieXbQ7H4NcugdytrLt5LF+UuK1BgYCAJCQlFdkjKycmhf//++Pr66rFlgvB84h1U0Oh5S5su3Ynjn2tOcSQwgjxVfhW2rGZpXd+M+S/Xo5WHKadup/Hp1gj2X0niwr10DlxJYsFvDzklCnyUS3R0NNbW1piamgKQlJTEtm3baNpUTJcI1Z8IwoJGhUubCvbZlYrstztpoDcv9WhEakYO3+68zvv/PcHBC2Hk5KnU94+KT2dTwB1Wbb3KpoA7RMWn6/HZ6JeVmQGjuz7Z7lMW5S4rlKurKwkJCZw/f56MjAxmz56NoaEhkZGR+m6aIDyXSMx6TCRmaaZraVNGdh6HLjxk95lQUjNysbEwYkjXhhgZKPhp3y2xtOkp284lcOBKEvka/toUEvT3sWFkp5q5A5m+ybLMihUr2LFjB5mZmYwePZpBgwYxbdo0Ro4cyYcffqjvJgqCViIIPyaCcNll5eRx+HIEu06HkpiarfW8ury06fuAGC7cS0fbX1uHRua83c+5ahslCILeieFoodxMjAwY1LkBK2f406aRvdbz6vKuTTrLXQJ3o7OIScqpsvYIglA9iCAsVBgjAyXmpoZiaZMGfl6WOtLcIClDxSdbIjh0NYl8TWPWglDNxCTnsu1cAt8HxLDtXAIxNSyvISUlhb///e/4+vri7+/PTz/9pPG8nJwc5s6dS58+ffD19WXIkCHs3LmzwtohligJFep5S5si4tKIik+vc/Wnna0Neb2HY/Fyl8DrPRxxtDJg3dE4fj+TwOXQDCb2dMRJLF0SqqmTwan8/Mxref+VJF7v4Yifl6WeW1cyixYtIicnh+PHjxMREcGkSZPw9PSkZ8+eRc7Ly8vDycmJdevWUa9ePS5evMjUqVOpV69ehSyBE3PCj4k54YoRFZ/OP9ec1Dr3CQVzw52buzC8mycNXGrGH2xF0VTusjDYZufms/18An8FpWColBjZyY7eraxQaBtaEAQ9iEnOZcFvDzX+jUsSLBrjUaEfIPfv38+aNWvw9vbm9OnT5OTk8P777zN69OgyXzMjI4NOnTqxdetWmjVrBsDKlSu5f/8+X3/99XPv//bbb9O1a1fefPPNMrehkOgJCxVK165NU4a1xNnWjB0n7nHmRjRnbkTj29SBF7t5qnduqu27NjlZG2rNgjY2VDDWzwHfhuasOxrHb6fjuRyazsSejjhaGRKTnMuppwK4KPQhVKRdFxO5UILKbimZKq0fsmUZlu6IwMpUqfX+HRpZMLS9rdbbnxUcHMzdu3eZOHEiS5YsYc+ePcyZM4eBAwdiafnkQ/zUqVO5ePGixmu4uroWGUIODQ1FlmV1AAbw9vbmwIEDz21PRkYGQUFBvP766yV+DrqIICxUuJ4+7nh52Gpd2uT9anvuR6Xw58n76m0Uvevb4ulqxb5zD4osbdp5+n6dW9rUzM2Uj16ux7azCRy5kcKiP8LxaWDGhXvpNXr4T6gdVCrdg6fPu720goODGTNmDCNHjgRg6NChzJ07l/DwcJo3b64+b+3atSW+ZkZGBhYWFkWOWVlZkZ6uu55Bfn4+s2fPpnXr1vj7+5fiWWgngrBQKVzszHRu8ODpasXfX/Yh8lE6O0/d5/jVSG6FJQKo55MLv3678zpeHrZ1ammTiaGCV/wdaNfInP8djuV8SMGbQ+HbW+HX9cfiaOJiIuaPhXIb2t62RD3U561579XSukLXvAcHBzNhwgT1z6mpqWRlZWFnV/bHMDMzKxZwU1NTMTfXPuomyzILFiwgNjaWH374AamCpolEdrSgV24O5kx9sRV92tXTuoSnLi9t8nIzpUMjC623SxQkyQhCVdGV6S8D3SpwZCY9PZ3w8HAcHR3Vx/bt20eLFi1wdi66rn7y5Mn4+vpq/PfsVpYNGzYE4M6dO+pjN2/e1FrqVJZlFi5cyM2bN/n+++91BuvSqtY94aCgIBYuXEh6ejrm5uYsW7aMRo0aFTvvpZdeQqV6UjLxwYMH9O/fn+XLl1dlc4VySM/Ke5Iu/AxZhuiEulv2MjkjD0lC4zycDMSn1qylIULN9rxM/4oclbl9+zYKhYKdO3cyc+ZMLl26xJdffsmKFSuKnfv999+X+LpmZmYMGDCAlStXsnz5ciIjI/njjz9YsmSJxvMXLVrElStX+Omnn4oNY5dXtQ3COTk5zJgxg8WLF+Pv709AQAAzZsxg9+7dxYYBtm7dqv7+xIkTLFy4UGxhVsPoXNpEwYYRvx2+w4BO9bE2N676BupRYaEPTb0PWQZLHUkwglAZ/LwsaeJiojXTv6IEBwczaNAgwsLC6NChA+7u7ixevBg/P79yX3vBggXMmzeP7t27Y25uzttvv11kedLkyZPp0KEDw4YN49dff8XIyIhevXqpb586dSrTpk0rdzuq7RKlixcv8q9//YuAgAD1sd69e/Pll1/i4+Oj8T5paWkMGTKE//znP3Ts2FHjObGxscTFxRU7HhISwqxZs8QSJT3RtbRJAuytTXiUnIWhgYJebd0Z0qUBTrZmRe5fW7OqdS0JATAzkni5qwN+zSwqbJ5KEKqDhQsX4urqypQpU/TdlEqj957w0aNHeeedd4odnz59Oq6urkWOOTs7ExkZqTUI//DDD7Rv315rAAbYvHkzq1evLl+jhQr3vKVN3Vu7cSE4lj9P3ufghYcEXAynS0tnXvTz5F5UCt/tul5rs6p1Df/5e1kS+CCD9UfjOBWcyqv+DrjbGem3wYJQQYKDgyssC7m60nsQ7tmzJzdu3Ch2fNeuXZw9e7bYcYWWTeRzcnLYsGEDGzdu1Pl4Y8eOpU+fPsWOF/aEBf153tKmTs2d6ejtxPXQBHaeCuVUUDSngqLV96/NWdW6hv9GZqvYcT6RYzdSWLwlnH5trBnazhZjQ5F3KdRst2/fpkmTJvpuRqXSexDWxs3NjZiYmCLHYmNji/WOC504cQJPT08aN26s87pOTk44OTlVWDuFivW8pU2SJNHK055Wnvbci0rhm21XiYzP0HwuBVnVuq5Xk2gr9GFurORV/4Lh6A0nHnHgSjLn76YxrpsDbRuaiyIfQo114cIFfTeh0lXbINymTRuys7M5deoUfn5+HDlyBEmStM7Xnj17li5dulRxKwV9auRqRQMXK6ISMrRmDtelDSMaOpkwZ4Q7R26ksON8AmsOxFDPzoiIhJyC7GpEkQ9BqG6q7XiVgYEBa9as4auvvmLIkCGsXr2aVatWoVQWZIK+/fbbRZK2QkNDqVevXpW3MzI2iXU7TrPixwOs23GayNikKm9DXVaYVa2JLENCShapGXVni0CFQqJPK2sWjfWgtYcp4Qk5yEC+XPD/Ufh1/bE4YmvYrjeCUBtV2+zoqlaWDRwOnb7J178eUa/hLPz6t/G96dfFu5JbLEDJNowwNlTQs607gzo3wNm25s8Pl9S2cwnsv5Kk8f9GIUF/H5sKrWwkCELpVdvh6OouMjaJr389gizL6je5wq9fbzhMi8auuDla66+BdYSurOrJQ1pgbKhk1+lQDpx/yMELD+nk7cxQv4Y0div43dTmpU26injIsijyIQjVgQjCZXTw9E2tVYwkCQ6eusHE4V2rvmF10POyqru2dOFGaAK7zjzg7M0Yzt6Mwbu+LQ2cLTlwIazWLm3SWeQDiEjIISEtDzsL8TYgCPoi/vrKKDYhVfsQqFxwu1B1dGVVS5JES097Wnra8zA2ld2nH3DiWu3fMMLPy5L9V5K03h6ZmMv8zQ95oY01A9raYCKWNAlClRN/dWXkZGeJtuJE+bKMuWndKq1YU3g4WTJteCv6tveo9RtGFBb5kKSCOeCnv07s6cjMQS44Whmw53ISH216yIlbKeQ/tTVOTHIu284l8H1ADNvOJRAjErkEocKJnnAZvdC1OVsOBWq9/cDpGxgaKhkzoD3WFqZV1zChRNIyc7VuGJEvQ1hs7RjJeF6N3+buppy8lcqfFxL4+dgjDgel8HIXOxLSVfz8TIUusbRJECqeyI5+rEzZ0Wdu8fWGwwU94sfvVLIMI/v6cDU4grsP4zAzMWLUC7682LsNJkaiQEJ1sSngDrtOh5Kv4+Xv09ieQZ0b0LqRfa2vyZyZk8++wCQOXUsmT8em7JIEi8Z4iP2LBaGCiCD8WFmCMEBkXDIHT90gNiEVJztLXvBrgZujNfn5Micv32X9zrNEP0rBztqcV4d0pF9nb5RKBZGxSRw8ffPJ/bo2x83JpvKeoFCEzg0jJGjdyJ5r9+KRZajnaM6gzg3o1toVIwOl+v61Mas6PjWXr/dGE52keehZLG0SCtX0v4GUlBQ++ugjjh07hrm5OZMnT2bSpEnPvV9CQgKDBg2iQYMG/Pbbb+VuhwjCj5U1CD9Pbp6KfSeus2nfBVLSsvBwscXHqx67jwWJ9cV6dvRKhNYNI3r6uBObmMH+8w85cjmczBwVVmaG9GvvgbmpIb8cDC6SVS0j15qs6u8DYrgQkq4xq1qSoEMjcyb3ddZwq1BXHAmMKLZpSk37G/jnP/9Jeno6K1asICIigkmTJrF06dIi2xlqMmfOHMLCwsjNza2QICzmhCuZoYGSYb3a0LezN1sDAtkacJmHR68BiPXFeva8pU1OtmZM6O/FqJ6NORoYwb5zD9h6/J76/rU1q9re0lDr8jtksDITbxt1WVR8Ot/tuo4sV83fwP79+1mzZg3e3t6cPn2anJwc3n//fUaPHl3ma2ZkZLBv3z62bt2KhYUFXl5ejBkzhi1btugMwufOnSM0NJSXX36ZzZs3l/nxnyb+mqqImakRrw3tRHpmNruPXtPayxDri6vW8zaMADAzNmBQ5wb07+jBV39c5UJwrMbzasuGEbqWNsnAqeAUrEyV9G5pJXZqqkW2Hgvh9PXo556XnJ6jdXmmLMP8/53F2lz7dppdW7rwUg/dG+08LTg4mLt37zJx4kSWLFnCnj17mDNnDgMHDsTS8kmS4NSpU7l48aLGa7i6urJz5071z6GhociyTLNmzdTHvL29OXDggNZ25OTk8Mknn7BixQqNO/+VlQjCVSwlLRNt3Yz8fJmLN8IY6N8SZ3srPbRO0EWpUGBooNDaS5RliEpIr/qGVTBd+xd3bWrBrcgstp1L4NC1ZAa1taFHc0sMDUQwritUqvxy3V5awcHBjBkzhpEjRwIwdOhQ5s6dS3h4OM2bN1eft3bt2hJfMyMjAwsLiyLHrKysSE/X/vf77bff0rVrV7y9vUUQrskK1xdr+yR5PyKeyQt+oVkDJ7r5Nsa/XROc7J582hMJXfpVuGGErGEsQwYuBsfy/e4b9O/gQX3nmruUR9fSplyVzMlbKey5nMRvp+M5cDWJIb62+HlZYqCUxNaJNdRLPRqXqIeqa2WBQoIXOtSv0NGg4OBgJkyYoP45NTWVrKws7OzKnhxoZmZWLOCmpqZibq45sezBgwds27aNHTt2lPkxtRFBuIrpWl8sSRLTx/bgekgUZ6/d58ftp/lx+2m8Gjrj79uYfFnmpx1niiR0bTkUKBK6qlDPtm7sPH1f420S4O5gwV+XwvnrUjgtGtjSv2N92ns5olQU9BRrUkaptv2LDZUSvVpa4+dlybEbKewLTGLDiUfsv5KEl5spp26nivXFtZiuvwEZ6OVbcYlZ6enphIeH4+joqD62b98+WrRogbNz0eTAyZMnax2OdnNzY/fu3eqfGzZsCMCdO3do2rTgA8PNmzfV3z/r4sWLPHr0iAEDBgCQlZVFdnY23bp1Y//+/cV61aUhsqMfq6zsaE20rS9+Opjm5OZx6eZDTl6+y9mroWRma69WJEkS/zf/VZHQVUV0ZVX3aOPG3Yhk9p8P4+yNGFT5Mg7WJvRr74GRoYKfD9S+rOqs3HwOByWzLzCJrFzNbydifXHt8ryVBRXl8uXLjB8/nilTpjBz5kwuXbrEe++9x4oVK/Dz8yvXtT/44AMyMzNZvnw5kZGRvPHGGyxZskRjYlZWVhapqU8K+OzZs4cdO3awdu1aHBwcylVHQAThx6oyCIP29cWaZOfk8fm6g5y+ovnTp0KSGNmvLZNEQleViU7I0JpVXSgpLZuAi+EEXHpIUpr2PY0lCT6f7l/js6p/O/WIv4JSNCYdivXFtU9J/gbKa9OmTZw/fx5Zljl8+DDu7u588MEH9O7du9zXTklJYd68eRw/flzjOuHJkyfToUMHpk2bVuy+W7duZdOmTWKdcEWq6iBcWit+PMDxSyFo+3UZGxkwwK8FXdp40qKxK0pl0UQZMZesP3mqfL784wqXbsdpvF0hwdCunjU+q/r7gBgu3EvXmu/g29CMaf1dqrZRQo22cOFCXF1dmTJlir6bUmnEnHANoSuhS6Igc/fPI1f588hVLM2M6diqIV18PPH19uDEpbt8/esRMZesJwZKBcaGSq2/v9pSq1rX1okAV8My+O10PP1aW4vtE4USCQ4Oxt/fX9/NqFTiL6GG0LlhhCTxxb9eJicnjzNX73Pm6n3+OhfMX+eCMTRQkJtXsGRAFAfRH11Z1QCBdx/x6S8X6NuuHh28nDBQ1rwlP7rWF0uAg6UhAdeSORyUTIfGFvT3scbD/sluYyKrWnjW7du3adKkib6bUanEcPRj1X04GkqW0FUoNiGVM1fvsz0gkLjENI3XUygkXurbVhQHqQI6a1UDPk3suXYvAVW+jLW5ET3butPH1x0nW7Mi16jumdWnglM1ri9+vYcjXZtZcD08kwNXkgiOzAKgRT1T+rexJiEtj5+PP9J4P5FVLdRmIgg/VhOCMJQuoQueP5fs5mjN1NHdadXUDSPD4gMjYi654jwvozQ5LZujVyIJuBROXFImEtCmsT1923uQmpHD97tv1IjM6tjkXK1bJxZ6EJfNwatJOueQQWRVC7WfCMKP1ZQgXFrrdpxma0Bgkc3aNTExMqCttwcdWzWgQ8sG2Fmbc+j0zWJzyWKjifIpSUZpvixz7V48f10M5+LtOJ3bLdb0zOpHKbmsPRRD2CPN2eMiqxq2/HWVn7f+RVZ6GibmFkx4qQ+j+rTRd7OECiLmhGu55xUHWfDOYB5EJnDuWijnroVy5mrBMqgGbnY8iEwAxFxyRSpJrWqFJOHT2AGfxg4kpmbx9ZarBD9M0nhuTa9X7WBliLO1IQ/jNdcjluWCueK66u2Pf+b29cvqnzPTEvh6zU/sP+bLtx9P0HFPoaaoedkfQqm4Odnwt/G9kSQJhUJC8firJEn8bXxv2rdowEv9fFn6j5H8svQNPpjYjx7tmxIZm6z1moUbTQiVz9bSBDsrE7TVAsiX4X50itbphpqgMKtaExm4EprOd4diuBaWgUrDiE5Mci7bziXwfUAM284l1JqgveWvq0UC8NOCr19m61/XqrhFQmUQPeE6oF8Xb1o0dn3uXLKluQm9OjajV8dmLP/fAU5cvqt5SU2+zIUbYfRo35SG7vYaq8WIueSK87zM6mv34vnX/52ih48b/q3dsLU0LnJ7dU/o0pVVDVDP3ogL99K5cC8dSxMFHRpb0LmpBQ0djTl1O42fn0kEqy1lMn/e+pfO29dvDeClPq2rqDVCZRFBuI5wc7QuVRa0s70lkiRp7WGFRsTzt6W/YWdlRrsW9Wnfoj5tvT2wMDPWOJcs1iWXnc561RJ0a+VK4N1HbAy4w6a/7uDT2IEePm60b+bEyaCoYpuv7zx9v1oldOnatakwmD5KzeXc3TTO3knj8PUUDl9Pwd7CgPi0PODJ2uTCr+uPxdHExaRaJXSFR8Wx56+zRMUl4upoy+A+nann6ljsPJUqnxt3HpCSoHnLzEJZ6ZpXPQg1i0jMeqy2JmaVVWRsEtM+2agxCEuSxN/H9yY0Mp6LN8J4GJ0IFCx58nR3IOSh5spQosZ12T0vszo3L5/Ld+I4eiWCK3fjyZdlTI2UZOaoNF6vOiZ0lSSrWpZlwh7lcOZOKidupZKTp/ntqyQJXVW5LnnP4XMsX7P58XNAPb3w4TtjGdS7E4nJqZwNvMXuw5e5HhyCKu/5Q+q2zvXZvvq9SmmvUHVEEH5MBOHiSrouOSY+hUs3wrh4I4wL1x9onLcDsS65vEpaqzcpLZuT16L48+R9UjM1v5nXhlKZ3x2K4eK9dK0VuhytDBjsa0tDJ2NcbAxRPDVtcjI4tdgwdknXJZc2eIdHxfHa35dqXybo4khk9JMPrgoDYwyMzTGztCEp+p7W6/79nTfEcHQtIIajBa1KOpfsbG/FoO6tGNS9Fct+2M/JwBCtc8lnrt6nTbN6tGjsirGRWJdcGiXJrAawsTBmSNeG3ItM5szNGK1Zxw/janapTAcrQ517c8el5LHuaEFwMzGUaOhkgqejMbbmSjaejEem9MPYmoL38+agf9t9QmfiXGR0HAbGFhgYm2Fp64x/p1a8OaQNVubGTPn4Z4I1JGd5tfQVAbiWEEFY0Km0c8kuDlY655LDY5KY/9+dGBkqadHYFV9vD3y9PWjgZs9fZ2+JueQK5GhjpjWhSwYu33nEv787TZcWLnRu4YyzbdFedU1O6JIkeH+oK6mZKu7HZnM/NpuQ6CxuRWQ+97p7LyfyYgc7zE0UGBk8WUASk5zLz8fiyE5LIinyBrmZqRiaWmLj1oL1xygWvHNz87hy8z4Bp67ofDwDE0s6de3Bm8N8aepRdPj8248nsPWvtqzfGqBeJ/z6S32rRQCu6R+YU1JS+Oijjzh27JjGXZQ02b9/P6tWrSI8PBxbW1vmzJlD//79y9UOMRz9WFmHo0uabFFXPG8uee6UgYRHJxF46yHXQ6LIzSuYs7Q0NyY1PVvjNcVcctk8r1RmB28nrt9PICO7ILmpkasVnVs407m5MzceJBZL6KqOFbp0lcl8tmeqypeJTMjh52NxPNBSHORZRgYS5sYKLEyUZGTnExJ8legbhx/fWviI4NKyD028WqPISyXq4X3iokJJiXtIvur5c7vWDvX4c837JX3K1UJtKOTzz3/+k/T0dFasWEFERASTJk1i6dKlGvcTBjh9+jSzZ8/miy++wNfXl8TERDIyMvDw8ChXO0QQfqwsQXjP4XOsWLOZZ1+J/3qcbFFXlXQuOTsnj+shkQTeCuevs7dITsvSeD0xl1x2JUnoCrofz5kb0VwIjiPzcUDWpqYmdD1t27kEDlxJQlPqggQ0djGhvoMR6Vn5pGWrSMvKJz1LRVTMI0JO/qL1ugbGFuRlp6mvZGhqiYGxOcj5ZCbHaL2fS6O2bF72egmfrf4974N2RX9g3r9/P2vWrMHb25vTp0+Tk5PD+++/z+jRo8t8zYyMDDp16sTWrVtp1qwZACtXruT+/ft8/fXXGu/zyiuvMHLkSMaMGVPmx9VEDEeXUXhUHCvWbC4oKfhMSanlazbT2tuzzvaISzqXbGxkQLvm9WnXvD7xSWlaa1zn58vsO3kDVb5Mm2butGjsipmJUZFzavrQWGXp6eOOl4et1oQuQwMFvk0d8W3qSG5ePlfvPWLzX3cIj0vXeL3qWKHLydqwVGUtda5LlmBiT0eNQXzmklM86WsXp8rNwtTGFaWhKQrFk2FsY2NjZFkmK6X4kiMTKyfq13Mtcdsr08a95zl+8e5zz0tKzdQ63STLMv/8zxZsLE213r97+ya8MqhjidsVHBzM3bt3mThxIkuWLGHPnj3MmTOHgQMHYmn5ZLRj6tSpXLx4UeM1XF1d2blzp/rn0NBQZFlWB2AAb29vDhw4oPH+KpWKa9eu0bt3bwYMGEB6ejrdu3fn3//+d5E2lIUIwmW056+zaM0KkWBXwFmmvTZU6/1r+zB2aeeSde2XDJCTm8e2gEC2BQSiVCho1tCJ1k3dadPMnci4ZNZsPibmkrUoaUKXoYGC9s2cOBUUTcQjzRsr5Mtw7GoEtpbGtGlsj4udWbFiLdV9LrlwXfLaXTdICrtKTlYaRiYW2NRvw9ShLYoF4PjEFM5fDSbk9g2075YMBkamGBoXPE83Jxt6d2xGjw5NeZSSy9yVv2FgaEJOZgr5qjwUSgOMTK1QGBjx+hDfyny6FU6l0rzsraS3l1ZwcDBjxoxh5MiRAAwdOpS5c+cSHh5O8+bN1eetXbu2xNfMyMjAwsKiyDErKyvS0zV/+Hz06BG5ubns2bOHdevWYWZmxgcffMCSJUv47LPPyvCsnhBBuIyi4hK1Roz8fJnfdx3hXlgkbZo3orV3I7wbe2BsVPDHrWkYe+OOw3V6GPt5Na6/nj2WnLw8rt2O4OrtCK7dieTmvWh+2//kk6+ocV0xnlehKzkth3X7bxWca2NKm8b2tGnkQEtPO87djKn2xUEAzp45ReiZ/RT2bDOQSIq8yTnHAbRv2Jurt+5z/kowZy7f4kF4dImuaWBoxIg+PvTs2IzG9RzUH07cHGFIPz92HzqFiaVDkfsM6edH60YOmi5X5V4Z1LFEPVRdm8IoFBKDu7eq0Kmj4OBgJkx4Uic7NTWVrKws7OzKvqmHmZlZsYCbmpqKubnmD4umpgU9+/Hjx+Pi4gLAtGnTePfdd8vchkIiCJeRq6Ot1p6wJIGNtSUXr97h7OWCNytDAyVejT1oWM+F3X+dLRjOEcPYaoU1rrXNJbs72wDg6e7Ai719UKnyCQmP48ftpwm6E6nlqjIbdp1l5qu9MTHWPEcohrGLe16FrsWTuxCfnMXVkEdcvRdPwMVwAi6Go5BQz7MWBvDCr9/uvI6Xh63OueSy9qBLe79LN8LYvms/hS19+uu2nfvZsfsg+fn5ACgMjDC1ccPYwh45N5vkmNtar9usmRdvvdRN423TRrSlW5t6bDpwjfikVOxtLBnXv3W1CcCloesDsyzDC34tKuyx0tPTCQ8Px9HxyXvivn37aNGiBc7OzkXOnTx5stbhaDc3N3bv3q3+uWHDhgDcuXOHpk0LRolu3ryp/v5ZVlZWuLq6aizRW14iCJfR4D6d2bjjsMbbJCS++ng6jnbW3LoXztWb97h26x5Bt0IJCg7VflFJYs9fZ5kyvm4OY5d0LhlAqVTQrIEzdlZmWpdEyTIcu3iXk5fv0aS+I62auNGyiSvNG7mK8po6uNqbM2VYS60JXZ4uVni6WNHBywkoKCJyNeQRe88+ICZR+xKgP47cZeJAbyzNjIrddiQwokw96LLcb/32ozqfv6Q0wsqpAW7u7ng4WZOTnU5YRAzxSfmYWDlpndt1sNX9gaF1IwdaT+ut85ya4HkfmCty5On27dsoFAp27tzJzJkzuXTpEl9++SUrVqwodu73339f4uuamZkxYMAAVq5cyfLly4mMjOSPP/5gyZIlWu/z8ssvs2HDBnr27ImpqSnfffcdffr0KdPzeprIjn6sLNnRex+XolMHgedkR6tU+cxe+j3nr9zSOvdpY2XOwJ4dad60Ac2b1MfJwUb96UtkYxenc2hMkmjZxBVTYyOuh0SRnlmwBEqSwN3JhvCYJI3XFEuiCpS0QlehVVuvcuZGtNbXdiE7S2M8nC2p72RBfSdLzEyU/GdzoOb0Ch3Z2IVLsHKz0klPCEeVk4nSyBRzu3oYmpizYFInZFkmJiGTh9GJXL8dSuiDMCIf3kfO154FbuPgQrdOvgTeDCM1o+A142RniZWFCSEP41Dl5hSb21UaGjOqX93K4I+MSy7RB+by2LRpE+fPn0eWZQ4fPoy7uzsffPABvXuX/8NMSkoK8+bN4/jx4xrXCU+ePJkOHTowbdo0APLy8li6dCl//vknSqWSXr16MXfu3GJzy6UlgvBjVbVO+NsNu9j45xH1cNezDA2U6rWzAHY2ljRvUh83Fwf+2H1MY49PIUn8/NXsWtMjLo2SLpfIz5d5EBXP9btRBN2N5ML1B2TnaH4jVkgSI/u2ZdII7W+oYhi7uE0Bd9h1OrRgxcAzJAlaNLDD3tqEsJhUIuLSyVVp/hsocj+gsbs1rTztAAmFVHBQIUlcuxfP5StBJIYH8SRrueCrjXtLDEzMyU6NJzstnpyMZAqHnCWFEjlfe/KQkZkNJpYOeLrb06WNJ118GuHpbk9UXHKVLs0RYOHChbi6ujJlyhR9N6XSiCD8WFXVjg6PimPC35dqfKNSSBLrv/wQlSqfm3fDuHEnjFt3wwgJi0Sl4w1LoZB45cXedXYYu6Trkp+ma6tGAKVConVTd7w8XfD2dMbb0xkLM5OCx6sFhQoqg87iIM/0aFX5+UTHZxAWm8bWYyFEPNKclapLXnY6McHHn3ueibERLZo1oouvF13aefPzznMcDNA8lQTg07YDc6cNxdneqthtZXmtCWX36quv8tZbb9G3b199N6XSiDnhKlbP1ZF/vTNW6zC2h1vBPFtDDxf1EHN2dg5zlv3ApaA7Wmsy/7HnOKHhMTRp6E5Tz4J/zg62SJJU67OxSzOXXEjXVo0SYGtlxo17UQQGh6uPe7jY4uFiy+nAewV1h0U2dhHPm0t+ekhZqVDg7miBu6MFD6JTiYoPJScrrdiwspGJOb3aujOka8OCpC9ZRpYhX85n4dcb0V4CA1xcXPhoxhg8PVy4/SCOa3ciWL3pBMH3Y7TO7ZpaOeHh5qAxAEPZXmtC2d2+fZsmTZrouxmVSvSEH6vqXZQqchhbAqytzElNzyzSY7Y0N6WeqyM374ZpvGZJhrFraw+6JMPYTnYW3A+P59b9aG7dj+FWaDSx8do3PVBIEiP7tWWSjnnBujCMXdq55Kj4dN75dDOJD4sPK9t6tOL/5o3Dxc6MuPgkLgXd5VLQHS4H3SXmUaLOdri5ueHRoBG3Q2NRPf67MTc1wsrclOj4ZI1zuwZGxqI6m1ClRBB+rLpvZfi8Yeyfv5qNo70ND8JjuHM/nDv3I7gbGsnNu2HkaVk8L0kS/h1b8varQ3B3ccBAqSxye21PBCvL0OLib/dw7lqo1mFsA6WCVk3daOzhSBMPR5rUd8TZvmBTCzGMrdnztvrr082XO/fCeRj1ZLu/+m5OpGXmkJCYpPW6RmY22DkWZMS3bupO66ZueNZzIOZRipjbFaoNEYQfq+5BGEqfjQ2wcOV6Dp++onMrNQADpRIPN0c8PVxo6OGClYUZX/1vW5kTwWpKD7q0GZ66srElwNbajLSMHHJynyR9WZgZU8/Zhlv3NQ+elvSNv7b2or/dsItfdxzW+Rp1drClXasmNPash6WlNY+SM9lzJJCIB8Fa79OpU2eWvj8apVJR7DYxtytUF2JOuAYZ1LsTrb09SxXcXJ3stM59KhQSndt609DDhdCHMdx/GM1fpwJL1Jat+07wtzdGarytJs1Bl7a8pq5CBUgSn703Emc7S8JjErkbFsfdh3GEPIzjdqiu2UuZ9X+eYdLwrjjZWaJQFC8IUJPWNF8MCuGnPwJ4lJCEg50Nk17uS/tWjYudl5KaztVb9/nrlO4PiS7OTjRu6s2VB/GcvHlNfVyS0D63a+1Mq2b1NQZgEHO7QvUhesKP1YSecFmUZBj76SCekZlNWGQsX/6whZt3NM8lF7KyMKOeqwP1XJ0ef3XEyNCAj/7zU63uQVdGNnYhI0MD6jnbqJPAPFxsMTYyZOGa3WUePi1rD7os9/vPd9vZeeBYseMvDujJ5LH9uHIjhMDrIVy5EUJIWNRzR2igYFjZxt4Fz3oONKrnQGMPRxrVc8DAQMnMJZtR5WZrXLcrhpWFmkAE4cdqaxCGsg1j60wEkyQ8PVywMDclPDKOhGTtyUpPU0gSL/bvyt/eeEljD6U8c9BVHbwrchhbIUn4eNfDw9mWsOgEwqMTeZRUsiU7kiQxwK85b43qholR8dKcZZ2HLsv9LgaF8P7C/5ao3QaGhhgamSErjJCUBmQmRWk9t0vnLix5/2WUCjGsLNQ+Igg/VpuDMJQ+SJWmB52ekUV4VBzh0Y/4ectB7j/UXfDeQKnE2cEGFyd7XJxscXWyx8hQyZqfd5WpB10TEshKuwdrRmYO4TGJhEUnsuXgJa3VvZ5mbmqEnbU59jbm2FmbY2SgZP/JGxq3YZAkiS9mjaKesy2GBkoUCkldma2wrXm52eQ+1cM0NLXCwNCYT2YMw9zUiPTMHNIysknLzCYtPYtN2w8RG6v9dy8pDTE2s0FpZIqdjRVuTja4OFgR8yiFK9fvkJlSfMje1NqZV4b10DllUBWVmwShsogg/FhtD8JlURk96Kae7tjbWhEdm0B0bAKZ2TnPbYckQfvWzRg50B9HO2sc7W2wsTJHoVCUerj9WVXZgy5rr01nMpgEzRo44+poRUJyBvFJ6cQnpZGlpRqYNpJUUK3NQKlEpconNSVBa41kI1MrZFkmPy8HVW4meTlZqHIzdVahArB3cOCLeVNxdrDC9KkNNQqDvhhWFuoikZglaFWWRDDdG1vAgvcmqO8vyzLJqelExSaw+qftXA8O1dhrk2W4cPU2F64+2cHGQKnEwc4KlSpfYwAufMDdf51lqpZKYuVJICtL8O7XxRtbSyN++v3Qk6Sl0f1o37KRzvvpTAZD4v2J/YoFqYzMHFb8eICLN8PIy80p1qNVGhjh4mBFi8au5OWpyM3LJzdPRZ5KRfC9SI0BGCArJRYpL6Ogl5ybW9ACCdycHUnPzCYpOUXr83B3dqChu32x409vCGBm5FipGwIIQnUjesKPiZ5wxanoHrRCIdGra1s6t/UmLj6JuIRk4uKTiUtIIvRhDLl52nt9kiThYGuFrbUltjYW2FpbYmdtCRJs3PGXxkSpyhr+Ls+w+aEzt1i5bv/jYJqLQmmIoakV/5g4QGsvet2O02zceZTMZM3Zw5qGeWVZ5p2P1nAz+K7Wtkig3ie7tXdDWjZriKWF2XPnhFd+PJN2LT213i6GlTX7+uuv2bt3L3v27NG6ld6uXbv45JNP2Lt3b7n22RWqnugJCxWuonvQyPDW2IEa7/+8SmJuzvZYWZiRmJxKaHh0kfW72uTLMjM+WkWj+q5YWphhZWGm/pqbp+KHTXtLvR90eFQcK9ZsLui1l2Ef6ZzMFNLjwx7fTUaSJHIyksjN1N7zbNXYSWMABshMjsHJ2ojj564RFhFb8C+y4F9auvbtCAE6t2vBsjmTix1v36oxLw7oyZ/7i28V+OKAnjoDMJR+uVhdEBkZybfffsucOXN07mU7YMAAli5dypdffsmiRYuqpG01YSWDLikpKXz00UccO3ZM4y5KzwoPD2fRokUEBgaiVCrp3r078+fPL/cuSiIIC5Winqujzg0lNJ2vq6a2tj9uncPfksTyf79dZPg7IzObxORUPv/uDy4H3dW6RCYrK4fge+HPDUhPy5dl3vznf7C3tcLUxAgzExNMTYwwNTUhLDJGvcH9s2Tgm/V/MrBXR5RKpXpu1sCg4F98YgrL12wu0tbC75et2YyxsRH2tlaoHg8n56lU5OWp2H/0grrTrcnS/24o8rOttQWN67uSnplFSGik1oSuxvVdtP4ffDB5OL27tuHH3w+qh9zfGP3CcwOwoNm6deuQZZlBgwaRmprKggULOHnyJOnp6bi7uzN9+nSGDx+OoaEhAwYM4Pfff+e9996r9N5wTaoFoM2iRYvIycnh+PHjREREMGnSJDw9PenZs6fG8+fPn4+9vT3Hjh0jOzubmTNn8tVXXzF37txytUMEYaHaKEsPujTBW5IkzM1MMDczoXljDwKvh2gpYqLgpUH+TBk/FJUqn7SMTFLTMkhJy+D/ftnF1ZshWgObiYkRZqYmZGZlk5iURmZW9nOTz2RZ5uSF65y8cL1k/1HP3HfhyvWlvh9AA3cnXhneh/ruTtR3c8LSoqC+c2Gym7bNLQb36azzuu1aetKuZe3deq4qHTx4kFatWmFnZ8f8+fPZvXs3gwYNwsPDg99++4158+YxcOBAjI2N6d69O7/88guHDx9m1KhRldam8o7qlNb+/ftZs2YN3t7enD59mpycHN5//31Gjx5d5mtmZGSwb98+tm7dioWFBV5eXowZM4YtW7ZoDcLh4eFMnDgRExMTTExMGDBgAAcOHChzGwqJICxUK6XtQUNlDH/L6kCjVCqwtjTH2tIcgFbNGnDt1n1kWdPctYKhfToXa39+fj5rft7JH3uOaclwlujWoSW9uvqgUuWTp1KRm6dClaciNy+PA8cucj8sSkvPFBo3cGdQ747q3rNSocDAQMnh01c4deG61g8a/h1baey1lHVUQqhYMTExRERE4OfnB8Ds2bMZM2YMDRs25MaNG5w6dYqgoCASExNxcXHBy8sLgEuXLpUpCP/0+34Ol6BiXmJyqtZkyHxZZvrcr7C1ttR6/95+bZk0ekCJ2xUcHMzdu3eZOHEiS5YsYc+ePcyZM4eBAwdiafnkcaZOncrFixc1XsPV1ZWdO3eqfw4NDUWWZZo1a6Y+5u3trTOoTpw4kZ07d9KhQweys7PZt28fvXr1KvHz0EYEYaFWqA7D308H76cpFAqG9/fjj93FK0lBQe/ynQnDtD5maloGoeExyBqXfSno3NaLlwf3KHZb8yb1OX3huuZBcC1tLVSWDzZCxYqNLZjPt7W1BeD+/fssWrSI27dv4+3tTWpqQZGcvMeJiYVD0IX3qyx5ebqXoj3v9tIKDg5mzJgxjBxZUCZ36NChzJ07l/DwcJo3b64+b+3atSW+ZkZGRrG5XCsrK9LTtRfJ6dSpE1u2bKFDhw7k5+fj7+/PhAkTSvlsihNBWKizKnv4uyLuB2UL/OV9zML7l3ZUQqg4hSMYhUF2+vTppKamcuzYMaysrHjnnXd48OBBsfspNFQWK4lJoweUqIeqeyWDghEDulXo6yY4OLhIsEtNTSUrK6tc895mZmbFAm5qairm5uYaz1epVEyePJlRo0axceNGcnNz+fTTT5k1axZfffVVmdsBIggLdVxVDX+X537lCaaiR1tzubq6AhAXV7CFY35+Punp6SxZsgRZljl8uOCDWeF67YSEBKBgH+XKVNYPhWWRnp5OeHg4jo5PXq/79u2jRYsWODs7Fzl38uTJWoej3dzc2L17t/rnhg0bAnDnzh2aNm0KwM2bN9XfPys5OZno6Ghee+01jI2NMTY25pVXXmHixInleXqACMKCUCZl7SWW9X7lCaaiR1szOTo6Ur9+fe7duwfAZ599xieffMKePXto2LAhI0aMYNu2bVy/fh1PT0+Cgwu2dfTx8anUdlVlzsDt27dRKBTs3LmTmTNncunSJb788ktWrFhR7Nzvv/++xNc1MzNjwIABrFy5kuXLlxMZGckff/zBkiVLNJ5vZ2eHh4cHv/76K1OmTCEvL4/ffvtNPQ9fHqJYx2OiWIcgCNXNF198wXfffcepU6fUc8PaLFu2jI0bN3LixIlyr10tiapYJ7xp0ybOnz+v7vm7u7vzwQcf0Lt373JfOyUlhXnz5nH8+HGN64QnT55Mhw4dmDZtGgC3bt3is88+4+bNm0iShI+PD/PmzaN+/frlaocIwo+JICwIQnUTFRVF//79WbBgAS+//LLW82RZ5oUXXqB79+4sWLCgCltYuRYuXIirqytTptTeJW9lm8EXBEEQKp2rqytvvvkmv/zyi87zAgICSElJYebMmVXUsqoRHBxM48aN9d2MSlWtg3BQUBCjR49m8ODBjB49Wj038qzIyEimTp3K4MGDGT58OBs3bqzilgqCIFSOf/zjH2zfvl3nOf369ePcuXO1rm707du3adKkib6bUamqbWJWTk4OM2bMYPHixfj7+xMQEMCMGTPYvXt3sRqqH374Ie3atWPt2rWkpaXx2muv4e7uTo8exddOCoIgCDXDhQsX9N2ESldtg/C1a9dQKpX4+/sD0LdvXxYvXszVq1eLZf8FBQWxcOFCACwsLOjatSt79+7VGIRjY2PVKf9Pu3nzJgAhISEV/VQEQRD0qlGjRpiamuq7GYIGeg/CR48e5Z133il2fPr06ep1coWcnZ2JjIwsFoTbtGnD9u3bee+990hKSuL48ePF1pAV2rx5M6tXr9banlmzZpXhWQiCIFRfIuG0+tJ7EO7Zsyc3btwodnzXrl2cPXu22HFN1WCWLVvGZ599xosvvki9evXo0qUL4eHhGh9v7Nix9OnTp9jxlJQUQkJCaNGiBcbGxqV6DiEhIcyaNYsVK1bU+iQCofzE60UorfK+Zho1alQJrRIqgt6DsDZubm7ExMQUORYbG1usdwyQlZXFp59+ql4b98knn6grojzLyckJJycnjbd17Vq+vUwbN24sPm0KJSZeL0JpiddM7VNts6PbtGlDdnY2p06dAuDIkSNIkqTxBbh69WrWry/Yzu3Bgwfs2bOH4cOHV2l7BUEQBKG0qm0QNjAwYM2aNXz11VcMGTKE1atXs2rVKpRKJQBvv/02AQEBQME87vnz5xk2bBjvvvsuCxYsKLK7hiAIgiBUR9V2OBqgRYsWbN68WeNt3333nfp7Z2dnfvzxx6pqliAIgiBUiGrbExYEQRCE2k4E4Qrg6OjIjBkzimy3JQjaiNeLUFriNVN7iQ0cBEEQBEFPRE9YEARBEPREBGFBEARB0BMRhAVBEARBT0QQLqHjx48zfPhwBg4cyBtvvEFsbGyxcxITE5k2bRqDBw9mwIABHDp0SA8tFaqLkrxmzp8/T7t27Rg+fLj63927d/XQWqE6+fLLL5k7d67G28T7TC0jC88VHx8vd+rUSb5165Ysy7K8bt06+c033yx23t///nf5888/l2VZlh8+fCj7+fnJkZGRVdpWoXoo6Wvmm2++kZcuXVrVzROqqYcPH8rTp0+X27RpI//73//WeI54n6ldRE+4BE6cOIGXlxdeXl4AjBs3jrNnzxbZEjEvL4/Dhw8zbtw4AOrVq4e/vz87d+7US5sF/SrJawbg4sWLBAUFMWrUKF5++WX27dunj+YK1cTmzZvx8/PjjTfe0Hi7eJ+pfap1xazqIjo6usjGEUZGRtja2hIVFaVet5eYmEhWVhYuLi7q81xcXIiKiqry9gr6V5LXDIClpSXDhg1TD0NPmDABZ2dnfH199dFsQc8++OADAFatWqXxdvE+U/uIIFwCsiwjSVKx409vqyg/Xm797Hmatl4Uar+SvGYAVq5cqf6+SZMmDB48mEOHDokgLGgk3mdqH/GbK4Fnt1XMyckhMTERNzc39TF7e3uMjY2LJN/ExMRo3HpRqP1K8ppJT0/nm2++IS8vT31MlmUMDMRnY0Ez8T5T+4ggXALdunXjxo0b3L59G4Dff/8dHx8f7Ozs1OcolUr69u3Lr7/+CkBERATHjx+nX79+emmzoF8lec2YmZmxZcsWtm/fDsDDhw/Zt28fgwYN0keThRpAvM/UPqJsZQmdOnWKFStWkJ2djb29PUuXLsXd3Z3hw4ezePFiWrduTUJCAvPnzyc0NJS8vDzeffddhg0bpu+mC3pSktfMrVu3+OSTT0hJSUGlUjFjxgwGDx6s76YLerZq1Sqio6P59NNPAcT7TC0mgrAgCIIg6IkYjhYEQRAEPRFBWBAEQRD0RARhQRAEQdATEYQFQRAEQU9EEBYEQRAEPRFBWBAEQRD0RARhQRAEQdATEYQFQRAEQU9EEBYEQRAEPRFBWBCq0MGDB3nppZdo164do0aNIjo6milTprB79259N00QBD0QZSsFoYrs2LGDpUuXsnbtWry9vZk8eTKGhoakpaWxefNmfTdPEAQ9ED1hQagCKpWKFStWMGvWLNq0aYORkRFt27blxIkTzJ49W9/NEwRBT0QQFoQqEBgYSEJCQpFtCnNycujfvz++vr56bJkgCPokgrAgVIHo6Gisra0xNTUFICkpiW3bttG0aVM9t0wQBH0SQVgQqoCrqysJCQmcP3+ejIwMZs+ejaGhIZGRkfpumiAIeiSCsCBUAV9fX9566y3ee+89/P39adCgAatXr+bIkSMsW7ZM380TBEFPRHa0IAiCIOiJ6AkLgiAIgp6IICwIgiAIeiKCsCAIgiDoiQjCgiAIgqAnIggLgiAIgp6IICwIgiAIeiKCsCAIgiDoiQjCgiAIgqAnIggLgiAIgp6IICwIgiAIeiKCsCAIgiDoyf8DIZ3LNlaLqbcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 400x270 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(3.4, 2.35))\n",
    "\n",
    "# set gradient palette (might need to re-run to make it effective)\n",
    "palette = sb.dark_palette(\"#69d\", reverse=True)\n",
    "sb.set_palette(palette)\n",
    "\n",
    "# plot curves\n",
    "for i, p_1 in enumerate(p_1s):\n",
    "    ax.plot(alphas, lyaps_HG_2[i, :, 1], \"o-\", label=f\"$p={p_1:.1f}$\", ms=5)\n",
    "\n",
    "ax.set_ylabel(r\"$\\lambda_2$\")\n",
    "ax.set_xlabel(r\"$\\alpha$\")\n",
    "\n",
    "ax.set_yticks([-0.5, -0.7, -0.9])\n",
    "ax.set_xticks([0, 0.5, 1])\n",
    "\n",
    "ax.legend(frameon=False, loc=\"center left\", bbox_to_anchor=(1, 0.2))\n",
    "\n",
    "# add inset with bound\n",
    "# Create inset of width 30% and height 40% of the parent axes' bounding box\n",
    "# at the lower left corner (loc=3)\n",
    "axins = inset_axes(ax, width=\"40%\", height=\"40%\")\n",
    "k = 3\n",
    "axins.plot(\n",
    "    alphas, lyaps_HG_2[k, :, 1], \"o-\", c=f\"C{k}\", label=f\"$p={p_1s[k]:.1f}$\", ms=5\n",
    ")\n",
    "axins.plot(alphas, bound, \"--\", label=\"bound\", c=\"grey\")\n",
    "\n",
    "axins.set_ylabel(r\"$\\lambda_2$\")\n",
    "axins.set_xlabel(r\"$\\alpha$\")\n",
    "# axins.set_yticklabels([])\n",
    "\n",
    "axins.legend(frameon=False, loc=\"center left\", bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "sb.despine()\n",
    "\n",
    "ax.text(\n",
    "    0.99, 0.02, \"(a)\", transform=ax.transAxes, va=\"bottom\", ha=\"right\", weight=\"bold\"\n",
    ")\n",
    "axins.text(\n",
    "    0.015, 0.03, \"(b)\", transform=axins.transAxes, va=\"bottom\", ha=\"left\", weight=\"bold\"\n",
    ")\n",
    "\n",
    "fig_name = f\"figure_2_phase_diagram_lines_p2_{p_2}\"\n",
    "# plt.savefig(f\"{fig_name}.pdf\", dpi=250, bbox_inches=\"tight\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50cb0782",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_xgi",
   "language": "python",
   "name": "venv_xgi"
  },
  "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.9.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "78db47bee367e3fae825d59af43d695909ce208402094df58551906768c6e91c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}