tutorials/case_study_zhang2022.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "417f4911",
   "metadata": {},
   "source": [
    "# Case study 1"
   ]
  },
  {
   "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": 1,
   "id": "8da1197d",
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "\n",
    "sb.set_context(\"paper\")\n",
    "sb.set_theme(style=\"ticks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "af5bece2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# utility functions to compute the Lyapunov exponents from the Laplacians\n",
    "def compute_eigenvalues(H, order, weight, rescale_per_node=True):\n",
    "    \"\"\"Returns the Lyapunov exponents of corresponding to the Laplacian of order d.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    HG : xgi.HyperGraph\n",
    "        Hypergraph\n",
    "    order : int\n",
    "        Order to consider.\n",
    "    weight: float\n",
    "        Weight, i.e coupling strenght gamma in [1]_.\n",
    "    rescale_per_node: bool, (default=True)\n",
    "        Whether to rescale each Laplacian of order d by d (per node).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    lyap : array\n",
    "        Array of dim (N,) with unsorted Lyapunov exponents\n",
    "    \"\"\"\n",
    "\n",
    "    # compute Laplacian\n",
    "    L = xgi.laplacian(H, order, rescale_per_node=rescale_per_node)\n",
    "    K = xgi.degree_matrix(H, order)\n",
    "\n",
    "    # compute eigenvalues\n",
    "    eivals, _ = np.linalg.eig(L)\n",
    "    lyap = -(weight / np.mean(K)) * eivals\n",
    "    return lyap\n",
    "\n",
    "\n",
    "def compute_eigenvalues_multi(H, orders, weights, rescale_per_node=True):\n",
    "    \"\"\"Returns the Lyapunov exponents of corresponding to the muliotder Laplacian.\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=True)\n",
    "        Whether to rescale each Laplacian of order d by d (per node).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    lyap : array\n",
    "        Array of dim (N,) with unsorted Lyapunov exponents\n",
    "    \"\"\"\n",
    "\n",
    "    # compute multiorder Laplacian\n",
    "    L_multi = xgi.multiorder_laplacian(\n",
    "        H, orders, weights, rescale_per_node=rescale_per_node\n",
    "    )\n",
    "\n",
    "    # compute eigenvalues\n",
    "    eivals_multi, _ = np.linalg.eig(L_multi)\n",
    "    lyap_multi = -eivals_multi\n",
    "    return lyap_multi"
   ]
  },
  {
   "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 = 50  # 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",
    "alpha = 0.5  # ratio between coupling strength of 1st and 2nd order interaction (must be in [0,1])\n",
    "# rescale = True # whether to rescale\n",
    "\n",
    "n_repetitions = 3  # 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 = [xgi.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_1 = compute_eigenvalues(HG, order=1, weight=1 - alpha)\n",
    "        lyap_2 = compute_eigenvalues(HG, order=2, weight=alpha)\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_1 = compute_eigenvalues(MSC, order=1, weight=1 - alpha)\n",
    "        lyap_2 = compute_eigenvalues(MSC, order=2, weight=alpha)\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": "iVBORw0KGgoAAAANSUhEUgAAASAAAADECAYAAAA/IfygAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvuklEQVR4nO3deVxU5f7A8c8wzIAwbCrIKqCm4i4uqImiZYpLbnVLcWnz57Us06Q09Zql15vmUm5ttyQl0nIhF8y0UG+aoqWoKbiyDyCyCLIMw/z+ICanGWRxYEZ43q8XL+GcM+d8YeTL83zPeZ5HotFoNAiCIJiAhakDEASh8RIJSBAEkxEJSBAEkxEJSBAEkxEJSBAEkxEJSBAEkzGbBLRnzx6GDx/OkCFDCA8P19t/6dIlxo8fz9ChQ1mwYAGlpaXVPndpaSnJyck1eo0gCHXPLBJQeno6a9as4euvvyYyMpJt27Zx9epVnWNCQ0NZtGgRP/zwAxqNhu3bt1f7/EqlksceewylUmns0AVBeABmkYCOHz9Onz59cHR0xMbGhqFDh3LgwAHt/pSUFIqKiujWrRsA48aN09kvCMLDydLUAQBkZGTg7Oys/drFxYXY2NhK9zs7O5Oenm7wXHl5eeTl5elsEy0fQTBPZpGADI0GkUgk1d5/r7CwMNavX2+84ARBqDNmkYBatGjB6dOntV9nZGTg4uKis//WrVvarzMzM3X232vq1KmMHTtWZ5tSqSQkJMTIUQuC8KDMogbUr18/Tpw4we3btyksLOTgwYMMGDBAu9/DwwMrKyvOnDkDwO7du3X238ve3h5PT0+dD1dX13r5PgRBqBmzaQHNnj2bKVOmoFKpeOqpp+jSpQvTpk3jtddeo3PnznzwwQcsXLiQgoICOnTowJQpU4weR3R0tPbzoKAgo59fEARdksYwHUdycjKPPfYYhw8fxtPTs9LjlixZov188eLF9RGaIDRqZtECEho30fJsvEQCEkzuyJEj2s9FAmpczKIILQh14eTJk0yePNnUYTyQdevWsW7dOlOHUWdEC0gwqegzScTe6UCJRoZcoiL6TBJBPbxMHZZQT0QCEkwm+kwS6789R4lGDkCJRs76b88BGC0J3b59m2nTppGYmIivry+PPPIIEomEOXPmADB//nwCAwM5evQoEomE+Ph48vPzmTFjBmPGjKGgoIB3332XK1euoFarmTZtGiNHjmTnzp3s2rWLnJwcBg0axMSJE5k7dy65ubm0bduWmJgYjh49yrp16zh79ixpaWmEhITwyCOPsGbNGoqKisjNzSU0NJTg4GDmzZtn8PoAsbGxPPvss6SnpzNu3DheffVVo/xszIFIQEKd+Ol0Ij+eSrzvMXEJ2ahKy3S2FavUfLT9LD+cTKj0dUN6t2Rwz5bViiM1NZWPP/4YDw8P/vGPfxAYGMh///tfZs+eTWFhISdOnGDJkiUcPXqU9PR0vvnmG7Kyshg3bhyPPvooYWFhdOzYkffff5/8/HyeffZZunbtCpQPot6/fz+Wlpa8+uqrBAcHExISwo8//sjevXu1MZSUlLB//34AXnvtNZYuXUrr1q05ceIE//73vwkODtae7+/XB8jKyuKbb74hPz+fwYMH8/zzz6NQKKr1/Zs7kYAEk/l78qlqe220b98eL6/y1lTr1q2xsbHBw8ODmJgYUlNTGThwIHJ5eQts3LhxyGQyXF1d8ff358yZMxw/fpyioiJ27NgBwN27d7ly5QoAHTp0wNKy/Ffol19+Yfny5QAMGTIEe3t7bQxdunTRfr5y5Up+/vlnDhw4wLlz5ygoKNDuM3R9gMDAQORyOU2bNsXJyYnc3FyRgAThfgb3rLqV8sLSg2RmF+ptd3ZqwvKX+xsljooEAeXjBzUaDePHj2fv3r2kpqbqdGekUqn287KyMiwtLSkrK2PlypV07NgRgFu3buHg4MCePXuwtrbWeW1lj9Tde9zEiRMJCAggICCAvn37Mnfu3Ptev7LvoaEQd8H+VFEMPZ3Xldg7HYg+k2TqkBq8KcF+WMmkOtusZFKmBPvV6XWHDRvGiRMnuHXrlrY7BRAVFYVGoyElJYXY2Fh69OhBnz59iIiIAMrHKD755JOkpaXpnbNfv37s2bMHKH+s4O8zMgDk5ORw8+ZNZs2axcCBA/nll19Qq9X3vX5DJ1pA1E8xVNBX8bNdF/Gr9i7YzKf71PnP3Nramm7dutG2bVud7UVFRYwfP56SkhLeffddnJycmDlzJu+88w4jR45ErVYTGhpKy5YtdQZPA7z99tu89dZbbN++nfbt2+t0wSo4Ojry9NNPM2LECBQKBd26daOoqIi7d+9Wev2GTgzF4P5dgS8WPlEfITZq9TkERqPRUFBQwDPPPMPmzZu180zNmzeP3r17M27cuFqd96uvvqJfv360adOGixcvsmjRInbu3Fnt1z/o9R9WogUE3DKQfO63XTCugQMH1tu1zp8/z0svvcQrr7yiM8ndg/L29mbOnDlYWFhgZWXFe++9Z7RzN2SiBUTlLaBmDtZs/tfQ+ghREBolUYTGcDEUoKi4lD9uZJkgIkFoHEQCorwYOvPprsglJYAGuaSEycF+2NtaMX/jL+yKvtqgbn0KgrkQNaA/BfXw4sjeL7Rf/+PxpxnxqC8fbvudL/Zc5I8bWcx61h9FE1m1zymmmRCE+xMJ6D5sm8iYP7UXkUevs3nvRWavieatKb1o4+lYrdeLaSaqRyTqxkt0waogkUgYM7A1y1/uT2lpGW+uO8aBEzdFl8yIjhw5ov0QGheRgKrJz7cpa+cE0alVMzZ8d47VEb9RVCyWen5QxxJOEd86g4vtlMS3zuBYwqk6vd60adMqXVOuJtq1awdARESE9klpQ6raDzB48GCSk5MfOKba2rlzJ/PmzTPJtUUXrAYcFFYsntaX7YfiiTh4mWvJucyf2guvFnamDu2hdCzhFJ/EhKOSlQ8+VcnK+CQmHIBA7951cs3PPvvMqOebMGHCA+1v7EQCqiGphYQJT7TDz8eJD8LPMGftEV55uhtB/pVPdt8YHbnxKz/fOH7fY+KzblBaptuKLFGXsOnUFg5f+1+lrxvk24+Bvn3ue26lUsncuXO5e/cuFhYWLFy4kG7dujF48GC++uorTp06RXR0NBkZGSiVSqZOnUpqaiq//vorjo6OfP7552RmZjJjxgy8vLxISEjA3d2dlStX4ujoqL1OxWyFr776Knv27GHTpk1IJBI6d+7Me++9x8cff6zdv3XrViIjIyksLEQikbB27Vpat25tMP7i4mKWLFnCmTNnkMlkvPzyywwfPpyzZ8+ybNkyiouLcXJy4t1338Xb25vJkyfj5+fHiRMnKCoqYuHChWzZsoWrV6/y3HPP8dxzz7Fu3Tpu3rxJYmIiOTk5PPPMM7z00ks6142NjWX58uUUFRXh5OTEkiVLcHJy4sknn2TZsmX07duXF198kcGDBxtlrT3RBaulbm1d+HBOEK08HFgVfoaNO86hKlVX/UJB6+/Jp6rtNfHdd98RFBTEzp07CQ0N1U5tca/z58/z+eefEx4ezn/+8x8GDBigHVB67NgxAOLj45k6dSr79u2jdevWla66m56ezvLly/niiy/Yt28farVap6aVn5/PoUOH2LJlC3v37uXxxx/n66+/rjT+LVu2cPfuXaKiovjyyy/ZsGEDJSUlzJkzh0WLFvH999/z7LPPaidWq7Bnzx5Gjx7N0qVLWbduHeHh4WzYsEG7Pz4+ns2bN7Nz5062bdvGxYsXtftKSkpYuHAhq1atYteuXTz//PMsWrQIhULBsmXLeOeddwgPD0cikRhtoU/RAnoAzRyasGzGo3y1/xK7oq9yJTGbt6b0wrWZralDM7mBvn2qbKW8vGcBt+7e1tve3KYp7wyeY+AV1de3b19effVVLl26xMCBA5k0aZLeMf7+/igUCu3cOn379gXKF8KsGM3u4+NDQEAAAGPGjNGZPuNev//+O/7+/tpFMFeuXAnApUuXAFAoFKxatYp9+/Zx8+ZNjh07hp9f5aP+Y2Ji+Mc//oGFhQXOzs7s27eP+Ph47O3ttfMLBQcH869//Ys7d+4AaBfrdHd3p2vXrjRp0kTnewEYOXIktrbl/z8HDx7Mr7/+qh30evPmTZKSkpgxY4b2+Pz8fO3Ppk+fPqxevZqoqKhK464p0QJ6QJZSC14Y1ZEFz/cm7VYBr685wskL+tM1CPomdBmNXCrX2SaXypnQZfQDn7tHjx7s27eP/v37s3//fv75z3/qHSOT6T7Tde+8O4a2aTQanTl77vfa27dvc/v2X8k1LS2NZ555hjt37jBgwADGjh173zupfz9fQkICZWX6E7VpNBrtlB73fj+GvhfQn3Po7197enoSGRlJZGQkO3fu1LbSNBoNN27coEmTJty8ebPSuGvKLBJQamoqISEhDBs2jBkzZujMElchIyODF198kdGjRzN27FhOnDhh9DgGDhyo/aipPp3cWDsnCNdmNiz98hRLPj/BOTG/0H0Fevdmeq8QZCoL0IBMZcH0XiFGKUCvWLGCyMhIxo4dy7/+9S/++OOPWp3nxo0b2lbMjh07Kl0SvHPnzpw7d47MzEwA/v3vf3P48GHt/vPnz+Pt7c1zzz1H165dOXr0qM5cQH/Xq1cv7fxAWVlZTJo0CQ8PD3JycoiNjQVg//79uLu769SkqnLo0CFKSkrIzc3l559/pn//vyZ+a9WqFbm5udqpRnbs2KFt8X399dfY2NiwceNGFi5cqJ1C5EGZRRdsyZIlTJw4kREjRrBhwwY2btxIaGiozjErVqxg0KBBTJo0ievXrzN58mSOHj1a6V+k2njQh+Bcm9myYmYgSz7/ldOXMgAxv1BVAr1789Pmv5r0gZOMc/dr8uTJvPHGG+zatQupVFrraT4cHBz46KOPSExMpF27dixdutTgcS1atGDBggW8+OKLlJWV0a1bN8aNG8fGjRsBePTRR4mIiGD48OHI5XK6dOmindrVkIkTJ7J06VKefPJJABYtWoSdnR1r1qzhvffeo7CwEAcHB9asWVOj78fKyoqJEyeSn5/P9OnTadOmjTahyeVyPvzwQ22RW6FQ8P7775OUlMSmTZv49ttvcXNzo3///qxYsYJ33nmnRtc2xOSj4VUqFQEBAZw6dQpLS0vS0tKYNGmSzl8PgIMHDxIQEICDgwNqtZqAgAB+/vln7Ox0b4Hn5eXpzUanVCoJCQmpcmlmYxHzC9WMuS6JnZyczJQpU/jpp59MHYpR3HvHzlyYvAWUnZ2NQqHQ9lmdnZ0NPij2xBN//eL+97//xc/PTy/5AISFhVV6p6K+iPmFaqY+5wMSzEu9JqCoqCjtygEVfHx89I6TSCSVnmPz5s1s27aNrVu3Gtw/depUxo4dq7OtogVUX5o7NTHYArK0tCD7ThFOdtYGXtV4mev4L09PzwbT+gHzavlUqNcEFBwcrF0DqUJFF0ytViOVSsnMzMTFxcXg61esWMGRI0cIDw/X3u78O3t7e4Pz8danKcF+rP/2HMWqv4qMllIJZeoyZq2KJnRyTzq3bm7CCAXBPJj8LphMJqNnz57ahdt2795t8E7D5s2bOXnyJBEREZUmH3NhaH6hWc90Z82cIGysLVm46Re2H4qnrEwMaBUaN5MXoQFSUlKYN28eWVlZuLm5sXr1ahwcHIiIiCAjI4PXXnuN3r17o1AodFo3n376KS1atKjy/FVNyVpXDBVX7xap2PDtOY6eTcG/vQtzJvjjoLCqt5gEwZyYvAgN5U+ebtmyRW/7vQP5YmJi6jOkOmNjLWPupB50at2MT3df4PXV5V2yDr7NqnytmDdHaGhM3gVrjCQSCcH9fFn5WiCWlhbM3/gLO3+uetpXMW+O0NCYRQuosWrj6cja2UF8tP13vtx7kYvXs3h9QnfsbORVvlYQTMlYrXHRAjIx2yYy5k3pxf+N6cxvcem8vjqa+MRsU4clCPdlrNa4SEBmQCKRMCqwFe/PDATgrfXH+P7YNTHtq9DgiQRkRtq2dGLtnCD827Xgs90X+M9XMRQUqkwdliDUGVEDqkO1GWJgZyNn4Qu92X3kGpv3/cH1lJqtxCEIDxORgOpQbYtzEomEsUFtaO/dlBVbYgj96BhB/h7E3ulAiUaGXKIi+kySGFkvVIs5P74humBmrGIlDg8XWw7FJFGikQMS7fQeYo4hoTrM+fENkYDMnIPCiruF+nMkF6vUfBV1yQQRCYLxiAT0ELiVI6b3EBomkYAeAs2dmhjcLpdZkC/ukgkPMZGAHgJTgv2wkulOPSuVSihRlTFLPLjY4ERHR2s/GjqRgB4Chqb3eP2Z7qx4LRA0Gt5cd4xd0VfF9B4NhDkXjY1NJKCHRFAPL7rY/UFP+3N0sfuDoB5etPduyodzgujd0ZUv9lzkvS9OkptfbOpQBaHaRAJ6yCls5Myf2ot/ju3M2fhMZq2O5sK1W6YOSxCqRSSgBkAikTCifys+eC0QuUzKgk2/sO3HONSiSyaYOfEkdAPS2tORtbMHsmlHLFsPXCb26i3eCOlBU/vKJ8E356dkHxbiZ1h7IgE1MDbWMuZM9KfrI83ZtPM8s1ZFM3uiP/7tDE/0f2+hU/zy1I74Gdae6II1QBKJhMd7e7P69QHYK+Qs/vQEYfv+oFStv7a4INTUsYRTxLfO4GI7JfGtMziWcKrW5xItoIdITUfXe7vas2rWAD6PvMB3P13h4vUs5k7qgYuTTR1FKDR0xxJO8UlMOCpZ+R8zlayMT2LCgfJltmtKJKCHSG2a99ZyS2Y+3Y0ubZqz/ttzzFoVzevPdiegk5vxA3xIiJpNzZWVlZGan87m37ZToi7R2VeiLiEiNlIkIKFyA7p70sbLkRVbTrP0y1M8GdiKVu72jXKKD1Gzub8yTRnK/Eyu307g2u1ErmcnciM7kaLSyp8xy7p7u1bXEgmoEXFvrmDlq4F8ufcPvj92HYkENJryCfArpvgAGkUSaiwq6jUqyzJkpRYcSzil01LRaDSk52dyLTuB638mm+vZiRSqigCQSWX4OHoS5NOXVk1bEhEbSXZRrt51mtk0rVV8ZpGAUlNTCQ0NJSsrC19fXz744ANsbW0NHpufn8+YMWNYtmwZAQEB9Rzpw09mKeX/xnTmyG/J5BXoNqUrpvgQCahhMFSv+ThmK1eybiCXyrl+O4Eb2YkUqMpnVZBZWOLt6Emgd29aOXnTumlLPO3dkFr8NQ5RaiHlk5hwnW6YXCpnQpfRtYrRLBLQkiVLmDhxIiNGjGDDhg1s3LiR0NBQg8e+99575OXl1XOEDc+dvyWfCmKKj4ZBo9Gw9dxOvXqNSq3iwJVopBZSvB086NuyJ62dWtKqqTde9m5YSu+fEipaTx8fC9O2qqb3CalV/QfMIAGpVCpiYmLYsGEDAOPGjWPSpEkGE9D+/fuxtbWlXbt29R1mg9PcqQmZBpKNg515LRMtCsZVq+hGlXefkriRncD17CQKSu5W+pqvxq1BJpXV6nqB3r35aXPUX19Pql3yATNIQNnZ2SgUCiwty0NxdnYmPT1d77jU1FTCwsIICwtj2rRplZ4vLy9Pr4WkVCqNG3QDMCXYj/XfnqNYpdbZnnOnmG9+jOPpx9oitZCYKLq/NLaCcVU1mzJNGen5t7h+T83mRnYSd//sRlW0bPp4+nMy+TfyDSSh5jZNa518jK1eE1BUVBTLly/X2ebj46N3nESi+x+/rKyMBQsWsGjRIqytKx9WABAWFsb69esfONaGrqLOsy7iV+1dsP8b34sL17MIP3CZ3+MyeGNiD1yaimeG6ovBms2prcRnXkMmlZUnm5wkbYHY0sISbwcP+rXsSSunlrRyaklLB3dtN6qDyyNGrdfUhRolILVazeHDh7GwsGDQoEFIpeXFqaioKIKDg6t8fXBwsN5xKpWKgIAA1Go1UqmUzMxMXFx0hw1cv36d69evs2DBAgASExNZuHAh7733Hn369NE5durUqYwdO1Znm1KpJCQkpCbfaqMQ1MOLI3u/0H49tO/TDO3rg397FzbtiOW1VT/z8lNdGdDd04RRNg7qMjVbzxqo2ZSp+OHa0b8KxC1706ppebLxrKJmY+x6TV2oUQJ688038fDwQCaT8fnnn/P+++/j7e1NREREtRKQITKZjJ49e7J//35GjRrF7t27GTBggM4xbdq00WmKT548mZkzZxq8C2Zvb4+9vX2tYhHKDerhhZ9PUz4IP8PKrWc4czmD6WM7Y2NddbNd1GyqVqouJSkvjRvZf3WhEnKSKVFXPr1u2Pi1WFpIK91fGWPWa+pCjRJQRkYGq1atAsqLxW+99RYzZ8584CAWL17MvHnz2LRpE25ubqxevRqAiIgIMjIymDVr1gNfQ6gZ12a2vP9Kf775MZ7th+K4dOM2b4T40877/s97iJqNbs2muLSEhJxkbmQncePPZJOYl4q6rLz21sTSGl8nLx5vHcixmye5U1Kgd43mNk1rlXweBjVKQCqVipKSEuRyOR4eHnz88cfMmTOH+Pj4BwrCw8ODLVu26G2fMGGCweMNHSsYn1RqQciw9nRr68yqr8/w5vr/MXFoO54abB4FalMzXLPZwtnUi0gkEm5kJ5J8R4lGUz4vk53cFl+nloxs+xi+Ti1p5eSFi6I5FpLyMeGtm3qbfc3G2GqUgN5++23y8vJo3rw5AAqFgk2bNrF37946CU6oe9UZ4NqxVTM+emMQm747x9aoy/wel8mcif6NelBrTlEeYb9/Z6BmU8qxxFM4NXHA16klAV7d8XVqia+jF81snPRusNzrYajZGFuNElCXLl30tkmlUkaPbrgZuqGrbjdJ0UTG3Ek96OHnwsc7Y3ltVTSvPNWVwG4edRugiZVpysgoyOJmdhI3spO4mZPEzexkg8MR7vXJk/+p1fXMvWZjbDVKQIWFhcTHxxMXF8fly5eJi4sjPDy8rmIT6kFNisYSiYTBPVvi59OMVeFnWLHlNL9dzuD/xnamiZXJHymrlvvVbErL1CTnpnEz555kk5Osve1tIbHA096Nzq7t8XH0IvLSD+QW39G7RvNajotqjO77v+batWscOHCAuLg44uLiSEpKQqFQ0K5dO/z8/HjqqafqK06hjtSmaOzW3Jb/zOxPxME4vj0cz8UbWcwN6UHblk51FKVxGKrZbDz5FYeu/o+i0iKS8tIoLStfBttKKtfe9vZx8sLXyQsvB3fk9zzA52Bt1+hqNsZ23wT04osv4uPjg7+/P3FxcfTq1Yt169aJ29wCllILJgf70b2tM6u+/o031x2jb2c3s5veQ6PRkF2Yy82cZP575hu9mo1ao+byrat0atGO4W0H4eNYnmzcFC5YWNx/wtDGWLMxtvsmIBsbG7788kskEgnTp09n9erVjB49miVLlug9qyM0Tp1aN2fdG0Es/vQE/zuXCphueo/SMjWpeUpu5iRzMyeZhJwkbuakcKc4v4pXalgUVLtHPRpbzcbY7puA9u/fr/3cysqK+fPnM2zYMN5++21++OEH5s+fj0KhqPMghXLm+pCfwkZO9h39yaoedHqP+9VrCkrukvBnormZk0xCdrJOF0pmYYmXgzu93Lvg7eiJj5MnH534kqxC/WWsazuXjfDgalw57N69O5GRkaxdu5YxY8Zw6NChuohLMMCcH/K7lWN4Go/aTu9RWb1m7+VD5KvuklmQpT3WwcoOHydPhrsOwsfRE29HT9ztWujMYwMwsesYUbMxM7W6dSGXy3nzzTcZNmyYseMRHlKVTe8hkUiIvZpJlzbOVZ6jUFVEUm4qCTkpbDm3w2C9JjE3hQAvf4a0DsTH0RMfR08cmzhUK0ZRszE/D3Tv1NBzQULjNCXYj/WH9oPbZSTyIjQl1mhS22FT5M3Cj48zZmAbJge3R2YpLX+2Jv8WCbkpJOSkkJiTQkJOMukFVS8pXaYp4/W+L9Y6TlGzMY6artBSmYfj4Q3B7EmbpyH3vUCpprwGI7EqwrLVRSZ29eP38xK+/+MwRzL24+xWSnphOsV/TnAukUhwU7jQqqk3g1r1o6WDB96OHiz+aTW3DEx0Luo15sFYJQCRgOqQuRaN60JE7G5t8qlQqinly7PfACD3gUK1jJupdrRv0YFB3Tri4+SJp70bVpZyvfNN6DJa1GsaAZGA6pA5F41rS6PRkFWYTVJuKok5qSTmppCYm8qtu/p3lyrMC3wFH0dPLNTWrP/2HCejlWiS4fVnXQwmHxD1GmMyVnepLogE1IhVNZVEQcnd8kSTm6JNNkm5qdpVFACaNnGkpYM7GZbWFJYW6V2juU1T/N07ab9e8HxvfjyVyGe7z/Pqqp+ZMa4LA/0NT3gm6jXGYc5//EQCaqQM3ebedOorfkmIQYOGxJxUnWdmbGRN8HJwp1/LnrR08KClozteDu4o5LY656uqyySRSHgiwJvOrZuz+uszfBB+hlN/KJkxrgsKG8OtIaHhEgnoHg29ZlNapkaZn0FybprBYQmlZWp+S7uAt6Mnfi6P0NLBXZtsmjUx7lQSbs1t+c8r/fnupyt8fTCOP65n8foEf7o+UvXt+obOnLtMxiYS0D3MuWZTVXfpXmVlZaQX3CIpN7X8Iy+N5Nw0Uu+ka58UrowEWDl0Qa1irGmXSSq14Jkh7ejezoXVX59h4cfHGT2gNVOG+yGXNcwZAKvD3P7v1SWRgB4ChrpLn8SEo9GU0b55G5Ly0rTJJjk3jeQ7SlT3zC/sbNsML3s3urt1xMvBHU97N1b+72OzGZbQtqUTa2cH8cXei0QevcbZ+AzeCOmBr3v1HjAUHl4iAT0EImJ363WXStQlrD8ZprOtWRMnPB3cGOrSFi8H9z+TjSvWMv2ljMxtWIK1lSUvj+9K7w6ufLjtd+asPUo/MxxdLxiXSEBmpKLrlJKXRnKesrw1k5d231vc/9czBC8HNzzt3bCVV3+KVHO9zd3TrwXr5w5i8WfHOXo2BVOOrq+uxlSzMTaRgOrIfWfeU5eizM8kuSLR5KWR8meNRnVPjaZZEyc87F2xtrSiqFR/tHlzm6Y83rp/rWM019vcDgor8vL1l6h50NH1daUx1WyMTSSgOmB4JHcY++N+okhdjPJOBmpNmfZ4F9tmeNq70cXVD097Nzwd3PCwc8VG3kTnfObSXaoPlY2uNzTgVXh4iQT0gDQaDXeK80m5oyQlL53UPCU/XD2Kqkz3L7haU8aNnCR6uHemt0c3POxd8bR3w92+BdaWVve9hrl2l+pS5aPr4XhsKv26uNf63KLLZD5EAvpTVbe51WVqMguySLmTTkqektQ8pfbz/HsWk5NLZXrJp4JGU0Zo/3/WKj5z7S7VlSnBfqz/9hzFKrV2m8zSAic7K5aHxRDUw5PpYzrX6uFF0WUyH2aRgFJTUwkNDSUrKwtfX18++OADbG1tdY4pKSlhxYoVnD59GpVKxfz58+nfv/b1j3tV9lTwicQzWFpYkpKXRlp+ps4zNA7W9njYtaCPlz8edi3wsHfFw96VZjZOzNy7SIzkfkAVdZ6N38ZQqLKgiayMl5/2p383D7YfimfboXhir9xi1jPd8W/vYuJohdoyiwS0ZMkSJk6cyIgRI9iwYQMbN24kNDRU55jPP/+c7Oxsdu3axdWrV3nhhRc4evTofZ/Ora6I2EiDTwWfTo3Fzc4FDztXurt31iYad7sWKKxsKzmbGMltLEE9vAwWnCcObU/vDq6sjviNxZ+dYFhfH14Y1fGhWRpI+IvJ3zGVSkVMTAwbNmwAytecnzRpkl4CioqKYuXKlUgkEh555BG+/PJLNBqNXgLKy8sjLy9PZ5tSqbxvDFkGWitQ/lTwh8OX1PA7apw1m/rWxsuRtbMHsiXqkvbhxdef9adjq2amDk2oAZMnoOzsbBQKBZaW5aE4OzuTnp6ud1xCQgIxMTG8/fbbyGQyZs+eTZs2bfSOCwsLY/369TWKoZlNU6N3mRpbzaZCfRZ45TIpLz7ZiYCOrqz95nfmb/wfYwa2YdKw9o16KMfDpF4TUFRUFMuXL9fZ5uPjo3ecoW6VWq1GqVSyY8cO4uLieOmll4iKisLOzk7nuKlTpzJ27FidbUqlkpCQkErjEl0m4zFFgbdT6+asmzuIL/ZcZFf0VU5fSmfOBH/aeDnWeyxCzdRrAgoODiY4OFhnm0qlIiAgALVajVQqJTMzExcX/aJi8+bNGTFiBBKJhPbt2+Pq6sqNGzf05qW2t7ev8cKJosv08GtiZckrT3WlTydXPtp2lrkfHeWZx9vy9ONtsZTef4FBwXRM/s7IZDJ69uypXYNs9+7dBhc9HDRokPaYpKQk0tLS8PX1NVocgd69aXvNhY5xrrS95tJoks/AgQO1Hw1Bj/Yt2BA6iMBuHnx9MI7Qj46SqMyr+oWCSZg8AQEsXryY7du3M3z4cE6fPs3rr78OQEREBB9++CEAc+fOJSMjgxEjRvDPf/6TpUuX6nW/hJoLCgrSfjQUChs5b4T0YN6UXmRkF/L6miPs/Pkq6jKNqUMT/sbkRWgADw8PtmzZord9woQJ2s8VCgUrVqyoz7DMTkNppdSXR7u606FVUzZ8e44v917k1B9K+nRy5ftj17mVXUhzpyZMCfYzu7FljYlZJCChehpSK6W+ONlZs+D53vx0OokN353l4vW/VlTNzC402xH2jYVZdMEEoS5JJBIe69USe1v9MXcVI+wF0xAtoDokukzm5Xau/qodUPv164UHJxJQHRJdJvNS2Qh7CwsJcQm3aectxurVN9EFu0dDuyUt6JoS7IfV356QtpRaYG0lJXTdMT7dfZ67RYZnMhDqhmgB3UO0WBq2ikLzV1GXdO6C9e7oylf7L7H3f9c5cT6Nl8d3oVcHVxNH2zhINBpNg384Ijk5mccee4zDhw/j6Wl4FU5BuHTjNuu/O0ui8g6B3TyYNqYTTnb6E/oLxiO6YILwJz/fpqydHUTIsPacOJ/GjPd/4uDJBBrB32iTEQlIEO4hs7Tg2SHt+OiNIHzc7Fm3/SwLNh0nJTPf1KE1SCIBCYIBXi3s+PeMR5n5dFeup+Tw6gc/s+1QHKrSsqpfLFSbSECCUAkLCwlD+/iw8a3H6N3Bla1Rl5m9Jpq4BMMT2Ak1J+6CCUIVmtpbM29qL05eSOPjnbGErjvGiEd9mRzsx6mLSr27amJYR/WJBCQI1RTQyY3ObZqzZf8l9v1yg+gzyRSXqFGpy7tlYmxZzYkumCDUgI21jOnjurDi1UAKi0u1yaeCGFtWMyIBCUIttPduSlkl8wuJsWXVJxKQINRSc6cmBrc3sbYkv1AM6agOkYAEoZYMjS2zkEi4W1TKS8t+ZNuhODG2rAqiCC0ItVTZ2DJvN3vCD1xma9Rlvj96nacGP8LwR331kpUgxoIJQp2JT8xma9Qlfo/PpKm9Ff94rC1P9PFGZikSUQWRgAShjl24doutBy5z8XoWzk5NeHZIOwb39BLLBSESkCDUC41Gw9n4TLYeuER8Yg5uzW2Z+EQ7Art7IrXQX4izsRA1IEGoBxKJhO7tXOjW1pmYP9LZeuASq77+je2HrxAyrD19O7lh0QgTkUhAglCPJBIJvTu60tOvBcfPpxJ+4DL/CYuhlYcDk4a1J/+uii0HGs/QDpGABMEELCwk9O/qQd/O7hz5LZmIg5d5978nkUigoijSGIZ2mEUVLDU1lZCQEIYNG8aMGTMoKCjQO6akpIQ33niDUaNGMXr0aI4fP26CSAXBuKQWEgb39GLTW49hZyPj7xXZhj60wywS0JIlS5g4cSIHDhygU6dObNy4Ue+YyMhIysrK2LNnDytWrGDevHkmiFQQ6oal1IL8u4YfWszMLmTbj3EkKPMa3OyMJu+CqVQqYmJi2LBhAwDjxo1j0qRJhIaG6hxXVlZGYWEharWawsJCrK0Nz9Wbl5dHXl6ezjalUlk3wQuCEVW2bJCl1IKtBy6z9cBlPJxt6dvZnb6d3XjEyxGJxDSF6+joaO3nD7KYg8kTUHZ2NgqFAkvL8lCcnZ1JT0/XO27s2LHs2rWLwMBA8vLyWL16tcHzhYWFsX79+jqNWRDqwpRgP9Z/e45ilVq7zUomZebTXencpjknLyo5cT6NXdFX+e6nKzR3sKZPZzf6dXang29TpPX0XFH0mSTW7c2kRCNDLlGBXVKta1T1moCioqJYvny5zjYfHx+94wxl9fXr19OtWzciIiK4efMmzz33HB07dsTDw0PnuKlTpzJ27FidbUqlkpCQkAf/BgShDlU2tKNi+/B+vgzv50v+3RJO/ZHOifOpHPw1gb3/u4GdjZw+nVzp09mNbo84I/9z2Ef0mSSjTpgWfSaJ9d+eo0QjB6BEI3+gQnm9JqDg4GCCg4N1tqlUKgICAlCr1UilUjIzM3FxcdF77eHDh1mzZg0SiQRfX1+6du1KbGysXgKyt7fH3t6+Tr8PQagrQT28qvxFVtjIGdzTi8E9vSgqLuVMXAa/nk/jl9hUfjyVSBMrKT39XHGwlXPwVAIlqqonTFOVqsnNLyE3v5jcghLy/vw3N7+YvD//zc0vIT4xG/XfpiGpKJSbfQIyRCaT0bNnT/bv38+oUaPYvXs3AwYM0Duuffv2HDp0iLZt23L79m0uXLjAnDlzTBCxIJgPaytLHu3izqNd3FGVlnH+6i2On0/l5AUlOfnFescXq9Rs+O4cR8+mkJdfQm5BeWIpLC41eH4LCwn2tnLsbeU42FrpJZ8KtZ0DySyGYqSkpDBv3jyysrJwc3Nj9erVODg4EBERQUZGBrNmzeLWrVssWrSIxMRELCwsmD59OiNHjqzW+cVQDKGxUZdpGBP6faX7W3k44GArx97WCgeFHHtFeYJxUPy1zUFhha21TOcJ7ReWHjRYKHd2asIXC5+ocZxmkYDqmkhAQmNk7GQBf9WADBXKa9MFM4vngARBMD5DE6ZZyaRMCfar9TmDengx8+muyCUlgAa5pKTWyQfMoAYkCELdqOqu2gOd9841vevUhkhAgtCAVeeuWq3O+wAPH95LdMEEQTAZkYAEQTAZkYAEQTCZRlEDUqvLbxmKQamCYBqurq7a8Z73ahQJKDMzE0CMBxMEE6nsGbxG8SBiUVERFy5cwNnZGam08iVRKgathoeH4+rqWo8RCpUR74l5qun70qhbQNbW1vTs2bPax7u6uoonps2MeE/M04O+L6IILQiCyYgEJAiCyYgEJAiCyYgEdA97e3tmzpwpJjQzI+I9MU/Gel8axV0wQRDMk2gBCYJgMiIBCYJgMo02Ae3Zs4fhw4czZMgQwsPD9fZfunSJ8ePHM3ToUBYsWEBpqeE5cwXjqeo9Wb9+PYMGDWL06NGMHj3a4DGC8eXn5zNy5EiSk5P19j3w74mmEVIqlZpBgwZpsrOzNQUFBZpRo0Zprly5onPMiBEjNL///rtGo9Fo5s+frwkPDzdBpI1Hdd6T6dOna3777TcTRdg4nT17VjNy5EhNx44dNUlJSXr7H/T3pFG2gI4fP06fPn1wdHTExsaGoUOHcuDAAe3+lJQUioqK6NatG1C+Wuu9+wXjq+o9Abhw4QKfffYZo0aN4t1336W4WH/VB8G4tm/fzuLFiw0ulWWM35NGmYAyMjJwdnbWfu3i4qKzGuvf91e2WqtgPFW9JwUFBfj5+fHWW2+xa9cu8vLy2LhxoylCbVSWLVtW6TAmY/yeNMoEpDHw5MG9q7FWtV8wvqp+5ra2tnz22Wd4e3tjaWnJCy+8wJEjR+ozROFvjPF70igTUIsWLbh165b264yMDJ0m5t/3V7Zaq2A8Vb0nqampfPfdd9qvNRqNwdHVQv0xxu9Jo0xA/fr148SJE9y+fZvCwkIOHjyosxqrh4cHVlZWnDlzBqDS1VoF46nqPbG2tmblypUkJSWh0WgIDw9nyJAhJoxYMMrvyYPVyB9e33//vWbEiBGaJ554QvPpp59qNBqN5qWXXtLExsZqNBqN5tKlS5rx48drhg0bppkzZ46muLjYlOE2ClW9JwcOHNDunzdvnnhP6tGgQYO0d8GM+XsihmIIgmAyjbILJgiCeRAJSBAEkxEJSBAEkxEJSBAEkxEJSBAEkxEJSBAEkxEJSBAEkxHPsgtmqbS0lE8++YQdO3ZQUFDAwoULUSqVlJaWMmPGDFOHJxiJSECCWVq7di0XLlwgMjKSmJgYVq5ciYWFBdu2bTN1aIIRiQQkmJ38/HzCwsLYv38/dnZ2dO3alevXrzN79mwUCoWpwxOMSNSABLPz66+/4uPjg5eXFwAqlQo7OzsmTZpk4sgEYxMJSDA76enpOtM6bNu2jRYtWojWTwMkumCC2XF1deXy5ctkZGSQlpZGZGQkd+/epaSkBLlcburwBCMSo+EFs1NSUsLChQv56aefcHBwYN26daxcuZKioiIiIiJMHZ5gRCIBCYJgMqIGJAiCyYgEJAiCyYgEJAiCyYgEJAiCyYgEJAiCyYgEJAiCyYgEJAiCyYgEJAiCyfw/8iSbZsWfhIgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 288x194.4 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(4, 2.7))\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=\"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=\"simplicial 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\"lambda2_HG_SC_N_{N}_ps_{ps}_nrep_{n_repetitions}\"\n",
    "# plt.savefig(f\"{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_1 = compute_eigenvalues(HG, order=1, weight=1 - alpha)\n",
    "        lyap_2 = compute_eigenvalues(HG, order=2, weight=alpha)\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": 12,
   "id": "b04ab8c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAADDCAYAAACS2+oqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABMGklEQVR4nO2dd3iUZfq2z+npCSU9gQRCCb0HUDqoKEUQRaoFFd1lLftbXFTWhi7rqny6IGJZEVlEEBSkV0GkBIy0dBICpHfSp8/3R8yQyUwggSQzSZ7zODgO8s5b7mnXPO/9XM99S0wmkwmBQCAQNHuk9g5AIBAIBA2DEHSBQCBoIQhBFwgEghaCEHSBQCBoIQhBFwgEghZCqxB0vV5PWloaer3e3qEIBAJBo9EqBD0rK4tx48aRlZVl71AEAoGg0WgVgi4QCAStASHoAoFA0EKQ2zsAR8FoNBGdWk5qvpbgdkp6BbsglUrsHZZAIBDUGSHoVIr5x7szScnRoNGbUMklhPqoeOF+fyHqLQC9Xk9WVhZ+fn7I5eIjL2i5iJQLEJ1abhZzAI3eREqOhujUcjtHJmgIxKR4w2IwGDkRFcO6Lfs5ERWDwWC0d0iCPxDDFSA1X2sW8yo0ehOp+Vr6dHS1U1QCgeNhMBj52zufEZd0FbVGi5NKSXhYRz5YuhCZTIwP7Y14B4DgdkpUcsvUilQC/l4KO0UkEDgmkefiiEu6SoVai8kEFWotcZeuEnkuzt6hCRCCDkCvYBdCfVRmUZdJwWiCveeuU1gmFiMJBFVcupxOhVprsU2t1ZJ0JcNOEQmqI1IugFQq4YX7/S1cLnqjia+P5LL8x3SeneBLJ18ne4cpENgFg8FI5Lk4Ll1O59TZWKvHnZRKwkIC7BCZoCZC0P9AKpXQp6OrRc7c11PJ6n1ZfLgjg7kjvRnW1d2OEQoETU/1nHnVyFwhlyGTydBotTgplYR36UhEv3A7RyoAIeg3JbCtklemBfLFoRy+PpJLap6GaUPaEpdeIfzqglZB9Zx5FTKplEcmj0IhlxMWEkBEv3AxIeogCEG/BW5OMp6f6MeWU/kcii7mZGIpeqMJrfCrC1oBtnLmGp0OhVzO/Icm2CkqQW0IQa8DMqmEmcPbgwkOxxSbt1f3qwt7o6ClIHLmzRch6PXAzVlmtU341QUtCZEzb94IQa8HVX716ouQZFIIaqu0Y1QCQcMhcubNGyHo9aDKr15VJkAqAYMRTiaW0C3QGSeF+JALmjeJl9NEzrwZIwS9HtT0qwe1VZJTpGXr6UL+vT2D5+7xxdtDrC4VNC+qcuaJl9M4duai1eM1c+Ymk4kdO3bQv39/goODmzJUwS0Qgl5PrP3qrgS2U/HFoRyW/5jO0+N9CQ90tmuMAkFdsZUzVyrkSKXSWnPmEokEV1dX0tPThaA7GELQG4AeQS68Oi2Q1fuy+M/uTGYMbcfYXh5IJMLKKHBsbOXMpRKJzZx5QUEBarWagIAAxo0bZ8eoBbUhBL2B8PZQ8Pepgaw9ksPmk/lcy1PTt6Mrmdd1YgGSwGFJSEqtU8782rVrfPfdd7i5ufHcc8+JwYqDIgS9AXFSSlk4wZedvxWw62wRkZfKMIFYgCRwKKpy5nFJ1zj861mrx2vmzC9evMj27dvx8vLi0UcfFWLuwAhBb2CkEgkhPk7IZUXoDZXbxAIkgaNQlTOPvVRZzxxApVQgkUiscuYmk4ljx47x888/07FjR2bOnImzs5gfskVKSgpLlizh+vXreHl58d577xESEmK1388//8zHH3+MyWTCZDKxaNEi7rnnngaLQwh6I5CarzWLeRViAZLAEYg8F2ch5gASsJkzNxqNpKen06dPHyZPntyk7fvMq1VT0ukSGujw3vc33niD2bNnM3XqVLZv387rr7/ON998Y7GPyWTi5ZdfZsOGDXTt2pX4+HhmzZrF+PHjkUob5rkJQW8EbC1AAricrUZvMCGXiVtWgX24GJ9iIeZgnTOvqKhAr9fj7u7Oww8/jEwma9I0S1N0RVqxYgWxsbGUlpZSWFhIcHAwq1atwsmp/mWy8/PziY2NZe3atQBMmjSJZcuWUVBQQNu2bS32lUqllJSUAFBSUoKPj0+DiTkIQW8Uai5AUskluKikRKdW8P92ZbJwgi8eNsoICASNQdVo93zsZfYeOW31ePWceWFhId9++y0qlYoFCxY0+Kh879Ez7DlsHUN1ikrKuJqWjdFUOSCqUGs5F5PEgsUf4Ole+x3uxLFDuG/U4DrFERMTg8Fg4KuvvkKlUjFv3jz27NnDtGnTzPs8//zzXL161ebxmzZtMot/ZmYmvr6+yGSV32mZTIaPjw+ZmZkWgi6RSPjoo4/405/+hIuLC2VlZXz++ed1ireuCEFvBGw1zOgV7MKZ5FLW/5LHP39I47l7/OjorbJ3qK2O5nYrf6fcyJlfQa3RAeDspMRkwipnnpqaynfffYfJZGLSpEl2m/ysUGvMYl6F0WSiQq29qaDXh+joaL799ltcXFwA6NKlC9evX7fY5z//+U+DXKsKvV7PZ599xurVqxk4cCBRUVG8+OKL7Nq1C1fXhnleQtAbCVsNMyK6uOPfprJpxvs/ZTBvZHsiuoimGU2F0WhqdQ2OI8/FEZN4BY1WZ95mMpp4ZMpoi5x5fHwcP/74I56ensyePZt27do1Sjz3jRp8y1H0iagY3v5ovYWd0lml5IUF0xg+sOcdx5CWloZWq6Vz587mbRcvXmTSpEkW+9V1hO7v7092djYGgwGZTIbBYCAnJwd/f3+LY+Li4sjJyWHgwIEADBw4EGdnZ5KTk+nTp88dPy8Qgt7kdGiv4tXpQXx+MJuvfs7laq6GrgHOpBeIhhmNzbnYJItFNNUbHDeEUDgK1e9CohNSLMQcrHPmBoOBo0ePEhAQwKOPPmoetdqLiH7hhId1JO7SVdSNUOExOjoajUZDWloaQUFBbNq0CVdXVwYPtvyhqesIvV27doSHh7Nz506mTp3Kzp07CQ8Pt8qf+/n5kZWVxeXLl+nUqRPJycnk5+fToUOHBnleIATdLng4y3jpAX82ncjjUHQxh2OKMZmEX72xSUnNspoQVGsqGxy3FEG3tZS/JlU5c4PBgMlkQi6XM2/ePJydnZvUyVIbMpmUD5YuJPJcHElXMhq8wmNMTAwzZszgpZdeoqysjG7durFy5co7Ouebb77JkiVLWL16NR4eHrz33nvmx55++mmef/55evfuzZtvvskLL7xgTmf985//xMvL646uXR37v3utFJlUQq9gF47Hl6A3Vm5rzX51g8HAoUOHkEqljBkzxjzBtGfPHiZOnFjn8xQXF1NcXGyxLSsrC4DQYD+cVEoLoTMBmTn5rNuyv0Xk1G0t5ZdIJCjkMnR6vXm027d7KBs2bMDFxYWHHnoId3fHSv3JZFKGD+zZKD+00dHRPPPMM7z99tsNds7OnTvz/fff23zsiy++MP9/ypQpTJkypcGuWxMh6HYkNV9rFvMqWqtf/eWXXyYwMBCFQsGXX37Je++9R8eOHdm4cWO9BH3dunWsWrXK5mP9eoQRHhZvvpVXKOTodHp2/+G6cHZq/jn1SynWLeNMJhOjhvalQ6APYSEBdAvxZ926r8nPz2fy5MmtbuVnbGwsvXr1sncYjYIQdDtSm189NV+D0WhqVWmXnJwcPvzwQwCmT5/O3//+dxYtWlTv8zz22GMW1jOoHKHPmTMHqVRicSuv0+vZvOOI2f3RXHPq1VvGRUUnWj3urFIy9q5+DB/Yk/T0dNau/QqDwcDcuXMJDQ21Q8T2JTIy0t4hNBpC0O2ILb+6s1LK2ZRyVu7N4qmxPrg6tQ6/uk6nQ6vVolQqCQwMZM2aNfz1r38lMdFaoG6Gh4cHHh4etT5e/VZ+3Zb9VhOGFRoth4+faza2Rls5c7lMhlxu3TJOr9ezadMmFAoFjz/+OO3bt7dz9IKGRgi6HanNr348oYSNx/P454/p/OlePwJbQYu7V199leLiYrPIuLm58emnn7Jz585Gu2aX0ECrnDrA4eNnMZqMzcLWaGspv1xm2TJuSN/uSKUSJBI5M2fOxMvLq8F8zwLHQgi6nbHlVx8R7kFAWyWfHcjmvW3pPDbam4Gd3OwYZeNjy4crk8mYOnVqo12zpj1OLpOj0+sxGCsnNhw1BVPdlph8JeOmS/mNRiO7d+/G09OTESNGEBgYaKeoBU2BEHQHpbOvE69OC+Szg9l8fjCHe/tq6OzrRFoL9atXVFSQmJhIQkIC8fHxJCQksGHDhka9Zk173NX0bA79+jvVFymqtY5la7Soc6LWYrKxT5UtUaPR8P3335OcnMzdd9/d5LEKmh4h6H9gNJo4l5THlaxiQvw86BfW3u6C6eUq56+TAvju11z2nS9CIilqMX715ORk9u7dS0JCAgkJCaSmpuLm5ka3bt0IDw9nxowZTRJH9Zz6iagYjp+JtrQ1mqCktNxhbI11tSV2Dw3gq6++Ii8vj8mTJzNgwAC7xSxoOoSgUynmyzdEkZRehFZnQKmQERboyStzBtpdMBUyCX1DXDl1qbRF+dUXLFhASEgIAwYMICEhgcGDB7Ny5cqbTmg2JHFxcYSGhlpU17NOwcjQ6Q18v/MoAE4OYGusiy1xQM8wPv10NRqNhjlz5tCpUye7xCpoeoSgA+eS8khKL0KjqyxirtEZSEov4lxSHgO6ets5upbpV3dxcWHt2rVIJBIWLlzIihUrmDp1Km+99RYjR45s1GuXlJSwdetWPD09eeSRR/D19QWsUzA6vZ5NPx0xO2HslVOvypknXk7jzLkEq8er2xKrGDduHH5+fvj4+DRZnAL745hT903MlaxitDrLjhRanYGr2SV2isiSKr96TVLzNRiMtrKojs/u3bvNC1pUKhWvvPIKK1asYPny5bz22muUlpY22rXd3d2ZP38+Wq2WL7/8kgsXLpgfq0rBzH9oAnKZDK3Otq1x3Zb9nIiKwWAw1jx9g1KVM3/r/61n7eZ9RCdeQSGX4aRSIpFUinl4l44M6dudU6dOER8fD1ROMgsxt01kZCTz5s1rsut169atya7lECP0jIwMFi9eTH5+PqGhoXzwwQdWtqqMjAweeOABcyGb9u3b89///rdBrh/i54FSITOP0AEkEvBt6xjttm7qV9+TxVPjfHBrAX71/v37s337dj766CMefPBBDh482GjX6tChAwsXLmTLli38+OOPpKenW61Irc3WeOj4WUxNZGu0VS1RJrW0JQ7u0439+/dx5swZ+vTpQ/fu3RslFoHj4xCC/tZbbzF79mweeOABPvnkE1avXs3ixYst9rl48SKTJ09u0PoLVfQLa09YoKc5hy6TSdAbTGw8mEh7T2e6Bnk1+DXrQ21+9ZOJJXz7a6Vf/dkJvnRo3/zrqyuVSl5++WXuu+++Rr+Wm5sb8+fP59ChQzZ92bXZGo2NbGusbkv87ULiTaslVjpZNnPp0iWGDRvGhAkTGiyOxuLrr7+22tazZ08GDx6MTqez6W7q168f/fr1o7y8nM2bN1s9PmjQoHot5y8sLGTBggXk5OTQp08f3njjDY4fP85HH32E0WgkODiYt99+m/bt2zN27Fi++eYbgoKCiIyMZNWqVaxfv5558+bRu3dvoqKiKCgoYOnSpYwaNYq0tDQWL15MeXk5ffv2rddrc6fYXdB1Oh1nzpzhk08+ASqXfc+dO9emoCcmJjJ9+nTc3Nx47bXXbN7K3Kw4U21IpRJemTOQc0l5XM0uoaOvO65OclZvv8jbX59hxujOTBkeatcJUlt+9bu63/Cr/3t7BnNHtmdoC6mv3lD1oW+FVCq1EMGkpCSkUimdOnWqk62xoVeW1qdaokajYe3ateTk5PDAAw8waNCg275uayMtLY1Vq1bRsWNHXnrpJT7//HM2bdrExo0bCQoK4ssvv+Ttt9++ZQldnU7Hpk2bOHz4MB9//DGjRo1i2bJlTJ8+nYcffpht27axadOmJnpWDiDohYWFuLm5mct2ent7k52dbbWfSqXiwQcf5NFHH+Xo0aP8+c9/Zvfu3SiVlqsob1ac6WZIpRIGdPW2mARd/vQwvtwVx+afk4hJKeBPD/aijXv9ew42JqE+Trw6PYgvDmaz9o/66tOHtCUuvcJiNG9vt05zwGQy8csvv5CWlsaYMWO4++67b2lrBPj5xDkMRkODpGAqV37e6C4Etm2JEf3CkUoldOrUifHjxxMWFnZHz70pefzxx2t9rKosQW24uLjc9PG6MmjQIEJCQgCYPHkyS5YsYciQIQQFBQEwc+bMOrWHGzFiBGDZ8ej06dPmukRTpkxh6dKldxxvXWlSQd+zZw/Lly+32Fb1olbHVvW3v/zlL+b/jxo1ig8//JDLly9b5QtvVpypvrg4KfjL9N706dyOdXvjWPL5SZ6Z1BOJROJQfnUPZxkvPuDP1lP5HIou5mRiCQYjaP/Itzd3z3pTIZFImDt3Ljt27ODw4cOkp6fz4IMPmq2NNVMwMpkMvd6A3lA593K7KZjqLpYTUbEWYg7WtsS2rnLy8/Pw8fHhnnvuabgXoBVRve67yWSy0hyTyYRer7f4G7DYBpUDTbDWrKr9JRJJk1azbFJBnzhxotXEk06nIyIiwty+KTc31+bs/Pr165k0aRJt2rQBMBfmr8mtijPVF4lEwuh+gXQJ9OQ/W8/z4eZzyKQSjEaTQ/nVZVIJjwxvj0QCBy/eSDm1BM96U6JUKpk+fTpBQUHs37+fzz//nKeeegoXF5c6p2ASL6cD1CkNU5Viibl0BU0NIa+iui0xMjKSzZv30a1bN2bOnNngz7+1EBUVRUZGBn5+fmzbto2nnnqKb775xqKLUUREBABt2rQhKSmJ4OBgDh06dMtzDx8+nJ9++ok5c+awf/9+tFrbqbPGwO4pF4VCwaBBg9i9ezeTJ09m27ZtNn3IZ86cQa1W8/TTT3P69GmMRmOTLpgI9HbjoVFhrPzhgtkq6Gh+dQAXlbXbpbl71psaiURCREQE/v7+xMfH4+x8w+1UlxTMll1H0ekNldUOVUq6d+7AjPtHknwtgy6hgQzq043fLiRwKSWd/MJizsclW9gfpRIJ8hoplsF9urFnzx5Onz5Nt27drO5CBfUjLCyMV199ldzcXIYOHcqCBQsICwtj0aJF6HQ6AgICePfdd4HK3qLLli1j1apVdSqh8Prrr7N48WK+++47evfu3aSF0CQmk8nuRub09HSWLFlCfn4+/v7+rFixAk9PTzZu3EhOTg4vvPAC2dnZLFmyhNzcXFQqFe+++26d7VlpaWmMGzeOQ4cOmXNkt8MPvySz9WiyVf2MGaM6M31kZ5vHNDUXrpbx5aEcixrrUgk8Pc6HAS28wFdtNMT7n5eXx2+//cb48ePNd4bmCcw/UjAqpRInlYLrxWUWx1bPgSsVij9WoOrR6vS2LgXAhBEDzSmW/j06s23bjyQmJjJ06FAmTJiAVCqWkAissfsIHSAwMJD169dbbZ81a5b5/76+vqxdu7Ypw7LCll8dIDoln/uGdMDFSWGnyG5Q07Muk4LBCHvOXSfEx4m2bg7xljc7kpKSiIyMJC0tjYcffhhPT0+bvS8Tk9P4+vt9Fj/6JpPJLN4arQ4NtlMrVdRc+anX69HpdEycOJEhQ4Y01lMUtADEt7se1PSrKxUyPN2UJFy7zitfnGTRtD50cUDPut5oYt2RXN79IY2nxvkSHugYC6aaE0OHDsXT05Nt27bx2WefMWPGDLO1sWbvSycn68VIt0Iul2EwGCxcLNnZ2bi7u+Pi4sK8efNaXas4Qf0Rgl4PbPnV+4W1JzmjiFU/XuCtr88wY1RnptzleJ71wDZKPj2Qzce7M3lwcFvu7espBKKehIeH4+3tzebNm1m/fj2zZ8+mS5cuFvvUdMIoZHL0BgPGm2Q2nZQKHpky2rzyM6JfOCkpl/n+++/p0qULM2bMEO+VoE44RA69sWmoHPrNKFfr+O/uOE7GZNGjYxuendKLazmlDmVvVOuMfHM0l6jLZQwIdWXeiPYkZatbvF+9od9/rVbL8ePHGTFihE2nVZUNMelKBp06+LFl9zHiL12rzLMrFJXNqasmTf8YkVf3rp85c4Y9e/bg6+vLrFmzmqwCpaD5IwS9ATGZTPxyIYO1u+Mw/NHkWa83OpS90WQycfBiEVtPFaCQSUDS8v3qNd9/nU7HzJkzmTZtGu+88w6vvPKK1WKVhIQE5syZw/bt22/a5UetVrNlyxbGjx+Pn5+fzX2qC3xYSIDZ5VL1d5Wt0Wg0cuDAAU6dOmUemddcOCcQ3AwxVd6ASCQSRvUNZN493TAaTej0RkxY2hvtjUQiYUIfL6YOboPWYEL7hxumul+9pbNhwwauXbtmLvRmi27dutGnTx+zda02ioqKyMnJ4b///S/nzp2zuU/1Co7DB/ZEqZBb/F01Mler1cTHxzNkyBAeffRRIeaCeiMEvREoKrOeENM4UDleAFtVd6v86i0Zk8nE119/zYQJE8wrQGNjY5kyZQqDBw/m1VdfRa1WAzBp0iQOHTpESkpKrefz9fVl4cKFBAUFsX37dnbu3Gm1mvBWlJWVYTAYcHFx4ZlnnmHixInClii4LcSnphGosjfWpKBYjaNkuGzVWJdIoF0LtzXGx8eTmZlp0ZLt5MmTLFy4kGnTprF161a+/PJLAAYOHAjAkSNHbnpOV1dX5s2bx1133UVUVBT79u2rczw5OTl88cUX7N+/H8BiEZNAUF9a9rfXTljbG6Uo5DIO/Z5GYYmGZ6b0xMPFvrfTNf3q8j/86j9E5uPtoaCzn2MVIasrt6q2ee3aNQBzlyKoLMT0wAMPcP/99/PDDz8QGRnJokWLzPtUHXMzpFIp48ePJzg4mICAAMB2jZDqJCcn8/3336NQKJq8zKqgZSIEvRGwZW/s27kdB6JS+fZgIks+O8lzU3vRu1M7u8ZY06/u6Srji4M5fLAjg4eHtWNMT49mZ5era7XNqprmUFlPqAqJRIJCoTD/H6hX+qOqpLPRaOTbb7+lQ4cOjBgxwuJ1NJlM/Pbbb+zZswdvb29mz56Np6dnna8hENSGEPRGwlY53vuGdCS8Q1tW/XiB5RuimDQshBmjOhOdUmAXe6Mtv/qr0wJZ+3Mum07kczlbzdyR3jgpmk9m7lbVNkNDQ4HKVEfHjh0B2LRpE6GhoSQkJFBcXMzw4cMBzGWcb6dmkNFoxMXFhZ9//pm0tDSmTZuGUqlEJpOh1Wo5cOAAnTt3ZsaMGeaKfQLBnSIEvYnp6OfOO08N5X8HEth58goHo1IxmkzodI5hb3RRyXjuXl/2nrvOT78Vkl6g5ZlxPuSW6JuFX/1W1Ta7du1KYGAgUVFRZkGfO3cuX375JVlZWcycOZP58+cDcO7cOaRS6W01rZbL5UybNo2goCD27dvHypUr8fLy4umnn0alUjF//nwCAgLE5KegQRGCbgdUChkL7u+Bm5OC7cdvOCgcpXqjVCLh/v5tCPVW8cWhbN7amo5MCnoDLcKv/sQTT/DJJ5+wbNkyEhISAFi0aJHVfvv27TPnxW8HiUTCkCFDCAgI4OjRowQEBGA0GpHJZI26HkLQehHDAzuikFu//I5kbwwPcuGhoe3BVCnm0DL86rNmzSIgIIDt27fXus+VK1eIjIxkyZIld3y9oKAg5syZw5gxY5DJmn8zb4HjIkbodiTEzwOVjeqNjjTwvV6mtyoX3Nzqqxv+6ChU3e1S1SsyLS3N5jFyuZxt27ZhMplq3UdgjZ+fn81yCIKmQbzydqSmvbFqxL7p5yQqNAZmjO6M/A4aDjcEVX716vXVARzoN+eW5ObmAtxWG0JB/Wjs8hqCmyME3Y7Ysjd27+DFhoOJ/HQiheiUfP48rTf+7ew3Eq7pV1fKKqV8+2+FIIH7+nkhdXBrY69evdiwYQPe3t7IZDKz62XDhg211l8R3JzaXkPxetoXIeh2xpa98elJPenbuT1f7Izh1S9O8dh93RnVN8AunnBbfvUwPye+/TWP7WcKSc5S8+QYH1ydHDc37OTkxKBBg6y2+/n5idHkHSJeQ8dCCLqDMiTcl86Bnny6/SKf74jhbGIOET38yCoodwi/+oKxPoT5FbP5ZD7v/pDOM+N9CPFpnqtLBYKWghB0B6adhxOvzhnEjpMpbDqcxJmEylywygH86hKJhNE9PQnxVvHZwRze/ymDh4a2pZ2bgrQCx/erCwQtkXoJusFg4NChQ0ilUgsL1p49e5g4cWKjBNjakUolBHu7o5RL0eorl6s7il8dIMTHidemB/LV4Ww2nShAKqms5NgS/OoCQXOjXhaKl19+mejoaOLi4pgzZw5Xr14FYOPGjY0SnKCSK1nF6PRGi20anYGLl+1fXx3AzUnGqB6eyKQ3yvI6sl/dw8ODRYsWiU5Ad4B4DR2Teo3Qc3Jy+PDDDwGYPn06f//7322usBM0LFXleGv61Q9GpeHfzpUJg4LtXkQrrUCLwfI3p9KvnqdxOL+6h4cHf/nLX+wdRrNGvIaOSb1G6DqdDq22sgFCYGAga9as4auvviIxMbFRghNUUuVXVylkSKjMoXcL9qJnaFu+3hvPv787y/VSjV1jtFVfHSAhQ41GZ7RxhEAgaGjq1VP0woULBAQE0L59e/M2g8HAzp07mTp1aqME2BA0VU/RxsRoNFn41fuFtUcigQO/pbLhYCJOShnPTOrJwG4+dovv492ZZr+6Si7BzVlGfokefy8Fz4z3JaCtaKkmEDQmokl0CyAtt5RPfrzI1ewSxg4IYva4LsRfu97kJXmNRpOFX71XsAsJmWq+OpxDhdbI7LvbM7ybe6PHcSt27NjBp59+ik6n4/HHHxcrSG+D0tJSHn30UdasWdMiv1PNlXoJekVFBYmJiSQkJBAfH09CQgIbNmxozPgahJYu6AB6g5Hvf05ix8krKP8oIaDTO0ZJ3qJyPf89nENChpphXd2YOawdl7LUdinHm52dzaxZs/jhhx9QKpU8+uijrFixgrCwsCa5fkvg/PnzLF26lJSUFPbu3dtiv1PNkZtOiiYnJ7N3714SEhJISEggNTUVNzc3unXrRnh4ODNmzGiqOAW3QC6TMmt8V1yc5Gz6Ocm83REsjp4ucl6835+dvxey6/frRF0uAxNoDaYmtzeeOHGCoUOH4uXlBcC9997L3r17xeR+Pdi8eTNvvPEGL7/8sr1DEdTgpoK+YMECQkJCGDBgAAkJCQwePJiVK1cKq5IDYzCakIBFhcSqkrx2rbEulTBlUFtkEvgp6vqN2KrZG5vCDZOTk4O3943XwcfHhwsXLjT6dVsS7777rr1DENTCTV0uLi4urF27lueff56ffvqJ7t27M3XqVH755Zemik9QT6osjjW5XqLGIaZLbNgrq8rxNgW2XgN7Wz4FgobipoK+e/du84ddpVLxyiuvsGLFCpYvX85rr71GaWlpkwQpqDvWFkcprk5yDkSl8b4D2xttbWsMfH19ycu7sSArJycHHx/7OIMEgoam3sW2+/fvz/bt2/H09OTBBx9shJAEd0JVSd6/TO/DjNFh/GV6Xz7962geu7c7MVcKeHnNCU7HZdstvqpyvFUCrpBJkEtha2QBBy5cx9jIdxHDhw/n5MmTFBQUUFFRwf79+2+rZ6hA4IjckW3xwoUL9OnTpyHjaRRag8ulLqTnlbF620VSMosZ2SeAuRO6kphWZHd7Y6i3iv/9mse5K+X0Cnbm8dE+uDs3XjneHTt28Nlnn6HT6ZgxYwZPP/10o12rJTN27Fi++eabVv2dcjSED72VoTcY+fHYZX48dhmFTAISCXoHsDeaTCaOxBaz5VQBriopT47xoXugc5PHIRA0Z0T53FaGXCbl4dFhOKlkbDx4iSo/jL3tjRKJhDE9PQnzc+LLQzl8tCuTe/t6EurrRLooxytopqSkpLBkyRKuX7+Ol5cX7733HiEhIRb7FBYW8vLLL3Pt2jWUSiUdO3bk7bffpm3btvW+nn0bVgrshk5ntOoLWmVvtCfB7VS8Oi2QYV3d2Hu+iDX7s/npt0K+PJTDx7szMRpb/A2l4CYYjSYuXC1j1++FXLha5vCfhzfeeIPZs2ezb98+Zs+ezeuvv261j0Qi4amnnmLfvn3s2LGD4OBgPvjgg9u6nhiht1Jqq+CYllOC3mC0a3NqlUJK/1BXTieXov8jvKb2qwscD1v1ghp6UdqKFSuIjY2ltLSUwsJCgoODWbVqFU5O9e/GlZ+fT2xsLGvXrgVg0qRJLFu2jIKCAovRt5eXFxEREea/+/Xrd9slyYWgt1Kq7I1J6UVodQaUCilOSjknY7PJLCjnuam9CPaxX92V1HytWcyr0OhNXM5WC0FvgZxMLOFEws3vDkvVBjKv6zBVq7mfkKlm2dY03G7S03Z4N3eGda3bZzkmJgaDwcBXX32FSqVi3rx57Nmzh2nTppn3ef755829IGqyadMms/hnZmbi6+trbgQkk8nw8fEhMzOz1nSK0Whk48aNjB07tk7x1kQIeiulyt5Ys4JjVGIO/90Vy2tfnuLh0WE8MDTELnnrKr+6Rm95S30svoQeQS50DRATpq0Njc5ETQuHyVS53a2B2tlGR0fz7bff4uLiAkCXLl24fv26xT7/+c9/GuZiNli2bBkuLi7MnTv3to4Xgt6KkUolDOjqbTEJOri7L92C2/Df3bFsPHSJ3xJyeGZST7ILK5rU3ljlV69+e+3npaBcY2DFzkwm9PVkyqC2lU4dQbNnWNdbj6IvXC3jy0M5Fj/yKrmER+9q1yB3bWlpaWi1Wjp37mzedvHiRSZNmmSxX11H6P7+/mRnZ2MwGJDJZBgMBnJycvD397d57HvvvcfVq1dZs2YNUuntpTyFoAus8HBV8uKMvhyPzmLtnlheXnMCmUyCwWBqMnujVCrhhfv9rcrxag0mtpzMZ//5IuLSKnhyrA8BbUSd9daArR/5UB8VvYJdGuT80dHRaDQa0tLSCAoKYtOmTbi6ujJ48GCL/eo6Qm/Xrh3h4eHmfhE7d+4kPDzcZrplxYoVREdH8/nnn6NU3v7nWfjQBTfl6Ll0vtgZQ3UzgUoh4y/T+9i12Nf5K2V880suGp2JB4e0wdtdQVoLsDfq9Xo+++wztm7dSllZGUuXLiUrKwu9Xs9zzz1n7/Dsjq2a+w31Xn/44YcUFRURFxdHWVkZ3bp144033jBX5rwdkpOTWbJkCcXFxXh4ePDee+/RqVMnAJ5++mmef/55nJycmDRpEiEhIebRfVBQEJ988km9rycEXXBTfvglma1Hk6n5IZkxqjPTR3a2eUxTUVSuZ92RHGLS1EgklfnUpi7H29B88MEHREdHs3LlSs6cOcP777+PVCpl06ZNuLm52Tu8Fs0TTzzBM888w7Bhw+wdym3jECmXjIwMFi9eTH5+PqGhoXzwwQe4ulrmxJ599lkyMzOBypngxMREtmzZQu/eve0RcquhNnvjqdgshvX0w7+d/Rwnni5yRvXwJCFDjf6PtqXN2d5YWlrKunXr2L17N+7u7vTt25fLly/z0ksvCTFvAmJjY+nVq5e9w7gjHGJh0VtvvcXs2bPZu3cvvXr1YvXq1Vb7rFmzhu3bt7N9+3bGjx/PI488IsS8CbDVoDqwvSv5xWqWfH6SXaeu2HVxR1qB1izmVWj0JpKy1PYJ6A44deoUISEhBAcHA5VN2d3d3W/b8SCoH5GRkbi7279F4p1g9xG6TqfjzJkz5nzR9OnTmTt3LosXL7a5f3JyMtu2bWPHjh02Hy8uLqa4uNhiW1ZW1i3jMBiNRMVc43JaHp2C2jOwZwdktznT3JKozd5YVKbhv7vi2HAgkcjYbBZO7ol/O1fOJeU1qRumNnvjkZhi/NsoGdrFrdnUO8/OzrYo5btp0yZ8fX3F6FxQZ+wu6IWFhbi5uSGXV4bi7e1Ndnbt5V0//fRTFixYUOuHfN26daxatapeMRiMRt5YtYOEqzloNDpUKgXdOvrw1qLJQtSxbW9s4+7E/83sx4noLNbti+eVz0/SztOJ66XaPxYqNY0bxpbzIbCtEjDx9ZFcfksuZe4Ib9q42f2jfkv8/PyIj48nJyeHzMxMtm/fTnl5OVqt9o6cDw1NuVrL22v2kHwllc4hwbz+7ERcnBwnvtZMk06K7tmzh+XLl1tsCwkJ4cqVK+YuSHq9nv79+3Px4kWr44uKirjvvvs4cuQIKpXK5jVqG6HPmTOn1knR0xev8P7XB1BrdOZtTio5ix+/hyG9Q+r7NFsdRaUaVnx/jktpRRbbm8oNY8v5AHAktpgfTxcglcDDw9pxVzd3hx6ta7Vali5dyuHDh/H09GTlypW8//77qNXq214K3tCUq7VMefoddOoyKgu7SVA4ufLTF0uFqDsATTpsmThxIhMnTrTYptPpiIiIMJvvc3Nza+0gc/ToUUaOHFmrmAN4eHjUu+fp5bQ8NNXEHECt0XPgZBx9uwWhUjr+6M6eeLqp6Nu5vZWgN1UvU6lUQp+OrlaToGN7edK7gwvfHM1l/S95/JZcxuy725F1Xdcotrc7RalU8u9//9tiW1UdEEfh7TV7qok5gAmduoy31+zhXy9OtWdoAhwg5aJQKBg0aBC7d+9m8uTJbNu2rdYOMufOnWPQoEENHkOnoPaoVAqLEbpEAqcupPDEP77hvrt6MHFET7zbuItcey2E+HmgsuGGib1SwPBefvi2aZjFH/XF20PBS5P8ORZXwpZTeby+OQ2ZBPTG5m9xbEpSs6+z/Is9xEWfAysTq4nkK6l2iEpQE7sLOlSWmFyyZAmffvop/v7+rFixAoCNGzeSk5PDCy+8AEBqaiqjR49u8OsP7NmBbh19LHLoXTv68Mg9A9l17CJbD5xl68GzDO0TSnZ+Cek5hWg0epFrr4atYl9uzgoSUgv52+rjjO4XyIMjOtHOo4GKbtQDqUTCqB4eyKWw/lhei7A4NgY1c+OLHx/HF9t+45cTpym/novJZACJrc+5hM4hwU0er8AasbDoD6pG3inpeYQGWo68s/OL2fVLNLuPRaPR6i2Oc1IpWPz4BJFrpzKXbcsNs+3XFA7/noZUImHcwCAmDQvhSlZJk7e+2/V7IT/9Vmi1fXxvDx4e1r7Rr+/IWOfGASTm/zu7t2XU8ME8+8jdPPzn90QO3UERgl4P/rczkk17o6y2Tx7dm2dmjLiTEFs8udcr+PHYZY6eSwcq894GowlVE7a+s1XcCUAqgQl9PJnYvw3OytZ3p6U3GPnzu5uJv3iGmukUj7Z+/OcfTxIadOMHT7hcHBeHSLk0F7p29MWpRq4dYMeRi1zNKGDi3T0Z2jcUuUwmcu018PZy5pnJPekc4MHaPfEYjDda3yWkFnLkXDpjBzRuWQZbFsegdkrau8vZd76IE4mlTB3Uhru6VS4uaayaIU2JXm/kp8h0EtOK6RrkwZSIQORyKVq9gQNRGWw79Dupl+OpKM7FOjcOTs4uFmIO4OKkFBOgDooQ9HpgK9feKag9A8KD2X8ijve+2k8bDxfGDe1OTFIGKen5wtdeg6IyLTVvCvUGE1/uimXv6asM6OrDgK7ehAV4AjToQqXaKjhKpRLG9lKz+WQ+/zuWx8/RRUilEnKKdI3WGacp0OuNPPfxcXJSE9CWF3HSxZMff+mEi1RL2pUEKopzMGgrAJCrnNFrKmqcwfFy41VpvaZO190udekpWpNVq1axcuVKduzYQdeuXet1PSHo9UAmlfLWosk2c+0z7hnA2bhUdh+LZsv+3y2OU2t0JFzNISrmWqvPtduqDaNUSBne05+c6+XsPHGFn46n4O6iQCaVUFahR28wNthCpdosjiE+TiyeEkDU5TK+O55HifpGPQFHnDytS9XB73+5yuWz+9FrygET5YUZFGUkVE5uAkqVE6NHDuWpR8bh5elu01/++rMTra5tL4xGE8s3RFWbeG+6dN3tUtVTdOrUqWzfvp3XX3+db775ptb9Y2JiOHfuHIGBgbd1PZFDbwT++8Nxth0+b7V9yug+PD3j7ka/viNzqy9laYWOC8l5HIxKJf7adYtjZVIJM0Z35v6IEBRyqcU5G3LU9tNvBez6/brV9ru6uTF3hLfdxcNoNPHRrgwSrxWg1mhwUqno2qEtj4/2ISVXw4WUIi4m53EpOpLizARqplJc3L14+4VHGdSni8VCK3vmxn85n8HR8+k33aekXEt6XplF1yKJBALbu+LuUnuco/oGMrJvQJ3iaOieovfeey+RkZHmBhcRERHs37/fZk10rVbLvHnz+PDDD5k/fz5r1qwRI3RHoHeXQPYej7XKtf905AKXruUwfmh37h4QhouTstXl2murDVMlkm7OCob38ieroJyEa9ctpMhgNLHpcBI//HKZrsFe9OzYlu4d2rDlaBLJGcUNNmoL8VbZrA9zPKGUi9cqGBDqysBOroT5VX7J65Jrr2sd77rsd/5KGWdjktGU5mPQaShXqIgq8uJcfCra0nw0pQVUFGehLrGdF/fx9mFwX2uhcPTcuFprsN2CTmvAvYGWOdizp+jHH3/MlClT7mjQKQS9EbCZaw9sx+DeIRw6Fc/Kb4/w+ZZfuatfZ65m5pOec71V+dpt1Yapia3UjEohZWJERzQ6AzEpBWw+kmR1nEZnICm9iHNJeTbPX5fRvK3J0w7eKkZ0d+dsShnHE0o4EluMu5O0Mi2kMaIz1J5rr2u3+pr7KeUS/NsoGNvLk9wiHdl//EtJL6A4KwlNWT5GnQapQoVc5YZBp0ZXUYROXQomIyqVCo3B0mYLEkYNrt+orykY2TfglqPo3xNzWfnDhRqfCRmP3RfeYKuR7dVT9OzZs0RHR/O3v/3tjs4jBL0RuFmu/aHx/UlIyebgqXh+PpOAttqHU+Tab2C9UKly5D1jVJhZBIvLtazbE8/JWMtqmhqdgc93xjCwqzedAzzpHOBJkI8rEiR1ysHebPI0oos7ap2R6GvlHI4uIjlbc+O6ehPxGWr+sSkVbw8Fbk5S3JxklKkNJGVZ1mxPylLz9ZEcPF3lqLVG1DoT2UU6ruaq0WvKMeg06BUqruhcWPtzLiajAReZBolBTUn2NUpykjHo1NT0jCsUcsbd1Y9J4yLo0bUjj/zpXxQVF1UOZSUSPD08mfug7ZXYjk5tn4l+YQ2zhsCePUXPnDlDcnIy48aNAyrrTy1YsIDly5dz9911T9MKQW8kZFIpQ3qHWAmzRCKheyc/unfyw8vD2crXrtbo2H0sms7B3rTzcowJOHtwq9QMgIeLkrt6+/P7pVyLUZtMKqGdu4oz8Tn8fLYyL6tSSPH2ciYzv9zCMnmz0bxeU46mtBi9mwdw457eSSFlUGc3sot0FoJe/fpqnZH8Eh0laiMVWqPVPnojRCaVoZBJcFJKcVZUHlOWn462rDKVIlOokDu54+GipLi4lAJNGUZdOeWFWRh0NR0pJjp16sQnbz6Fi/ONfO+WNa/y7fZjxCRepWfXjsyeOgKlvHl+7evymbgT7NlT9JlnnuGZZ54x/z127FiRQ29u2PK1SyQQFXuNJ/6xjl5hAYwc2IXh/Trj6qJsVbl2qFtqprZR2ytzBiKRQHZhBUnpRSRnFHEmLtss5lVodAa+O5zI5cwigtq7Eejtio+XCx9sOnvLkXxwOyVKGZSXl5kF2MXFlRlD21q4Yc6mlPHfw9lUVNvP2cWVBWN96RfiwvVSLddySth16hqXLEbeEiTSXNRyOVpNBQZdBTKjGnWJ9WpXgC4dfCzEHEApl/P4Q2Pq/+I7KHX5TNwuMTExzJgxg5deesncU3TlypV3dM4333yTJUuWsHr1anNP0Sqqeoo2ZKMe4XKxI7XVYX/m4bv59fdkfolKIj3nOlKJBFcXJWqNHr3e0Gpy7XXFVskBW6M2WzlYqQQ8XFUUlWoskhdgOZ0ol0mYMjyUAd188HJV4uGqBBP8eeVJSsrKzSkNd1cXPvnLMOTVXDhanYG/rDxJSVmF+awKhZywAHfSckspKa/8QddVlFBRlIVeW27OjUukcmQmDSZtKaVl5bg6OxHetQO/X0zCaLwx8pdKpSz72+PcPbh5t1CzJy2hp6gQdDtzsxoyJpOJy2l5bNr7GyfPp1gcp5BLWTRrNGMjutsj7GbJzSyTeoORzPxy0vNKORiVRvw126Pg6jgpZWi0BgvhlwBebiokUtDqjGh0BnR6IyaTCYO23DxClyld8G/rgn8bBdqKMrJy80lJy6W8MKNGbrySof3DuXfUIO4a1BO5XM7/vbOGmISraHU6lAoFPbt15MOlzyKTiR/42yUiIoKDBw826zZ0QtCbAd/t+Y1vd522YUCDbiG+DOkdQkTvEDr4t0UikbQ6K2R9qMto3tZIXqmQMmNUZ/zaulJUpqWoVMPZS7kkZxTXvATBPm50DvBEqZCikstISr9O1LkYS6GWyFApJOZib6GB7SguzCElJcXqfDOnjudPc++32GYwGIk8F0fSlQzCQgKI6BcuxFwgcujNAVv12lUKOcP6hZKeU8T6HZGs3xGJbzsPhvTuSExSJhk519FoW48Vsq7cSV7+/ogQC/EP8fPgP1vPU15WciOH7urOzDFdzOcvKVPz6fdp6LXlFqkUudIFv/beDOzmTUZGJid+O09pue3G1i4q66+pTCZl+MCeDB/Y8w5fEUFLQgj6H1SNeC6lpNMlNNChRjy2fO3dOvrw4rxxyKRS8q+XcSb6CpEXr7D7lxgM1XKrao2O+JRsomKuMqR3qMV5xUjeNnV1U/Tu1BZdSSYVJaXmHLrcUEpFWShfbE0g+lIGKel5GI2mGqkUCRKpjJSyTC6eLcdZpWRERB+C/Nvz7bbDqDVa8zWcVEq6drq9ZeCC1odIuVAp5n975zPikq6i1mhxUikJD+vIB0sXOoyo3yzXXp31OyPZbKPEr7NKQUSfEPp2DaJPtyDaebmKxth3yOmLV/j32v1WNfIBFHIZ4Z386N0lkNS0DHYeOErNZY4Bft48+cg93D24F85Oqhufw0tXUWu1OCmVhHdxrM+hwLERI3Qg8lwccUlXqVBXjowq1FriLl0l8lycw9zS1uZrr0k3G1ZIhVxKp6D2nItP48iZSwC08XCmuFRttvGJRU2W2Lp7MZlMXMkoIOFyFvFXsvgt5ppNMR8/rDvPPTISnU7PsdMX2X7+opWYA9w7ciATRgw0/y2TSflg6UKRGxfcNkLQgUsp6Ra3uQAVGi1RFy85jKDXldrSM28tmoxUIuFaZgHnE9LYfSyGwmLLxSlqjY4dRy/gpFIQ1sHbXJipJaVm6vJcquyk8Vey0Wj1yGRSnJRytDoDOn3lRGkbDxeCfL24dDUHdUWpOTfu6uaBu8rEuys3cOK3GLQ6PW083ZDLZOgNNyZZa0uliNy44E4QKRfgRFQMb3+03jxCr073zsFMGDGQsXf1p62Xu0Pn2quoS3rm9MUrvL92P+pqI8wbDccqFzgF+7Wla0dv4lKyyS0oRae7+SSrowt/bb7/pc/eT2ZuESnp+aSk5XE+MY2UtHyLYyUSGNyzI6MGdaV7Jz+827ihMxiY8ew/byytr4aXhxtjhvdj/N0D6B4WzOJ3PhepFEGjIwQdbOYuO4cEctegnhw6fpakK+lIJRIG9O5CbkER2bmFaLSOmWuvK7WJ298en0ByWh6JV7JJvJpNTFImFTWqRkqlEu7q14mIPp0I9mtDoI8Xcrm0zjn5ugp/Q+936nwKH6w7YJEmkUgqyzEY/0g9KeQy3FxUFBaXWxwrAeZMGsLM+wYBlWsEvt/1C2v+twOD4cYktEwqZd5D45k3fQJyuexGjMJmKGgChKD/wc2+cFdSszhw7Hd2HjrF9eJSi+OcVEreeGles7xFrstIfuOeM2zcdcbKAy+R3BiUSiTg5e5MUYkaY7WPk1Ih408zRzF6cFfza1nbD0lN4b/d/ZRKOR382vDg2H7kFJaQlVdMdn4x2XnFZOWXWHVLAujRyZ+JI3oSGtSeQB9Pfo9N5d9r91NaUmROpbi5e/LyE/fg7eXE4eNnOXT8LOlZeVbnkkjgyZkTmf/QhPq+HQLBHSMEvR6s+34/X23ea7W9S2ggC+dMon+vMOR/1D5uDqmZunD64hXe//qAxSSrk0rBX+eNI8DHk2tZhaRlFfLr2SSuZdpeXSmRgKebM209XZFIJFxJz7OoqaKQS7n3rh6EBfsgk0mRy6RcTstj2+Hz5pw1gFwuZeSAMLw8XCgr11JaoSEj5zpXMvJtzTkC4OHmhF87D/zae2Awmjh9MQWd/saI2kmlYPHjEywmgrV6vWUqRSJBpVQR4ONFSmoWUomE/r260DHIh92HT1vMvzirlLzeTH/gBQ1PfVrQaTQa/vnPf3Ly5ElUKhX9+vVj2bJl9bqemBStB106BeLspLTItcukUq6lZ/O3dz6jjacbo4f1Y+zwvqzdvI+4pGsOa4OsK7VNsg7pE4JMKqVjQDsAOgd7Wwm/Qi5jwrDueLg5U1hUTn5RGcmpuVYFsnR6IzuPRt8yFr3eyOHTiSgVMlydVbi5qKhQa22K+QMjejF/ylBcnG90sqlt1D+wZweLY387n4BGXXbjFsRkQqNRYzSaeOHJaYwe1s88n3IlNdsqNx7RL7yuL6+gnjj6PE1N6tOC7v3330elUrFv3z4kEgl5edZ3gLdCjNDrQW0+4XdffoLfLiRy8NeznIyqdDbUpDmP3OqSmqlrisTWiF+llLPw4RH07hKA3mDEYDRxPiGNddtPWtSLVynl/HX+OIb363zT89kaddflueTkFXLk1AU27zhCbkGRxXES4MlHrVMpIjfedNT1M3Yn2KsFXVlZGaNGjeLo0aO4ut5+2WwxQq8HN/MJj4zow8iIPpSVq/nX6o38EnnR4tgKjZbT5+IZNqCHRR/H5kBdPPA3a+pRndpG/GMjulnsG+TrxekLKVb7RfQJrdP5ao66zZhAry1DXVqAXutMTl4hx05H8/OJc8Reqmxa4O/TFrlchl5vaTMMC7HuqCNshg3D4ch4DpyMv+k+xWVqUrMKzfMgao2OC5cyeH75ZjxcaxfcCcO617mInb1a0KWmpuLl5cWqVauIjIzE1dWVF154gUGDBtUp7iqEoNeTW32BXV2cmDhmCGfOJ1jZIH/ce5yoC5cYPawvY4b3IzTYD6PR1CJy7dCwwt/Q+8GNO6zYS1dQa3RIJBKzOISFBPL07PsZPbQv/j7tbN6JiVSKfVFrdFaT2iaTCbVGd1NBrw/2akFnMBhITU2lR48e/P3vf+f8+fM8++yzHDhwADc3tzqfRwh6IxDRL5zwsI4WgtClUyDj7hrA0VMX+N8PB/lm6wE6BPig0+spuF6CVqdr1rn2+lDXVa8Nud/V9Gz+9+NBzsYkmUXBZDIhl8t44clpTJkw3GJ/sWKzaRkbcetRdG3ptYUPj2iQ1c32bEHn7++PXC43X6tv3760adOGlJSUejXAEILeCNwsNfPgvXdRcL2Eo5EX2L7vONcyCszHVai1xCRe4dTZOO4aJG7h64otR5FUKiHhchrHIi/yy+kLXEvPqeVYA9eLy6y2i1SK41Hv9Fo9sWcLurZt2xIREcHx48e5++67SUlJIT8/n44dO9brOYhJUTuybst+1m7aa+Xxdnd15oFxQxk1tA/dOwcjlUpbjA2yoalZWE2hkOPp7gomyC0oQiaV0rdHZ0YMqSyA9fFXP1ikwprzZHVrpK5F6m6HDz/8kKKiIuLi4swt6N544w28vLxu+5zJycksWbKE4uJicwu6Tp06AdYt6FJTU3n11Ve5fv06crmcF198kVGjRtXrekLQ7YitkgMKuYyQYD8uX8vEYDDi3c6Luwf3JDrhCmmZuc3eBtnQHDp+ln99stHKWRQe1oGp9wxn+KCelQJP7S4l8ToKoGW0oBMpFztiK9deJTDlFWpORMXyS+QFdhw4ib7a8vIKtZZYB6sG2RjUlkq5mp7Dqd9jORkVy/m4y1YTZRIJDB/Uk4ljhlhsF9UMBTcjNjaWXr2ad09WMUK3M3XxMX+5cTfrfzhodWxwgDePTBrFsAE98G7nZXXO5pyesUqlyOV4ergil8nIzKmcd+jUwZ/gAB9O/R6LRntjokykUQStFTFCtzN1mXzr0bWj9QpVmZTSsgo+/HwLUFl+YPjAngwdEM4X3+5y6FWqt/rBMRiM/LDnGBfiL5u94Fqdntz8Irp3DmbW1LEM7R+Or3ebWtMowmIoaI0IQW8G1JaaeX/pM6Rl5HIiKoYTUbGs33qAdVv2WxzraM06bHeH6sBz8yZzPu4yZ6OTOB93mdKyCqtjJRK4a3Avpt5zw2Io0ihNg06nY+bMmUybNo158+bZ3Gf69OlMmTKFxx9/vGmDE5gRKZdmQl1SM9eLS/nw8++tVqlCZXpi6j3D6d09lNBgv0ZxztTlfCeiYnjr/623aihSRaBvO/r1CsPNxYnt+0+KwlcOwtdff82qVas4duwYzs7ONvfZsWMHr732GgcPHsTHx6eJIxSAGKE3G+qSmvHycLO5SlUqlZKbf53/9+VWANxcnenVNYTUjBzyCotvuaipLkJta+TdvXMH/vzYFK6kZXMlNZuU1EwuxF22KeYjI3rzl8cfxKd9G/P5LqVkiFSKA2Aymfj666+ZMGECSqWSd955h507d1JaWkpgYCBLlixhzJgxjB8/nn/84x9s2LCBl156yd5ht0rECL2FUVtO+f2lz5CTd50LcZe5EH+Z02fjrQpQAfj7tqVjgC9tvdxp6+VBG083dh46RUZ2PhqtDqVCTnCAN3OmjUOnN6DR6NBodVxKSePw8XMWbdaqI5NKCQ7wxt3NhbikaxZ1UmobeYvCV45BXFwcDz74IO+88w5hYWEsXryYKVOm0KNHD5YtW0a7du344YcfAJg1axalpaXs2LHDzlG3TsQIvYVxs5xygG87Anzbcd/owZWLmjbvtSo9K5fKKCwqIelKOoVFpRiMRovHtTo9yVczefuj/9UpntHD+jL/oQkEB/igVMjrNYkpVms6BteuXQPA19eX/v37s2bNGo4dO8aePXsoLS01F5+q2icuLs5eobZ6hKC3QOoihF1CA3FSKa1WTf7psSnm44xGI59/u4vvtv9stZp12n138/ADI1EpFahUSs5GX+Kfq761Ot+9owbRueONKoViErP5YjQa2b9/Py+++CKPPfYYs2bNIisri8zMTPM+Eomk2VUTbUmIb1Erpco546xSIpFUim/NkbJUKqVPeCecnJQWxzqrlAzp141Av/a0b+uJu6szdw3qdcvzVVH1gzP/oQkMH9hTiLmDExpaWbI4JyeH06dPYzAY8PLyIj4+nujoaAzV0mzZ2dkWxa0ETYsYobdS6jpSrs0yWVOoxci75dK1a1cCAwOJiori2Wef5ffff2fVqlUEBwfTv39/zpw5Q0FBAW5ubsTFxbFgwQJ7h9xqEZOiglsiJicF69ev55NPPuGXX35BqVTa3OfgwYP87W9/48CBA3h7ezdKHM1tFXRde4r+/PPPfPzxx5hMJkwmE4sWLeKee+6p9/UcQtAzMjJYvHgx+fn5hIaG8sEHH1i1YcrLy2Pp0qWkpaXh6urKkiVL6N+/f53OLwRdILgz9Ho9jzzyCLNmzeLhhx+2uc+cOXMYN24cTz75ZKPEYHtRmmOtgq7J/Pnzeeihh8w9Rbdu3WrVU9RkMjFkyBA2bNhA165diY+PZ9asWURFRSGtZyVJh0i5vPXWW8yePZsHHniATz75hNWrV7N48WKLff71r3/Ro0cP1qxZQ2pqKk888QQ7d+68rV5/AoGgfsjlcrM1sTY2bNhw2+ffe/QMew6fvuk+RSVlXE3LxvjHGLRCreVcTBILFn9grqhpi4ljh3DfqMG1Pl6dhu4pGhsby9q1awGYNGkSy5Yto6CgwKomulQqpaSkBICSkhJ8fHzqLebgAIKu0+k4c+YMn3zyCVC5fHju3LlWgh4XF8fChQsBCA4OxsvLi7Nnz1qVuiwuLqa4uNhiW3p6OgBZWVmN9TQEAsFN8PPzQy6/M7mpUGvMYl6F0WSiQq29qaDXB3v0FJVIJHz00Uf86U9/wsXFhbKyMj7//PPbit/ugl5YWIibm5v5zfb29iY7O9tqvx49erBr1y5efPFFEhMTSUpKIi8vz2q/devWsWrVKpvXmjNnTsMGLxAI6sSt0p33jRp8y1G0rf4BziolLyyY1mBrFezRU1Sv1/PZZ5+xevVqBg4cSFRUFC+++CK7du2ySj3fiiYV9D179rB8+XKLbbYmCGz5WF955RWWLVvG5MmT6du3LxERESgUCqv9HnvsMYtfUwCtVktqaiohISEWiyBqkpWVxZw5c9iwYQN+fn51fFaCxkS8J45Jfd+Xhnjv6uq4ul3s1VM0Li6OnJwcBg4cCMDAgQNxdnYmOTmZPn361Os5NKmgT5w4kYkTJ1ps0+l0REREmJ90bm6uzcI+5eXlLFu2zNwBe/LkyXToYN1L0MPDAw8PD6vtVW2f6oKfn5+YPHUwxHvimDTl+9LY1lh79RT18/MjKyuLy5cv06lTJ5KTk8nPz7epb7fC7ikXhULBoEGD2L17N5MnT2bbtm2MHDnSar///e9/tG/fnqeeeopff/0VnU5H9+437xIuEAhaFo1ZDiImJoYZM2bw0ksvmXuKrly58o7O+eabb7JkyRJWr15t7ilaRfWeom+++SYvvPCCOTvxz3/+87Z6mTqEbTE9PZ0lS5aQn5+Pv78/K1aswNPTk40bN5KTk8MLL7xAQUEB//d//0deXh6urq7mQkENibA3Oh7iPXFMWuL7InqKNhCBgYGsX7/eavusWbPM/2/btq3Z/iMQCAQNTUvoKeoQgu4oeHh4sGjRIps5eIF9EO+JY9IS35fIyEh7h3DHOETKRSAQCAR3jmOulxUIBAJBvRGCLhAIBC2EVivoO3bs4P7772fChAk2a1DExcXx0EMPce+99/Laa6+h1+vtEGXr4lbvyapVqxgzZgxTp05l6tSpd1Q7RFB3SktLmTRpEmlpaVaPie+Jg2FqhWRlZZnGjBljKiwsNJWVlZkmT55sunTpksU+DzzwgOns2bMmk8lkeuWVV0wbNmywQ6Sth7q8JwsXLjT9/vvvdoqwdXLu3DnTpEmTTD179jSlpqZaPS6+J45FqxyhnzhxgqFDh+Ll5YWLiwv33nsve/fuNT+enp6OWq2mX79+QGXBsOqPCxqeW70nULmS74svvmDy5Mm8/fbbaDQaO0Xbeti8eTNvvPGGzdXb4nvieLRKQc/JybEowO/j42NREKzm47UVDBM0HLd6T8rKyggPD+fvf/87P/74I8XFxaxevdoeobYq3n33XQYNGmTzMfE9cTxapaCbbDg1qxcEu9XjgobnVq+5q6srX3zxBR07dkQul/Pkk09y9OjRpgxRUAPxPXE8WqWg+/r6WpTezcnJsbilrPl4bQXDBA3Hrd6TjIwMtmzZYv7bZDLdcX1twZ0hvieOR6sU9OHDh3Py5EkKCgqoqKhg//79FgXBAgMDUalUREVFAdRaMEzQcNzqPXFycuL9998nNTUVk8nEhg0bmDBhgh0jFojviePRKgXd19eXl156ifnz5/Pggw8yadIk+vTpw9NPP83FixcB+OCDD1i+fDkTJ06koqKC+fPn2znqls2t3pO2bdvy9ttv89xzz3HfffdhMpl44okn7B12q0R8TxwXsfRfIBAIWgitcoQuEAgELREh6AKBQNBCEIIuEAgELQQh6AKBQNBCEIIuEAgELQQh6AKBQNBCEIIuEAgELQSxdlrgkOj1ej777DO2bt1KWVkZS5cuJSsrC71ez3PPPWfv8AQCh0QIusAh+eijj4iOjmb79u2cOXOG999/H6lUyqZNm+wdmkDgsAhBFzgcpaWlrFu3jt27d+Pu7k7fvn25fPkyL730Em5ubvYOTyBwWEQOXeBwnDp1ipCQEIKDgwHQ6XS4u7szd+5cO0cmEDg2QtAFDkd2drZFGdZNmzbh6+srRucCwS0QKReBw+Hn50d8fDw5OTlkZmayfft2ysvL0Wq1KJVKe4cnEDgsotqiwOHQarUsXbqUw4cP4+npycqVK3n//fdRq9Vs3LjR3uEJBA6LEHSBQCBoIYgcukAgELQQhKALBAJBC0EIukAgELQQhKALBAJBC0EIukAgELQQhKALBAJBC0EIukAgELQQhKALBAJBC+H/A2WB0l7UCnPpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 288x194.4 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(4, 2.7))\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\"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": "Python 3.10.6 ('venv': venv)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "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
}