{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "9b2ffbd7",
      "metadata": {},
      "source": [
        "# Sampling from the *Hy-MMSBM* generative model"
      ]
    },
    {
      "cell_type": "raw",
      "id": "69ba0207",
      "metadata": {
        "raw_mimetype": "text/restructuredtext"
      },
      "source": [
        ".. admonition:: Definition\n",
        "\n",
        "   Hy-MMSBM assigns each node a membership vector and samples hyperedges using an affinity matrix.\n"
      ]
    },
    {
      "cell_type": "raw",
      "id": "1368c6ed",
      "metadata": {
        "raw_mimetype": "text/restructuredtext"
      },
      "source": [
        ".. admonition:: You will learn\n",
        "\n",
        "   Sample synthetic hypergraphs from the Hy-MMSBM model and inspect communities.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Overview\n",
        "\n",
        "- Sample synthetic hypergraphs from the Hy-MMSBM generative model.\n",
        "- Explore parameters and inspect the resulting community structure.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Setup\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import matplotlib as mpl\n",
        "\n",
        "mpl.rcParams.update({\n",
        "    \"figure.figsize\": (6, 4),\n",
        "    \"figure.dpi\": 120,\n",
        "    \"savefig.dpi\": 150,\n",
        "})\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "84e3a3e9",
      "metadata": {},
      "outputs": [],
      "source": [
        "import sys\n",
        "from collections import Counter\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "sys.path.append(\"..\")\n",
        "from hypergraphx.communities.hy_mmsbm.model import HyMMSBM\n",
        "from hypergraphx.core.hypergraph import Hypergraph\n",
        "from hypergraphx.generation.hy_mmsbm_sampling import HyMMSBMSampler\n",
        "from hypergraphx.linalg.linalg import binary_incidence_matrix\n",
        "\n",
        "np.random.seed(123)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e4165b9a",
      "metadata": {},
      "source": [
        "## Vanilla sampling from the generative model"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "782341db",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAHqCAYAAAB/Z/s1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwsUlEQVR4nO3de1zUdb7H8fcMCkgI6qKgOEam5aoJJ7yElVqHMrcsc901u6Bs2VpiFp1OuaW41kZtm0u7WlqbuW6ammUXLfYYZ7HN7FiadvOSaYoaqCkXUUCZ3/nDZVYEdIYZGWa+r6eP32OX3+X7+8yQ8/Hzme/v97NZlmUJAIAgZ/d3AAAANAUSHgDACCQ8AIARSHgAACOQ8AAARiDhAQCMQMIDABiBhAcAMAIJDwBgBBIeAMAIJDwAgBFIeJAkFRYW6q677pLD4VDLli1ls9k0ZMgQt7YPGTJENptN8+fP92lM52pcnBnvO4JVC38HAN+bPXu2MjIyJEn33HOPnn/++TPuf+LECV199dXavHmzJKlt27YKDQ1Vu3bt3Nre1HJyclRcXKxx48YpISHBLzE0N7wnwNmR8ILQggULXP9/yZIlysnJUWhoaIP7//3vf9fmzZvVrl07ffLJJ+revbtH27t06aKLL75Y0dHRPn0dDY2bk5OjXbt2aciQIXy4/4sv35Nz9fsE/I2EF2S2bdumdevW6fzzz1f37t31wQcfaOXKlbr55psbPObrr7+WJF111VV1kpk7209NsL50rsbFmfG+I1jxHV6QqfmwuuWWW3TbbbdJkv72t7+d8Zhjx45JkiIjIxu1HQACgoWg4XQ6rfPPP9+SZG3cuNEqLi62wsPDrdDQUOvHH3+ss39WVpYlqcFl7NixZ9y+c+dOy7Isa/DgwZYk65VXXqlzjpp4/vGPf1gHDhywJk2aZHXp0sUKDQ21unTpYk2ePNk6fPhwva/n9HFfeeWVM8YzePBgy7Isa/fu3ZbdbrckWd9++22D79fmzZstSVZISIi1d+9et97jU1/Pnj17rLvuusuKj4+3WrVqZV1yySW13oPq6mpr1qxZVlJSkhUREWHFxMRY48aNs4qKiuode/fu3dZTTz1lXXPNNVbXrl2tsLAwKzo62kpJSbH+9Kc/WZWVlXWOcfc9OT32Xbt2WePHj7e6dOlitWjRwrrpppsafN8ty7IefPBBS5LVuXNnq7i4uE4cR48etXr06OH678YT5513niXJ+u677+psGzJkiCXJOu+88+psO3LkiBUVFWXZbDZr+/btHp0TZiLhBZH8/HxLktWzZ0/XupEjR1qSrOeff77O/s8884wVGxvr+sAJDw+3YmNjXcuMGTPOuH337t2WZbmX8BYsWGA5HA5LkhUZGWm1bNnS9aF86aWX1vthfvq4ixcvtmJjY13JrG3btrXiufnmm13HDh061JJkPfroow2+Xw8//LAlyRo2bJi7b7Hr9cybN8+KjY21JFnR0dGWzWZzvZ6nn37acjqd1qhRoyxJVlhYmNWqVSvX9l69elnHjh2rM/bPf/5z1z7h4eFW27ZtayWvIUOG1HmfPHlPamJ/4YUXrHbt2rl+F+Hh4WdNeBUVFVavXr0sSdbtt99eJ/b77rvPkmR16dLFKikpcfv9tCzLio+PtyRZmzZtqrX+yy+/rPX6jx8/Xmv7iy++6PHvD2Yj4QWRO++805JkPf744651y5YtsyRZKSkpDR5XU+k19C/zs213J+G1adPGSk5Otj799FPLsiyrsrLSWrBggRUeHm5Jsv785z+7Pe6plUpDli5d6voArq6urrP9xIkTVqdOnSxJ1pIlSxocp6HXEx0dbV1xxRXWN998Y1mWZZWWllqTJk2yJFmtWrWyHn30USsqKsp67bXXrKqqKqu6utpauXKlFRUVZUmynnvuuTpjT5061Zo9e7b13XffWU6n07Ksk5XT0qVLXUnh1N+tp+9JzT6RkZFWUlKStW7dOsuyTnYGTq2QGnrfN2zY4PqHyhtvvOFa/8EHH1g2m82y2WxWXl6eW+/jqWoS6Zo1a2qtHz9+vCXJat++vSXJOnDgQK3tycnJliRr5cqVHp8TZiLhBYljx45Z0dHRddp4x44dc33INtTea4qEFx0dXecDy7Isa/LkyXVab2cb150P98rKSusnP/mJJclatWpVne25ubmuiqiioqLBcRp6Pe3atavT2nM6ndZFF13kqkheffXVOsc//vjjDb7eM1mzZo2rpXimuNxJeG3btrX279/f4H5n+n0+8cQTliQrJibGKiwstIqLi12V+3333efRa6px+eWXW5Ks3Nxc17pDhw5ZERERVvfu3a1x48bV+e933bp1liTrwgsvdP3jADgbJq0EibffflslJSXq16+funXr5lofHh7umqF5tskr59KECRMUExNTZ/3w4cMl/XsmqK+EhobqjjvukKR6L6CuWTdmzBiFhYV5PP6ECRPqTNs/9WJ8h8OhMWPG1DnuP//zPyV5/noHDhyoNm3aaM+ePdq7d6/H8Z4qLS1N7du3b9SxjzzyiC677DIdPHhQd999tzIyMlRQUKCLL75YTz31VKPGbNOmjSTpyJEjrnUvv/yyjh49qokTJ7re58OHD7u2z5kzR5J07733ymazNeq8MA8JL0jUzM6s70O2Zt2rr77apDGdqk+fPvWuj4+PlyQVFxf7/Jy/+tWvJEnLly9XaWmpa31JSYneeustSVJ6enqjxu7du3e96zt06CBJ6tmzp+z2un+9arY39Ho/+ugj3X777brwwgsVEREhm83mWmqO+eGHHxoVc43LLrus0ceGhIRowYIFioiI0DvvvKNXX31VLVq00IIFC9SqVatGjVmT8MrKyiRJTqdTs2fPVmRkpNLT0xUVFSXp3+9ZSUmJFi9erIiICNfvGHAHCS8IFBUV6X/+539kt9s1evToOttTU1PVoUMH7dixQ2vWrPFDhFKnTp3qXR8eHi7p5N1cfO2SSy5Rv379dPToUS1dutS1fvHixaqoqFDv3r3Vt2/fRo3dsWPHeteHhIRIkuLi4s64vb7X++STT+rKK6/UwoULtWPHDp04cULt2rVTbGysYmNjXQm0vLy8UTHXaGx1V6N79+565JFHXD8/+OCD6t+/f6PHO73Ce/fdd/X9998rLS1NUVFRdRLeggULdPToUd12222uYwF3kPCCwGuvvaYTJ07I6XQqPj6+VlVgs9nUokUL7d+/X5J5FxXfeeedkmq3Nf/6179KksaNG+eHiOr35ZdfaurUqZJOtum2bt2qiooK/fjjjyosLFRhYaHrHw2WZXl1rpqk21iVlZVavHix6+e1a9fK6XQ2erzTE96f//xnSXLdHq8m4dW0NOfOnVtrO+AuEl4Q8CSJLV26VJWVlecwmuZlzJgxioiI0Jo1a7R9+3Zt27ZNa9euVYsWLXT77bf7OzyXN998U06nU4MHD9bs2bN10UUX1WqJVldX68cff/RjhP/26KOP6ptvvlF8fLzatGmjDz/8UDk5OY0e79SW5jfffKO8vDylpqbqpz/9qSTVqvD++c9/6uuvv9aVV17ZYJscaAgJL8B9/fXX+vzzz2W327V9+3YdPny4waVz584qLi7WihUr/B22V2oSgTuVTlRUlEaNGiXpZJVXU+kNGzZMsbGx5yxGT9VMRGmoxfp///d/rjve1MeT98QbH374of74xz/KZrPplVde0XPPPSfpZBKsubm4p06t8Gqqu0mTJrm2t27dWtLJhFczWYXqDo1BwgtwNdXd5ZdfrgsvvFBt2rRpcBkxYkStYwLV6d/pnE3NxIYFCxa4Zqo2drLKuVLzmrZu3Vpnm9Pp1PTp0906/lxM/qlx5MgRjRs3Tk6nU/fcc4+uueYapaWlacSIEaqoqFBaWlqjvoutSXh79uzR3/72N11wwQW64YYbXNtrXtu3336rN954Q506ddLIkSN98ppgFhJeAHM6nVq4cKEkufUBULPP+++/32zaY43Rq1cvSSe/u6yoqDjr/oMGDVK3bt1UUFCgPXv2KCYmptYHanNQc7nCihUrNHPmTFfb+fvvv9cvf/lLffjhhzrvvPMaPN7T96QxMjMztXPnTnXr1k3PPPOMa/3cuXPVvn17ffbZZ3ryySc9Hrcm4b3zzjsqLy/XvffeW6udW5Pw3nrrLVVWVurXv/61WrTgvvfwHAkvgP3v//6vqxXmTsIbNGiQYmJidPz48VqTDgJNzUSU119/XdHR0XI4HEpISNAtt9xS7/42m63W9PXbbrtNLVu2bJJY3XXddddp2LBhkk7OeoyIiFDbtm11wQUXaPny5frTn/5U73WMNTx9Tzz1/vvv66WXXqp1WUKNDh06uFqNTzzxhDZs2ODR2DUJ78SJE4qIiHC9lho1Ce/EiRNq2bKl7r77bi9eCUxGwgtgNe255ORkdenS5az7h4SE6MYbb5QU2G3Nq6++WsuXL9fgwYPVqlUr7d27V7t27VJhYWGDx5z6eKTm1s6UTiblt956S1lZWerWrZtCQkLUokULDRs2TKtWrTrrh3xj3hN3HTp0yJWEHnroIaWkpNTZZ+TIkbr99tt1/PhxpaWleTQx6tRLC2677Ta1bdu21vaahCdJo0aNavCSD+BsbNa5/pYbaAbmzJmje+65R0lJSfr888/9HQ4AP6DCgxFefPFFSdJdd93l50gA+AsJD0Fv1qxZ+vzzz9WmTRvX/TUBmIepTghKe/bs0RVXXKGysjIdOnRIkjR9+vRa3wcBMAsVHoLSiRMntGvXLpWUlOjCCy/UH//4R02ePNnfYQHQyRsYDB8+XJ06dXJN2Dqb/Px8XXrppQoLC1O3bt3qfQrK2VDhISglJCSc87uOAGic8vJyJSYm6le/+pVbl1Tt3LlT119/vSZMmKCFCxcqLy9Pd911lzp27KihQ4e6fV5maQIA/MZms2n58uWuO0HV5+GHH9bKlSv11VdfudbdcsstKi4uVm5urtvnCsoKz+l0at++fWrdujUPhwQQFCzLUllZmTp16lTvsxYbo6KiQlVVVT4Zy7KsOp+3YWFhjXrA8unWrl2r1NTUWuuGDh2q+++/36NxgjLh7du3Tw6Hw99hAIDPFRQUqHPnzl6PU1FRoVatfyKdOOqDqKTIyMhaT62XpKysrLPeB9YdhYWFdW72Hhsbq9LSUh07dszthw8HZcKrubt6aM+xsoWE+jkanG53/h/8HQIQcMpKS9XtAofr881bVVVV0omjCus5VvL2c7K6Ske++asKCgpqzYT2RXXnS0GZ8GrKaltIKAmvGeLSAKDxfP41TYtwrz8nLdvJFuupT6j3pbi4OBUVFdVaV1RUpKioKLerOylIEx4AwE02Sd4m0XM8VSIlJUXvvfderXWrVq2q976uZ8J1eACAJnXkyBFt3LhRGzdulHTysoONGzdq9+7dkqQpU6YoLS3Ntf+ECRO0Y8cO/fd//7e2bNmi559/XkuXLtUDDzzg0Xmp8ADAZDb7ycXbMTzw2Wef6aqrrnL9nJmZKUkaO3as5s+frx9++MGV/CTpggsu0MqVK/XAAw/oueeeU+fOnfWXv/zFo2vwJBIeAJjNZvNBS9Oz44cMGXLGG0PUdxeVIUOGeP2kE1qaAAAjUOEBgMn80NL0FxIeAJjMDy1NfwmMtAwAgJeo8ADAaD5oaQZI7UTCAwCT0dIEACC4UOEBgMmYpQkAMAItTQAAggsVHgCYjJYmAMAItDQBAAguVHgAYDJamgAAI9hsPkh4tDQBAGg2qPAAwGR228nF2zECAAkPAExm0Hd4gRElAABeosIDAJMZdB0eCQ8ATEZLEwCA4EKFBwAmo6UJADACLU0AAIILFR4AmMygliYVHgDACFR4AGAyg77DI+EBgMloaQIAEFyo8ADAaD5oaQZI7UTCAwCT0dIEACC4UOEBgMlsNh/M0gyMCo+EBwAmM+iyhMCIEgAAL1HhAYDJDJq0QsIDAJPR0gQAILhQ4QGAyWhpAgCMQEsTAIDgQoUHACajpQkAMIHNZpPNkIRHSxMAYAQqPAAwmEkVHgkPAExm+9fi7RgBgJYmAMAIVHgAYDBamgAAI5iU8GhpAgCMQIUHAAYzqcIj4QGAwUxKeM26pTl79mwlJCQoPDxcAwYM0Lp16/wdEgAgQDXbhLdkyRJlZmYqKytLGzZsUGJiooYOHar9+/f7OzQACB42Hy0BoNkmvJkzZ2r8+PFKT09Xz549NWfOHEVERGjevHn+Dg0AgkZNS9PbJRA0y4RXVVWl9evXKzU11bXObrcrNTVVa9eu9WNkAIBA1SwnrRw8eFDV1dWKjY2ttT42NlZbtmyps39lZaUqKytdP5eWlp7zGAEgGJx8OpC3k1Z8E8u51iwrPE9lZ2crOjratTgcDn+HBAABwSYftDQDJOM1y4QXExOjkJAQFRUV1VpfVFSkuLi4OvtPmTJFJSUlrqWgoKCpQgUABIhmmfBCQ0OVnJysvLw81zqn06m8vDylpKTU2T8sLExRUVG1FgDA2Zk0aaVZfocnSZmZmRo7dqz69u2r/v37KycnR+Xl5UpPT/d3aAAQPAx6PFCzTXijR4/WgQMHNG3aNBUWFiopKUm5ubl1JrIAAOCOZpvwJCkjI0MZGRn+DgMAgpcPWpIWLU0AQHPni+/gAuU7vGY5aQUAAF8j4QGAwfw1S9PThwPk5OTo4osvVqtWreRwOPTAAw+ooqLCo3OS8ADAZH64ebSnDwdYtGiRHnnkEWVlZWnz5s16+eWXtWTJEv3mN7/x6LwkPABAk/L04QAff/yxLr/8ct16661KSEjQtddeqzFjxnj8yDgSHgAYrKlbmo15OMDAgQO1fv16V4LbsWOH3nvvPf3sZz/z6LUySxMADObLWZqn37g/LCxMYWFhtdZ5+nAASbr11lt18OBBXXHFFbIsSydOnNCECRNoaQIA/MPhcNS6kX92drZPxs3Pz9eTTz6p559/Xhs2bNCbb76plStX6vHHH/doHCo8ADCYLyu8goKCWvcyPr26kzx/OIAkTZ06VXfccYfuuusuSdIll1yi8vJy3X333Xr00Udlt7tXu1HhAYDBfPkd3uk38a8v4Xn6cABJOnr0aJ2kFhISIkmyLMvt10qFBwBoUmd7OEBaWpri4+NdLdHhw4dr5syZ+o//+A8NGDBA27dv19SpUzV8+HBX4nMHCQ8ATOaHpyWc7eEAu3fvrlXRPfbYY7LZbHrssce0d+9etW/fXsOHD9fvfvc7z8K0PKkHA0Rpaamio6MVdsl42UJC/R0OTnP401n+DgEIOKWlpYr9SbRKSkp88szPms/JuF+9KntohFdjOauOqnDe7T6L7VzhOzwAgBFoaQKAwXhaAgAAQYYKDwAMZlKFR8IDAJP5YZamv9DSBAAYgQoPAAxGSxMAYASTEh4tTQCAEajwAMBgNvmgwguQWSskPAAwGC1NAACCDBUeAJjMoOvwSHgAYDBamgAABBkqPAAwmEkVHgkPAAxms51cvB0jENDSBAAYgQoPAAx2ssLztqXpo2DOMRIeAJjMBy3NQLksgZYmAMAIVHgAYDBmaQIAjMAsTQAAggwVHgAYzG63yW73rkSzvDy+qZDwAMBgtDQBAAgyVHgAYDBmaQIAjEBLEwCAIEOFBwAGo6UJADCCSQmPliYAwAhUeABgMJMmrZDwAMBgNvmgpRkgzweipQkAMAIVHgAYjJYmAMAIzNIEACDIUOEBgMFoaQIAjEBLEwCAIEOFBwAGo6UJADACLU0AAIJMUFd4u/P/oKioKH+HgdO07Zfh7xBwBoc/neXvENCUfNDSDJA7iwV3wgMAnBktTQAAggwVHgAYjFmaAAAj0NIEACDIUOEBgMFoaQIAjEBLEwCAIEOFBwAGM6nCI+EBgMFM+g6PliYAwAhUeABgMJNamlR4AAAjUOEBgMFM+g6PhAcABqOlCQBAkKHCAwCD2eSDlqZPIjn3SHgAYDC7zSa7lxnP2+ObCi1NAIARqPAAwGDM0gQAGIFZmgAAnEOzZ89WQkKCwsPDNWDAAK1bt+6M+xcXF2vixInq2LGjwsLCdNFFF+m9997z6JxUeABgMLvt5OLtGJ5YsmSJMjMzNWfOHA0YMEA5OTkaOnSotm7dqg4dOtTZv6qqStdcc406dOigZcuWKT4+Xrt27VKbNm08Oi8JDwBMZvNBS9LDw2fOnKnx48crPT1dkjRnzhytXLlS8+bN0yOPPFJn/3nz5unQoUP6+OOP1bJlS0lSQkKCx2HS0gQA+ERpaWmtpbKyss4+VVVVWr9+vVJTU13r7Ha7UlNTtXbt2nrHfeedd5SSkqKJEycqNjZWvXv31pNPPqnq6mqP4iPhAYDBamZpertIksPhUHR0tGvJzs6uc76DBw+qurpasbGxtdbHxsaqsLCw3hh37NihZcuWqbq6Wu+9956mTp2qZ599Vk888YRHr5WWJgAYzPavP96OIUkFBQWKiopyrQ8LC/Nq3BpOp1MdOnTQiy++qJCQECUnJ2vv3r165plnlJWV5fY4JDwAgE9ERUXVSnj1iYmJUUhIiIqKimqtLyoqUlxcXL3HdOzYUS1btlRISIhr3U9/+lMVFhaqqqpKoaGhbsVHSxMADFYzS9PbxV2hoaFKTk5WXl6ea53T6VReXp5SUlLqPebyyy/X9u3b5XQ6Xeu2bdumjh07up3sJBIeABit5sJzbxdPZGZm6qWXXtJf//pXbd68Wffcc4/Ky8tdszbT0tI0ZcoU1/733HOPDh06pMmTJ2vbtm1auXKlnnzySU2cONGj89LSBAA0qdGjR+vAgQOaNm2aCgsLlZSUpNzcXNdElt27d8tu/3c95nA49Pe//10PPPCA+vTpo/j4eE2ePFkPP/ywR+cl4QGAwfx1L82MjAxlZGTUuy0/P7/OupSUFH3yySeen+gUJDwAMBiPBwIAIMhQ4QGAwXg8EADACDweCACAIEOFBwAGo6UJADACszQBAAgyVHgAYDCbPH5+a71jBAISHgAYjFmaAAAEGSo8ADCYp4/3aWiMQEDCAwCD0dIEACDIUOEBgOECpEDzGgkPAAxGSxMAgCBDhQcABmOWJgDACLQ0AQAIMlR4AGAw7qUJADACjwcCACDIUOEBgMF44jkAwAjM0gQAIMi4lfB+97vfnes4AAB+UNPS9HYJBG4lvKysLL355pvnOhYAQBOrmaXp7RII3Ep44eHhSktL0/r16z0+QV5ensfHfPjhhxo+fLg6deokm82mt956y+MxAAA4lVsJb+HChaqoqNBNN92kffv2uT3473//ew0bNszjoMrLy5WYmKjZs2d7fCwAwH0mtTTdmqV50003KTs7Ww8//LBuuOEGffTRR4qIiGhw/6NHjyo9PV3Lli1TSEiIx0ENGzasUYkSAOAZZmnW46GHHlJ6ero2btyoW2+9tcH9vvvuO1122WV6/fXXFRMTo9zcXJ8EeiaVlZUqLS2ttQAAcCqPLkuYO3euBg0apHfffVcPPfRQne0rV65Uv3799NVXX6l///5av369rr76ap8F25Ds7GxFR0e7FofDcc7PCQDBwO6jJRB4FGeLFi20fPlyde3aVTNnztTLL7/s2jZjxgzddNNNKi4u1oQJE/TPf/5TnTt39nnA9ZkyZYpKSkpcS0FBQZOcFwACXU1L09slEHh8p5W2bdtqxYoVSklJ0b333qt27dpp/vz5WrFihcLDwzVnzhzdcccd5yLWBoWFhSksLKxJzwkACCxuJbwZM2YoKSlJSUlJ6tKliy6++GItW7ZM1113nUaNGiXLstS1a1e9+eab6tOnz7mOGQDgIzYfPPE8QAo89xLe9OnTXSVrdHS0EhMTlZiYqCuuuEL5+fn62c9+pkWLFikqKsonQR05ckTbt293/bxz505t3LhR7dq1U5cuXXxyDgCAWdxKePfdd5+++OILbdq0SYcPH9bq1au1evVqV+/2gw8+0ODBg9WnTx8lJSW5EuJPfvKTRgX12Wef6aqrrnL9nJmZKUkaO3as5s+f36gxAQB12X1Q4Xl7fFNxK+Hl5OS4/n9BQYE2btyoTZs2uf53x44d2rRpkzZt2qRXX33VtW/Hjh2VlJSkFStWeBTUkCFDZFmWR8cAADxn0nV4Hk9acTgccjgcGj58uGvdkSNH9MUXX9RKhF999ZX27dunH374wacBAwDQGD55Hl5kZKQGDhyogQMHutY5nU5t27ZNmzZt8sUpAADnAC1NH7Db7erRo4d69Ohxrk4BAPCSSU88D5QL5AEA8Mo5q/AAAM2fL55nFyjPwyPhAYDBfHEvzEBpFQZKnAAAeIUKDwAMZtKkFRIeABjMLh98h6fAyHi0NAEARqDCAwCD0dIEABjBpDut0NIEABiBCg8ADHbyAbDePi3BR8GcYyQ8ADCYSd/h0dIEABiBCg8ADGbSpBUSHgAYzPavP96OEQhoaQIAjECFBwAGo6UJADCCSQmPliYAwAhUeABgMJvNJpvXF54HRolHwgMAg9HSBAAgyFDhAYDBTLq1GAkPAAxmt/ngiecBkvFoaQIAmtzs2bOVkJCg8PBwDRgwQOvWrXPruMWLF8tms2nEiBEen5OEBwAGq5m04u3iiSVLligzM1NZWVnasGGDEhMTNXToUO3fv/+Mx33//ff6r//6L1155ZWNe62NOgoAEBxs//4er7GLp7fSnDlzpsaPH6/09HT17NlTc+bMUUREhObNm9fgMdXV1brtttv029/+Vl27dm3USyXhAQCaTFVVldavX6/U1FTXOrvdrtTUVK1du7bB42bMmKEOHTrozjvvbPS5mbQCAAazyya7l087qDm+tLS01vqwsDCFhYXVWnfw4EFVV1crNja21vrY2Fht2bKl3vE/+ugjvfzyy9q4caOXcQIAjOVtO/PUyxocDoeio6NdS3Z2ttfxlZWV6Y477tBLL72kmJgYr8aiwgMA+ERBQYGioqJcP59e3UlSTEyMQkJCVFRUVGt9UVGR4uLi6uz/3Xff6fvvv9fw4cNd65xOpySpRYsW2rp1qy688EK34iPhAYDBfHlrsaioqFoJrz6hoaFKTk5WXl6e69ICp9OpvLw8ZWRk1Nm/R48e+vLLL2ute+yxx1RWVqbnnntODofD7ThJeABgMH9ceJ6ZmamxY8eqb9++6t+/v3JyclReXq709HRJUlpamuLj45Wdna3w8HD17t271vFt2rSRpDrrz4aEBwBoUqNHj9aBAwc0bdo0FRYWKikpSbm5ua6JLLt375bd7vspJiQ8ADCYv+6lmZGRUW8LU5Ly8/PPeOz8+fM9P6FIeABgNLt80NL08rKGpsJlCQAAI1DhAYDBeDwQAMAIdnnf6guUVmGgxAkAgFeo8ADAYDabTTYve5LeHt9USHgAYLBGPN2n3jECAS1NAIARqPAAwGD+uLWYv5DwAMBwgZGuvEdLEwBgBCo8ADAYF54DAIxg0mUJtDQBAEagwgMAg3FrMQAAggwVHgAYzKTv8Eh4AGAwbi0GAECQocJDkzv86Sx/h4AzaNsvw98hoB5WddU5GZeWJgDACMzSBAAgyFDhAYDBaGkCAIzALE0AAIIMFR4AGIynJQAAjGCXTXYvm5LeHt9UaGkCAIxAhQcABqOlCQAwgu1ff7wdIxDQ0gQAGIEKDwAMRksTAGAEmw9madLSBACgGaHCAwCD0dIEABjBpIRHSxMAYAQqPAAwmEnX4ZHwAMBgdtvJxdsxAgEtTQCAEajwAMBgtDQBAEZgliYAAEGGCg8ADGaT9y3JACnwSHgAYDJmaQIAEGSo8ADAYMzSBAAYgVmaAAAEGSo8ADCYTd7PsgyQAo+EBwAms8smu5c9SW+fmN5UaGkCAIxAhQcABqOlCQAwg0EZj5YmAMAIVHgAYDAuPAcAmMEHF54HSL6jpQkAMAMVHgAYzKA5KyQ8ADCaQRmPliYAwAhUeABgMGZpAgCMwOOBAAAIMlR4AGAwg+askPAAwGgGZTxamgAAI5DwAMBgNh/98dTs2bOVkJCg8PBwDRgwQOvWrWtw35deeklXXnml2rZtq7Zt2yo1NfWM+zeEhAcABquZpent4oklS5YoMzNTWVlZ2rBhgxITEzV06FDt37+/3v3z8/M1ZswY/eMf/9DatWvlcDh07bXXau/evR6dl4QHAGhSM2fO1Pjx45Wenq6ePXtqzpw5ioiI0Lx58+rdf+HChbr33nuVlJSkHj166C9/+YucTqfy8vI8Oi8JDwAMZvPRIkmlpaW1lsrKyjrnq6qq0vr165WamupaZ7fblZqaqrVr17oV89GjR3X8+HG1a9fOo9dKwgMA+ITD4VB0dLRryc7OrrPPwYMHVV1drdjY2FrrY2NjVVhY6NZ5Hn74YXXq1KlW0nQHlyUAgMl8eFlCQUGBoqKiXKvDwsK8HLiup556SosXL1Z+fr7Cw8M9OpaEBwAG8+W9NKOiomolvPrExMQoJCRERUVFtdYXFRUpLi7ujMf+4Q9/0FNPPaUPPvhAffr08ThOWpoAgCYTGhqq5OTkWhNOaiagpKSkNHjc73//ez3++OPKzc1V3759G3VuKjwAMJg/bh6dmZmpsWPHqm/fvurfv79ycnJUXl6u9PR0SVJaWpri4+Nd3wE+/fTTmjZtmhYtWqSEhATXd32RkZGKjIx0+7wkPAAwmD/uLDZ69GgdOHBA06ZNU2FhoZKSkpSbm+uayLJ7927Z7f9uQL7wwguqqqrSqFGjao2TlZWl6dOnu31eEh4AoMllZGQoIyOj3m35+fm1fv7+++99ck4SHgCYzKCbR5PwAMBgJj3xvFnO0szOzla/fv3UunVrdejQQSNGjNDWrVv9HRYAIIA1y4S3evVqTZw4UZ988olWrVql48eP69prr1V5ebm/QwOAoOKPm0f7S7Nsaebm5tb6ef78+erQoYPWr1+vQYMG+SkqAAg+Bn2F1zwT3ulKSkokqcEbhVZWVta6SWlpaWmTxAUACBzNsqV5KqfTqfvvv1+XX365evfuXe8+2dnZtW5Y6nA4mjhKAAhQvnxcQjPX7BPexIkT9dVXX2nx4sUN7jNlyhSVlJS4loKCgiaMEAACl7+eeO4PzbqlmZGRoRUrVujDDz9U586dG9wvLCzsnNyVGwAQPJplwrMsS5MmTdLy5cuVn5+vCy64wN8hAUBQ8se9NP2lWSa8iRMnatGiRXr77bfVunVr141Co6Oj1apVKz9HBwDBw6RZms3yO7wXXnhBJSUlGjJkiDp27OhalixZ4u/QAAABqllWeJZl+TsEADCDQSVes0x4AICmwb00AQAIMlR4AGAyX9wLMzAKPBIeAJjMoK/waGkCAMxAhQcAJjOoxCPhAYDBmKUJAECQocIDAINxL00AgBEM+gqPliYAwAxUeABgMoNKPBIeABiMWZoAAAQZKjwAMJhNPpil6ZNIzj0SHgAYzKCv8GhpAgDMQIUHAAbjwnMAgCHMaWrS0gQAGIEKDwAMRksTAGAEcxqatDQBAIagwgMAg9HSBAAYgXtpAgAQZKjwAMBkBs1aIeEBgMEMyne0NAEAZqDCAwCDMUsTAGAEZmkCABBkqPAAwGQGzVqhwgMAGIEKDwAMZlCBR8IDAJOZNEuTliYAwAhUeABgNO8vSwiUpiYJDwAMRksTAIAgQ8IDABiBliYAGIyWJgAAQYYKDwAMZtLNo0l4AGAwWpoAAAQZKjwAMBj30gQAmMGgjEdLEwBgBCo8ADAYszQBAEZgliYAAEGGCg8ADGbQnBUqPAAwms1Hi4dmz56thIQEhYeHa8CAAVq3bt0Z93/99dfVo0cPhYeH65JLLtF7773n8TlJeACAJrVkyRJlZmYqKytLGzZsUGJiooYOHar9+/fXu//HH3+sMWPG6M4779Tnn3+uESNGaMSIEfrqq688Oq/NsizLFy+gOSktLVV0dLSKfixRVFSUv8MBAkrbfhn+DgH1sKqrVPnlSyop8c3nWs3nZOFB78crLS1VXEy027ENGDBA/fr106xZsyRJTqdTDodDkyZN0iOPPFJn/9GjR6u8vFwrVqxwrbvsssuUlJSkOXPmuB0nFR4AGKxmlqa3i7uqqqq0fv16paamutbZ7XalpqZq7dq19R6zdu3aWvtL0tChQxvcvyFBOWmlpmgtKy31cyRA4LGqq/wdAupR83vxdVOu1AefkzVjnD5WWFiYwsLCaq07ePCgqqurFRsbW2t9bGystmzZUu/4hYWF9e5fWFjoUZxBmfDKysokSd0ucPg5EgDwrbKyMkVHR3s9TmhoqOLi4tTdR5+TkZGRcjhqj5WVlaXp06f7ZHxfCMqE16lTJxUUFKh169ayBcoVkWdQWloqh8OhgoICvpNsZvjdNF/B9ruxLEtlZWXq1KmTT8YLDw/Xzp07VVXlm4resqw6n7enV3eSFBMTo5CQEBUVFdVaX1RUpLi4uHrHjouL82j/hgRlwrPb7ercubO/w/C5qKiooPiLG4z43TRfwfS78UVld6rw8HCFh4f7dMyzCQ0NVXJysvLy8jRixAhJJyet5OXlKSOj/glTKSkpysvL0/333+9at2rVKqWkpHh07qBMeACA5iszM1Njx45V37591b9/f+Xk5Ki8vFzp6emSpLS0NMXHxys7O1uSNHnyZA0ePFjPPvusrr/+ei1evFifffaZXnzxRY/OS8IDADSp0aNH68CBA5o2bZoKCwuVlJSk3Nxc18SU3bt3y27/90UEAwcO1KJFi/TYY4/pN7/5jbp376633npLvXv39ui8QXkdXrCprKxUdna2pkyZUm9PHP7D76b54neD05HwAABG4MJzAIARSHgAACOQ8AAARiDhAQCMQMJrpo4cOaKsrCxdd911ateunWw2m+bPn+/vsCDp008/VUZGhnr16qXzzjtPXbp00S9/+Utt27bN36EZ7+uvv9YvfvELde3aVREREYqJidGgQYP07rvv+js0NANch9dMHTx4UDNmzFCXLl2UmJio/Px8f4eEf3n66ae1Zs0a/eIXv1CfPn1UWFioWbNm6dJLL9Unn3zi8bVB8J1du3aprKxMY8eOVadOnXT06FG98cYbuvHGGzV37lzdfffd/g4RfsRlCc1UZWWlDh8+rLi4OH322Wfq16+fXnnlFY0bN87foRnv448/Vt++fRUaGupa9+233+qSSy7RqFGj9Oqrr/oxOpyuurpaycnJqqioaPBu/DADLc1mKiwszOMbo6JpDBw4sFayk6Tu3burV69e2rx5s5+iQkNCQkLkcDhUXFzs71DgZ7Q0AR+wLEtFRUXq1auXv0OBpPLych07dkwlJSV655139P7772v06NH+Dgt+RsIDfGDhwoXau3evZsyY4e9QIOnBBx/U3LlzJZ18esrIkSM1a9YsP0cFfyPhAV7asmWLJk6cqJSUFI0dO9bf4UDS/fffr1GjRmnfvn1aunSpqqurffbcNwQuvsMDvFBYWKjrr79e0dHRWrZsmUJCQvwdEiT16NFDqampSktL04oVK3TkyBENHz5czNEzGwkPaKSSkhINGzZMxcXFys3N9dmTqOF7o0aN0qeffsq1koajpQk0QkVFhYYPH65t27bpgw8+UM+ePf0dEs7g2LFjkk7+IwXmosIDPFRdXa3Ro0dr7dq1ev3115WSkuLvkPAv+/fvr7Pu+PHjWrBggVq1asU/TAxHhdeMzZo1S8XFxdq3b58k6d1339WePXskSZMmTVJ0dLQ/wzPWgw8+qHfeeUfDhw/XoUOH6lxofvvtt/spMvz6179WaWmpBg0apPj4eBUWFmrhwoXasmWLnn32WUVGRvo7RPgRd1ppxhISErRr1656t+3cuVMJCQlNGxAkSUOGDNHq1asb3M5fKf9ZvHixXn75ZX355Zf68ccf1bp1ayUnJ2vSpEm68cYb/R0e/IyEBwAwAt/hAQCMQMIDABiBhAcAMAIJDwBgBBIeAMAIJDwAgBFIeAAAI5DwAABGIOEBHrr//vtls9k0ffr0ere//fbbioyMlN1u1+OPP960wQFoEAkP8NCmTZskSUlJSXW2PfPMMxo5cqScTqeWLFmiqVOnNnF0ABrCzaMBD9WX8I4fP64JEyZo3rx5iouL09tvv63+/fv7KUIA9eFemoAHdu/erfPPP1/R0dEqLi6WJB06dEgjR47U6tWr1adPH61YsUIOh8O/gQKog5Ym4IGNGzdKkhITEyVJW7du1YABA7R69WrdcMMNWrNmDckOaKZIeIAHatqZiYmJysvLU0pKirZv364HHnjANVkFQPNEwgM8UFPhffnll7ruuutUVlamuXPnaubMmbLb+esENGd8hwd4oFu3bvruu+9cPy9btkw///nP/RgRAHfxT1LATWVlZdqxY4datmypYcOGSZLmz5/PE86BAEHCA9z0xRdfyLIs9ejRQ4sWLVLXrl21YsUK/fa3v/V3aADcQMID3HTqDM02bdrozTffVKtWrTRjxgy9++67/g0OwFmR8AA3nX7BeWJioubMmSPLsnTHHXfo22+/9WN0AM6GhAe46fRr8CQpLS1NEyZMUElJiW6++WYdOXLET9EBOBtmaQJuqK6uVuvWrXXs2DHt379f7du3d22rqqrSlVdeqXXr1mnUqFF6/fXX/RgpgIZQ4QFu2LZtm44dO6aOHTvWSnaSFBoaqmXLlql9+/ZatmyZnn76aT9FCeBMSHiAG069w0p9HA6HXnvtNYWEhOjRRx/VqlWrmjI8AG6gpQkAMAIVHgDACCQ8AIARSHgAACOQ8AAARiDhAQCMQMIDABiBhAcAMAIJDwBgBBIeAMAIJDwAgBFIeAAAI5DwAABGIOEBAIzw//aD7mP54vc0AAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 500x500 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Diagonal affinity matrix, which results in an assortative structure\n",
        "w = np.eye(3)\n",
        "\n",
        "plt.matshow(w, aspect='auto', cmap='Blues')\n",
        "plt.gcf().set_size_inches(5, 5)\n",
        "plt.title(r'Affinity matrix $w$', fontsize=17)\n",
        "plt.xlabel(r'$K$', fontsize=15)\n",
        "plt.ylabel(r'$K$', fontsize=15)\n",
        "plt.xticks(ticks=[0,1,2], labels=[1,2,3], size=12)\n",
        "plt.tick_params(axis='x', bottom=True, top=False, labelbottom=True, labeltop=False)\n",
        "plt.colorbar()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "0e19c1f1",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAKFCAYAAACuprKBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABB9klEQVR4nO3de1yUZf7/8fcMCngCVAxEQS1NLI9hGtrBXIqspYPaum7lYc39WeqatI8tt/JQbdRWxrZRZushK1fT7eChbI3U/WbYAcNsNzHTAjVQS0FQQOH+/WHMesfceIMzzIy8nj7uR3EfL5hhPrzv676v22EYhiEAAHBGTl83AACAQEHRBADAJoomAAA2UTQBALCJogkAgE0UTQAAbKJoAgBgE0UTAACbKJoAANhE0QQAwCaKJgAANlE04deGDBkih8OhxYsX+7opAYOfGeA9FE1JxcXFeuaZZ3TttdeqY8eOCg0NVatWrRQfH69x48bpvffeE+Pa+5/09HTNnj1b3377ra+bgkaG917j1cTXDfC1FStW6K677tIPP/zgmhcWFqaTJ08qNzdXubm5evnll9W/f3+9+eab6tixow9b2/jExcWpe/fuCg8Pr7EsPT1d3333nYYMGaLOnTs3fOP8VG0/M3gG773Gq1EXzfnz52vSpEkyDEO9evXSgw8+qOTkZNeHTWFhodasWaOnnnpKn332mXbt2kXRbGBLlizxdRMCDj8zwHsabdHcunWrpk6dKsMwdMstt2jZsmUKDg42rRMVFaUJEyZo3LhxmjNnjpxOzmYDQKNmNFLDhg0zJBlxcXHG0aNHbW1TVVVl+vrkyZPGvHnzjMGDBxsRERFGaGio0a1bN2PatGnG/v373e6jU6dOhiRjw4YNxt69e40777zT6NChg9GsWTOjV69exqJFi1zrVlZWGs8995zRt29fo3nz5kZkZKQxbtw4o7CwsMH3LcmQZOzZs8ft8g0bNhiSjE6dOtXaroMHDxpTp0414uLijODgYCMuLs6YNm2acfjwYbf7veqqqwxJprYvWrTI1R5301VXXWXk5eUZTqfTkGR8/fXXbvdtGIbx1VdfGZKMoKAgY9++fZbr/VxeXp7x+OOPG9dcc41x/vnnGyEhIUZ4eLiRmJhoPPvss0Z5ebnb7crKyoynn37aGDBggBEWFmY0bdrUiI6ONi655BLj3nvvNbZv337W27j7mVUrKioy/vCHPxidO3c2QkJCjLi4OGPy5MnGwYMHXT/Xq666qsZ2Z/MaevN9WW3Dhg3GyJEjjZiYGKNp06ZG27Ztjeuuu85YtWqV5Tb1+Z7svPeq1ee1rk15ebnRtGnTWt/Tq1evNiQZvXv3rtO+YV+jLJp5eXmuN/lTTz1Vr32UlJQYQ4cOde0nODjYCAsLc33dunVr4+OPP66xXfUv6sKFC42oqChDkhEeHm44HA7Xtk888YRRVVVljBw50pBkhISEGM2aNXMtv/jii43jx4836L49UTSXLFlixMbGGpKMli1buj4AJBmXXHKJ20LjrgAsW7bMiIqKchXF1q1bG1FRUa7plltuMQzDMJKTkw1JxgMPPGDxKhrGfffdZ0gyhg0bZrmOOyNGjHC1PTQ01GjdurXpw3PIkCE1vp+Kigrj8ssvd63jdDqN1q1bu74PScZ999131ttYFc2DBw8aPXr0cG3XrFkzo0WLFoYko0uXLsZTTz11xqJZn9fQm+/LqqoqIzU11fSzP/33UJIxdepUt69hfb4nu++9+rxuZ5KTk2NIMlq1alXjD/hqDz/8sCHJuP322+u0b9jXKIvmK6+84nrj7tixo177mDhxoiHJaN68ubF48WKjoqLCMAzD+Pzzz41+/foZkoyYmJgaf6lW/6KGh4cbl19+ufHf//7XMAzDKC4uNqZOner6MHvggQeMsLAw4x//+IdRUVFhVFZWGmvXrnV9IPz1r3+t0SZv7tsTRTMiIsJISEgwPv30U8MwTv3lvGTJEiM0NNSQZPztb3+rsW1tqen0pODO66+/7jqbUFlZWWP5yZMnjZiYGEOSsXz5crf7sPLQQw8ZGRkZxjfffOP6ADt27Jjx+uuvGx06dDAkGY888ohpm8WLFxuSjHbt2hlr1qwxTpw4YRjGqQ/YnTt3Go8//rgxf/78s97G6md26623GpKMNm3aGG+99ZbrZ/LBBx8YsbGxRkRExBmLZn1eQ2++L59++mlDktGxY0djyZIlRnFxsWEYhnH06FFj/vz5rm2XLFnile/J6r1Xn9ftTF5++WVDkjF48GDLdYYPH25IMv7yl7/Uad+wr1EWzQceeMD116zVX2y12b17t+svxldffbXG8sLCQtdf8I899phpWfUvW5s2bYwjR46YllVVVRkXXnihq0C52/cjjzxyxg82b+zbE0UzPDzcOHjwYI3l06ZNszzu2RTN8vJyo23btoYkY/369TWWr1u3zpUWysrK3O6jPjZv3uz6ID/dXXfdZUgyHn/8cdv7qs827n5mubm5rtdw7dq1NbbZunWr6z1d23urPq+ht96Xhw8fNlq0aGG0aNHC+Oqrr9z+LJYvX25IMnr06OGV78nqvVef1+1MqhP13XffbbnO+eefb0gy1q1b57HjwqxRXtlSfXtJ69at5XA46rz9m2++qaqqKnXp0kW33XZbjeXnnXeeJk6cKElauXKl231MmjSpxi0BDodDQ4YMkSTFxsZq9OjRNbb7xS9+IUn6z3/+Y9k+b+77bEyaNEmRkZE15qekpHjluMHBwbrjjjskye2N/tXzRo8erZCQEI8dd9CgQYqIiNDevXu1b98+1/ywsDBJ0vfff297X/XZxp23335bktSjRw9df/31NZb369dPQ4cOPeN+zuY19PT7cuXKlSotLdUNN9yg+Ph4t8ccPny4QkJC9NVXX1n+DL3xvvTU63a6bdu2SZL69u3rdnlxcbH27NkjSerVq5fHjguzRlk0z9bWrVslyfXL7s7VV18tSfriiy9UWVlZY3nPnj3dbnfeeedJki666CK3V+tWLz9y5Ijlsb2577PRu3dvt/M7dOjgteP+9re/lXTqD53i4mLX/KKiIr311luSpPHjx9dr3x9++KFuv/12XXDBBWrevLkcDodrqv5eTv/QvO666yRJzz77rO644w69++67Onr0aK3HqM827lR/4A4ePNhyncsvv/yM+zmb19DT78usrCxJ0jvvvKPo6Gi3U8eOHXXixAlJUn5+vse/Jyueet1OV/0a9unTx+3ynJwcGYahyMhIxcTEnNWxYK1RFs22bdtKkg4fPlyvkX4OHTok6X+/VO506tRJknTy5Em3v3Tt27d3u11QUJAkKTo6utblJ0+etDy2N/d9Nqx+kUNDQ7123F69eunSSy/VsWPH9Prrr7vmL1u2TGVlZerZs6f69+9f5/0+9thjuuKKK/Taa69p9+7dOnnypNq0aaOoqChFRUW5PvxLS0td2wwZMkSzZs2S0+nUq6++quuvv14RERG65JJLNHv2bBUUFNQ4Tn22caf6PWv13pCs3xenO5vX0NPvy+o/SEpKSlRYWGg5VVVVSZKOHTvm8e/Jiqdet2r79+/XoUOHFBQUZJkic3JyJJEyva1RFs3qUznl5eXauXNnvfdTXl7uqSbBiyZMmCDJfIr25ZdfliSNGzeuzvvbvn27HnroIUnS3XffrdzcXJWVlemHH35QQUGBCgoKXB/EP/+jbPbs2crNzdWjjz6qa665Rs2aNdPnn3+uOXPmqFu3bvrggw9qHK8+2zQG1cXwgQcekHHq+oxap9rODHmDJ1+36pTZrVs3NWvWzO062dnZkiia3tYoi+ZVV13l+v81a9bUeft27dpJkvLy8izX+e677yRJTZo0UURERJ2P4W+q/9ovKytzu7yoqKghm1Mno0ePVvPmzbV582bt2rVLO3fuVFZWlpo0aaLbb7+9zvt74403VFVVpauuukoZGRm68MILTacVKysrTcMy/twFF1ygBx54QP/61790+PBhvfvuu+rTp49KSko0ZswYt6fz67PN6ar77GrrY6tr+vG1qKgoSbX/Hvra2b5u1bZv3y7J+tTsyZMntXbtWknWp5vhGY2yaMbGxmrYsGGSTvU5lJSU2NquOjX069dPkrR582bLUzcbNmyQdOoNXF1wAll14T/9wpbTffbZZw3YmlOqC9WZTrGHhYVp5MiRkk6lzerEOWzYMNcHb11U/wysTut+/PHHOn78uK19NW3aVNddd53++c9/uvb9zTffeHyb6g/bjz76yHKdzZs322qzv7jsssskSe+//77XuhSs2H3vna4+r1u16vW6dOnidvk///lP1x9qFE3vapRFU5IeeeQRNW3aVHl5eRo7dqwqKios162srNTMmTP1f//3f5JOXZHndDq1d+9evfrqqzXWP3DggF566SVJcn1YB7rqUz6rVq2qsezw4cP6+9//3tBNcl2haOdCjeoLgpYsWaJXXnlFUv0vAKo+bm5ubo1lVVVVmj17ttvtanuPnX7K7fQ0X59t3LnxxhslSf/973+1bt26Gsu3bdumzMzMWvfhb2699Va1aNFC33//vZ544ola1z18+LBHj32m956nXrdq1X8UHDx4sMayvXv3avLkyZJOFfOLL77Y1j5RP422aCYkJCg9PV3SqdNtl156qVasWGG6wrKwsFALFy5Uz5499cgjj7j6UDp37uzqJ5s6dapeffVV1xV6OTk5uu6661RaWqqYmBjdddddDfuNeUl18Z8/f75eeeUV14fCZ599pmuuucYn/bvVHw7/+Mc/zvjhc+WVV6pr167Kz8/X3r17FRkZqV/+8pf1Om71LRBr1qzR3LlzXd/7t99+q1/96lf697//rRYtWtTYbsyYMZowYYLWr19vOruRm5vrKuAdOnQwfejVZxt3unfvrhEjRkiSbr/9dq1evdr1ft60aZNuuukmVyEIFJGRkXr00UclSQ8++KB+//vfu265kE5dILR+/XrdcccduvXWWz167DO99zz1ulXr3r27JGnp0qV65513dPLkSR0/flyvv/66EhMT1bRpU0lyXckNL/LBvaF+ZenSpUabNm1Mw26Fh4cbzZs3N80bNGiQaWzSkpIS4+qrr3YtDwkJqTGM3pYtW2oc70w3Rc+aNcuQZIwdO9bt8j179riO0ZD7Li8vN/r37+9a3rRpU6Nly5auEXeqR1k609izdT1ubYMbZGZmurYLDg42OnbsaHTq1MkYNWqU2+M89thjrvWnTZvmdh07qqqqXGMX66ch0qpH03E6ncaLL77o9nu+6aabXNs4HA6jdevWpqHimjdvbrz//vumY9VnG6uf2YEDB4zu3bu7tj19GL2uXbsaTz75pCHJuPbaa2t8z2fzGnrzfWkYhjFz5kzTkHytWrUyIiIiTPOGDBni0e/pTO+9+rxutSksLHS9x6RTYyVXf3833HCDMWfOHEOSMWLECNv7RP002qRZbfTo0dq9e7eefvppJSUlqX379jp+/LgcDofi4+M1fvx4vf/++9q8ebPp0vQWLVroX//6l1544QUlJiYqJCRE5eXl6tq1q37/+9/rP//5jwYOHOjD78yzgoOD9f7772v69OmKi4uTw+FQRESE7r77bmVnZ/vkkWlDhw7Vm2++qauuukrNmjXTvn379N1331le0HLLLbe4/r++p2alUzfkv/XWW5o1a5a6du2qoKAgNWnSRMOGDdP69ev1u9/9zu12aWlpeuKJJ3TttdeqS5cuKisrk2EYuvDCC3X33Xdr+/btrhR7NttYadeunT7++GPde++96tSpkyorK9W2bVtNnTpVH3/8sSutBNqFa3PmzNHWrVv129/+Vueff75OnjypY8eOqWPHjvrlL3+pv/3tb1qxYoVHj3mm954nXzfp1L2qH374oZKTk9WmTRu1atVKV155pV577TWtWbNGu3btkkR/ZkNwGEY9blQEAtC8efN01113qW/fvvr888993Ry/M3bsWC1ZskQzZ87UnDlzfN0cwC81+qSJxmP+/PmSpDvvvNPHLfE/eXl5rqs6k5KSfNwawH9RNNEoPPfcc/r8888VERHhGo+2sdmzZ48mTpyoLVu2uG6JOXHihN59910NHTpUpaWl6t+/v6644goftxTwX5yexTlr7969uvzyy3X06FH9+OOPkqT09HRNmzbNxy3zjR07dqhHjx6ur1u3bq2SkhLXld8xMTHKzMy0HPwcAEkT57CTJ0/qu+++U1FRkS644AI988wzjbZgSqcG9XjyySf1i1/8QnFxcTp+/LhCQkLUq1cvzZgxQ9u2baNgImD8+9//VkpKimJiYlwX553Jxo0bdckllygkJERdu3Z1+/SjMyFpAgACzrvvvqvNmzcrISFBw4cP15tvvqmbb77Zcv09e/aoZ8+emjRpku68805lZmbqnnvu0dq1a5WcnGz7uBRNAEBAczgcZyya9913n9auXasvv/zSNe/Xv/61jhw54naULCtNzqahgayqqkr79+9Xq1at6vUgagDwN4Zh6OjRo4qJiXH7bNK6Kisrq3VIQE8yDKPGZ3FISIjHHhCflZVV48rw5ORk3XPPPXXaT0AXzYyMDD355JMqKChQnz599Le//U0DBgywte3+/fsVGxvr5RYCQMPLz88/6wFHysrK1KxVW+mk++eQelrLli1rPDxj1qxZlmM511VBQUGNBzRERUWpuLhYx48ft3zk2s8FbNFcvny5UlNTNW/ePA0cOFDp6elKTk5Wbm6u60nvtWnVqpUkKfiisXIEBXu7uaijvI1P+boJQMA5Wlysrl1iXZ9vZ6OiokI6eUwhF42VvP0ZWVmhkv++rPz8fNMYyJ5KmZ4UsEVz7ty5mjhxoms4tHnz5mnt2rVauHCh7r///jNuX30awBEUTNH0Q4E2eDjgTzza5dQk1OufkYbj1KnksLAwr/3uR0dHq7Cw0DSvsLBQYWFhtlOmFKC3nFRUVCg7O9t0ftrpdCopKUlZWVk+bBkAwB8lJibWePzd+vXrlZiYWKf9BGTRPHTokCorK92en7YarLu8vFzFxcWmCQBwBg5JDoeXp7o3q6SkRDk5OcrJyZF06paSnJwc5eXlSZJmzJihMWPGuNafNGmSdu/erT/+8Y/asWOHnn/+eb3++uuaPn16nY4bkEWzPtLS0hQeHu6auAgIAALXZ599pn79+qlfv36SpNTUVPXr108zZ86UJH3//feuAipJXbp00dq1a7V+/Xr16dNHTz/9tP7+97/X6R5NKUD7NCMjIxUUFOT2/HR0dLTbbWbMmKHU1FTX18XFxRROAAhQQ4YMUW3DDLgb7WfIkCFn/YSjgEyawcHBSkhIMJ2frqqqUmZmpuX56ZCQEFcnszc7mwHgnOJwNswUIAIyaUqnovjYsWPVv39/DRgwQOnp6SotLT2rhwsDAFCbgC2ao0aN0sGDBzVz5kwVFBSob9++WrduXY2LgwAAZ6H6Yh1vHyNABGzRlKQpU6ZoypQpvm4GAKCRCOiiCQDwsobocwygPs3AaSkAAD5G0gQAWKNP04SkCQCATSRNAEAtGuI+ysDJb4HTUgAAfIykCQCwRp+mCUkTAACbSJoAAGvcp2kSOC0FAMDHSJoAAGv0aZqQNAEAsImiCQCATZyeBQBY40Igk8BpKQAAPkbSBABY40IgE5ImAAA2kTQBANbo0zQJnJYCAOBjJE0AgDWHowGSJn2aAACcc0iaAABrTsepydvHCBAkTQAAbCJpAgCscfWsSeC0FAAAHyNpAgCsMSKQCUkTAACbSJoAAGv0aZoETksBAPAxkiYAwBp9miYkTQAAbKJoAgBgE6dnAQDWuBDIJHBaCgCAj5E0AQDWuBDIhKQJAIBNJE0AgDX6NE0Cp6UAAPgYSRMAYI0+TROSJgAANpE0AQC1aIA+zQDKb4HTUgAAfIykCQCwRp+mCUkTAACbSJoAAGsORwPcp0nSBADgnEPSBABYY0Qgk8BpKQAAPkbRBADAJk7PAgCsccuJCUkTAACbSJoAAGtcCGQSOC0FAMDHSJoAAGv0aZqQNAEAsImkCQCwRp+mSeC0FAAAHyNpAgCs0adpQtIEAMAmkiYAwJLD4ZCDpOlC0gQAwCaSJgDAEknTjKQJAIBNJE0AgDXHT5O3jxEgSJoAANhE0QQAwCZOzwIALHEhkBlJEwAAm0iaAABLJE0zkiYAADaRNAEAlkiaZiRNAABsImkCACyRNM1ImgAA2ETSBABYYxg9E5ImAAA2kTQBAJbo0zQjaQIAYBNJEwBgyeFQAyRN7+7ek0iaAADY5JdF89///rdSUlIUExMjh8Oht956y7TcMAzNnDlT7du3V7NmzZSUlKSvv/7aN40FgHOYQw5Xv6bXpgCKmn5ZNEtLS9WnTx9lZGS4Xf6Xv/xFzz77rObNm6ePP/5YLVq0UHJyssrKyhq4pQCAxsQv+zSHDRumYcOGuV1mGIbS09P14IMP6qabbpIkLVmyRFFRUXrrrbf061//uiGbCgDnNK6eNfPLpFmbPXv2qKCgQElJSa554eHhGjhwoLKysiy3Ky8vV3FxsWkCAKAuAq5oFhQUSJKioqJM86OiolzL3ElLS1N4eLhrio2N9Wo7AQDnnoArmvU1Y8YMFRUVuab8/HxfNwkA/J+jgaYAEXBFMzo6WpJUWFhoml9YWOha5k5ISIjCwsJMEwAAdRFwRbNLly6Kjo5WZmama15xcbE+/vhjJSYm+rBlAHAO8vbtJqdGT/D1d2mbX149W1JSol27drm+3rNnj3JyctSmTRvFxcXpnnvu0aOPPqpu3bqpS5cueuihhxQTE6Obb77Zd40GAJzz/LJofvbZZ7r66qtdX6empkqSxo4dq8WLF+uPf/yjSktL9bvf/U5HjhzR5ZdfrnXr1ik0NNRXTQaAc1JD3HLi9VtaPMgvi+aQIUNkGIblcofDoYcfflgPP/xwA7YKANDY+WXRBAD4B5KmWcBdCAQAgK9QNAEA1vz4Ps2MjAx17txZoaGhGjhwoD755JNa109PT1f37t3VrFkzxcbGavr06XUes5yiCQAIOMuXL1dqaqpmzZqlrVu3qk+fPkpOTtaBAwfcrr906VLdf//9mjVrlr766istWLBAy5cv15/+9Kc6HZeiCQCw5PXHgtWzz3Tu3LmaOHGixo8fr4suukjz5s1T8+bNtXDhQrfrf/TRRxo8eLB+85vfqHPnzrr22ms1evToM6bTn6NoAgACSkVFhbKzs00P7nA6nUpKSrJ8cMegQYOUnZ3tKpK7d+/WO++8o+uvv75Ox+bqWQCApYa8evbnT58KCQlRSEhIjfUPHTqkyspKtw/u2LFjh9tj/OY3v9GhQ4d0+eWXyzAMnTx5UpMmTeL0LAAgMMXGxpqeRpWWluaxfW/cuFGPPfaYnn/+eW3dulVvvPGG1q5dq0ceeaRO+yFpAgAsNWTSzM/PNz1Mw13KlKTIyEgFBQXV6cEdDz30kO644w7deeedkqRevXq5RpZ74IEH5HTay5AkTQCAX/j5k6isimZwcLASEhJMD+6oqqpSZmam5YM7jh07VqMwBgUFSVKtI9D9HEkTABBwUlNTNXbsWPXv318DBgxQenq6SktLNX78eEnSmDFj1KFDB9cp3pSUFM2dO1f9+vXTwIEDtWvXLj300ENKSUlxFU87KJoAAEv+OozeqFGjdPDgQc2cOVMFBQXq27ev1q1b57o4KC8vz5QsH3zwQTkcDj344IPat2+f2rVrp5SUFP35z3+uW1uNuuTSc0hxcbHCw8MV0muiHEHBvm4Ofubwp8/5uglAwCkuLlZU23AVFRWZ+gbru6/w8HBFjXtFzuDmHmqhe1UVx1S4+A6PtNvbSJoAAGtnMcxdnY4RILgQCAAAm0iaAABL/tqn6SskTQAAbCJpAgAskTTNSJoAANhE0gQAWCJpmpE0AQCwiaQJALDGfZomJE0AAGwiaQIALNGnaUbSBADAJpImAMASSdOMpAkAgE0kTQCAJYcaIGkG0OWzJE0AAGyiaAIAYBOnZwEAlrgQyIykCQCATSRNAIA1htEzIWkCAGATSRMAYIk+TTOSJgAANpE0AQCWSJpmJE0AAGwiaQIALDkcpyZvHyNQkDQBALCJpAkAsHQqaXq7T9Oru/cokiYAADaRNAEA1hqgT5MRgQAAOAeRNAEAlrhP04ykCQCATRRNAABs4vQsAMASgxuYkTQBALCJpAkAsOR0OuR0ejcKGl7evyeRNAEAsImkCQCwRJ+mGUkTAACbSJoAAEsMbmBG0gQAwCaSJgDAEn2aZiRNAABsImkCACzRp2lG0gQAwCaSJgDAEknTjKQJAIBNJE0AgCWunjUjaQIAYBNFEwAAmzg9CwCw5FADXAikwDk/S9IEAMAmkiYAwBIXApmRNAEAsImkCQCwxOAGZiRNAABsImkCACzRp2lG0gQAwCaSJgDAEn2aZiRNAABsImkCACzRp2lG0gQAwCaSJgDAEn2aZiRNAABsImkCAKw1QJ9mAD3kxD+LZlpamt544w3t2LFDzZo106BBg/TEE0+oe/furnXKysp07733atmyZSovL1dycrKef/55RUVF1elYeRufUlhYmKe/BZyl1pdO8XUTUIvDnz7n6yYAPuGXp2c3bdqkyZMna8uWLVq/fr1OnDiha6+9VqWlpa51pk+frtWrV2vFihXatGmT9u/fr+HDh/uw1QBw7qnu0/T2FCj8MmmuW7fO9PXixYt13nnnKTs7W1deeaWKioq0YMECLV26VEOHDpUkLVq0SD169NCWLVt02WWX+aLZAIBznF8mzZ8rKiqSJLVp00aSlJ2drRMnTigpKcm1Tnx8vOLi4pSVleWTNgIAzn1+mTRPV1VVpXvuuUeDBw9Wz549JUkFBQUKDg5WRESEad2oqCgVFBS43U95ebnKy8tdXxcXF3utzQBwrmBwAzO/T5qTJ0/Wl19+qWXLlp3VftLS0hQeHu6aYmNjPdRCAEBj4ddFc8qUKVqzZo02bNigjh07uuZHR0eroqJCR44cMa1fWFio6Ohot/uaMWOGioqKXFN+fr43mw4A5wQuBDLzy6JpGIamTJmiN998Ux988IG6dOliWp6QkKCmTZsqMzPTNS83N1d5eXlKTEx0u8+QkBCFhYWZJgAA6sIv+zQnT56spUuX6u2331arVq1c/ZTh4eFq1qyZwsPDNWHCBKWmpqpNmzYKCwvT1KlTlZiYyJWzAOBB9Gma+WXRfOGFFyRJQ4YMMc1ftGiRxo0bJ0l65pln5HQ6NWLECNPgBgAAeItfFk3DMM64TmhoqDIyMpSRkdEALQKAxokB2838sk8TAAB/5JdJEwDgH0iaZiRNAABsImkCACxx9awZSRMAAJtImgAAS/RpmpE0AQCwiaQJALBEn6YZSRMAAJsomgAA2MTpWQCAJS4EMiNpAgBgE0kTAGDJoQa4EMi7u/cokiYAADaRNAEAlpwOh5xejpre3r8nkTQBALCJpAkAsMTgBmYkTQAAbCJpAgAscZ+mGUkTAACbSJoAAEtOx6nJ28cIFCRNAEBAysjIUOfOnRUaGqqBAwfqk08+qXX9I0eOaPLkyWrfvr1CQkJ04YUX6p133qnTMUmaAABrjgboc6zH7pcvX67U1FTNmzdPAwcOVHp6upKTk5Wbm6vzzjuvxvoVFRW65pprdN5552nlypXq0KGDvvvuO0VERNTpuBRNAEDAmTt3riZOnKjx48dLkubNm6e1a9dq4cKFuv/++2usv3DhQv3444/66KOP1LRpU0lS586d63xcTs8CACxV36fp7UmSiouLTVN5ebnbNlVUVCg7O1tJSUmueU6nU0lJScrKynK7zapVq5SYmKjJkycrKipKPXv21GOPPabKyso6/TwomgAAvxAbG6vw8HDXlJaW5na9Q4cOqbKyUlFRUab5UVFRKigocLvN7t27tXLlSlVWVuqdd97RQw89pKefflqPPvpondrI6VkAgCXHT/+8fQxJys/PV1hYmGt+SEiIx45RVVWl8847T/Pnz1dQUJASEhK0b98+Pfnkk5o1a5bt/VA0AQB+ISwszFQ0rURGRiooKEiFhYWm+YWFhYqOjna7Tfv27dW0aVMFBQW55vXo0UMFBQWqqKhQcHCwrTZyehYAEFCCg4OVkJCgzMxM17yqqiplZmYqMTHR7TaDBw/Wrl27VFVV5Zq3c+dOtW/f3nbBlCiaAIBaVA9u4O2prlJTU/XSSy/p5Zdf1ldffaW77rpLpaWlrqtpx4wZoxkzZrjWv+uuu/Tjjz9q2rRp2rlzp9auXavHHntMkydPrtNxOT0LAAg4o0aN0sGDBzVz5kwVFBSob9++WrdunevioLy8PDmd/8uFsbGxeu+99zR9+nT17t1bHTp00LRp03TffffV6bgUTQCAJX8esH3KlCmaMmWK22UbN26sMS8xMVFbtmyp17GqcXoWAACbSJoAAEs8hNqMpAkAgE0kTQCAJafDIaeXo6C39+9JJE0AAGwiaQIALNGnaUbSBADAJpImAMCSP9+n6QskTQAAbCJpAgAs0adpRtIEAMAmkiYAwBL3aZqRNAEAsImiCQCATZyeBQBYcvw0efsYgYKkCQCATSRNAIAlBjcwI2kCAGATSRMAYMnpODV5+xiBgqQJAIBNJE0AgCX6NM1ImgAA2ETSBADUKoCCoNeRNAEAsImkCQCwRJ+mGUkTAACbSJoAAEvcp2lG0gQAwCaSJgDAEn2aZiRNAABsomgCAGATp2cBAJZ4CLUZSRMAAJtImgAAS06HQ04vX6jj7f17EkkTAACbSJoAAEsOh/cHbA+goEnSBADALpImAMASgxuYkTQBALCJpAkAsESfphlJEwAAm2wlzVatWql3797q3bu3+vbtqz59+qh3795q3ry5t9sHAPAh7tM0s1U0S0tLlZWVpaysLFeHrcPh0AUXXOAqotX/7dChg1cbDACAr9gqmsXFxdq2bZu2bdumnJwcbdu2TV9++aW+/vprff3111q5cqVr3TZt2piKaN++fdWrVy+vfQMAAO+hT9PMVtFs2bKlBg8erMGDB7vmVVVVaefOna4iWv3fgoICffDBB9qwYYOkU4n05MmT3mk9AAANqN5XzzqdTsXHxys+Pl6//vWvXfMPHjyot956S2lpafr222890UYAgI9wn6aZx2452bp1q95++22tWrVKX3zxhQzDkCTFx8d76hAAAPhUvYvmiRMn9MEHH+jtt9/WmjVrtG/fPhmGoSZNmujKK69USkqKbrzxRnXt2tWT7QUANCCnvH9vYiDd+1inth4+fFivvvqqfvWrXykyMlLXX3+95s2bp+LiYt1666165ZVXVFhYqA0bNig1NbXeBfOFF15Q7969FRYWprCwMCUmJurdd991LS8rK9PkyZPVtm1btWzZUiNGjFBhYWG9jgUAgF22kuYzzzyjVatWafPmza6Lejp37qxx48YpJSVFQ4YMUZMmnhtcqGPHjnr88cfVrVs3GYahl19+WTfddJM+//xzXXzxxZo+fbrWrl2rFStWKDw8XFOmTNHw4cO1efNmj7UBAICfs1Xp7r33XjkcDnXv3l233367UlJSvHobSUpKiunrP//5z3rhhRe0ZcsWdezYUQsWLNDSpUs1dOhQSdKiRYvUo0cPbdmyRZdddpnX2gUAjQ0XApnVKR7m5uZq/vz5+uSTT0z3YXbp0sVb7VNlZaVWrFih0tJSJSYmKjs7WydOnFBSUpJrnfj4eMXFxSkrK8uyaJaXl6u8vNz1dXFxsdfaDAA4N9kqmr///e/1xRdfaNu2bcrLy1NeXp5WrVrl+uugepi90wc16NWrl0JCQurdsO3btysxMVFlZWVq2bKl3nzzTV100UXKyclRcHCwIiIiTOtHRUWpoKDAcn9paWmaM2dOvdsDAI2RwyE5GdzAxVbRTE9Pd/1/fn5+jQENdu/erQ8//FAffvihq5AGBQWpW7du6tevn1599dU6N6x79+7KyclRUVGRVq5cqbFjx2rTpk113k+1GTNmKDU11fV1cXGxYmNj670/AEDjU+erd2JjYxUbG2vqdywpKdEXX3xhKqZffvmlvvrqK+3YsaNeRTM4ONh19W1CQoI+/fRT/fWvf9WoUaNUUVGhI0eOmNJmYWGhoqOjLfcXEhJyVskXABojZwMkTW/v35M8cslry5YtNWjQIA0aNMg1r3qYvW3btnniEKqqqlJ5ebkSEhLUtGlTZWZmasSIEZJO9bXm5eUpMTHRI8cCAMAdrz2E+vRh9upqxowZGjZsmOLi4nT06FEtXbpUGzdu1Hvvvafw8HBNmDBBqampatOmjcLCwjR16lQlJiZy5SwAeBhXz5p5rWiejQMHDmjMmDH6/vvvFR4ert69e+u9997TNddcI+nUfaNOp1MjRoxQeXm5kpOT9fzzz/u41QCAc51fFs0FCxbUujw0NFQZGRnKyMhooBYBQONEn6ZZIA35BwCAT/ll0gQA+AceQm1G0gQAwCaSJgDAktPhkNPLUdDb+/ckkiYAADaRNAEAlngItVkgtRUAAJ+iaAIAYBOnZwEAlrjlxIykCQCATSRNAIAlpxrglhMFTtQkaQIAYBNJEwBgiT5NM5ImAAA2kTQBAJZ4NJgZSRMAAJtImgAASw6H9wdUp08TAIBzEEkTAGCJq2fNSJoAANhE0gQAWOLqWTOSJgAANpE0AQCWHD/98/YxAgVJEwAAm0iaAABL9GmakTQBALCJogkAgE2cngUAWOL0rBlJEwAAm0iaAABLDodDDq8P2B44UZOkCQCATSRNAIAl+jTNSJoAANhE0gQAWOLRYGYkTQAAbCJpAgAsOR0OOb0cBb29f08iaQIAYBNJEwBgiatnzUiaAICAlJGRoc6dOys0NFQDBw7UJ598Ymu7ZcuWyeFw6Oabb67zMSmaAABrjv9dQeutqT7PoF6+fLlSU1M1a9Ysbd26VX369FFycrIOHDhQ63bffvut/vCHP+iKK66o14+DogkACDhz587VxIkTNX78eF100UWaN2+emjdvroULF1puU1lZqdtuu01z5szR+eefX6/jUjQBAJaccjTIVBcVFRXKzs5WUlLS/9rpdCopKUlZWVmW2z388MM677zzNGHChHr/PLgQCADgF4qLi01fh4SEKCQkpMZ6hw4dUmVlpaKiokzzo6KitGPHDrf7/vDDD7VgwQLl5OScVRtJmgAAvxAbG6vw8HDXlJaW5pH9Hj16VHfccYdeeuklRUZGntW+SJoAAEsNOYxefn6+wsLCXPPdpUxJioyMVFBQkAoLC03zCwsLFR0dXWP9b775Rt9++61SUlJc86qqqiRJTZo0UW5uri644AJbbSVpAgD8QlhYmGmyKprBwcFKSEhQZmama15VVZUyMzOVmJhYY/34+Hht375dOTk5runGG2/U1VdfrZycHMXGxtpuI0kTAGDJXwc3SE1N1dixY9W/f38NGDBA6enpKi0t1fjx4yVJY8aMUYcOHZSWlqbQ0FD17NnTtH1ERIQk1Zh/JhRNAEDAGTVqlA4ePKiZM2eqoKBAffv21bp161wXB+Xl5cnp9PzJVIomAMCSPw/YPmXKFE2ZMsXtso0bN9a67eLFi+t1TPo0AQCwiaQJALDEQ6jNSJoAANhE0gQAWHKqAfo06zNiu4+QNAEAsImkCQCwRJ+mGUkTAACbSJoAAEtOeT9dBVJ6C6S2AgDgUyRNAIAlh8Mhh5c7Hb29f08iaQIAYBNFEwAAmzg9CwCw5Php8vYxAgVJEwAAm0iaAABL/vxoMF8gaQIAYBNJEwBQq8DJgd5H0gQAwCaSJgDAEgO2m5E0AQCwiaQJALDEMHpmJE0AAGwiaQIALPFoMLNAaisAAD5F0gQAWKJP04ykCQCATSRNAIAlnnJiRtIEAMAmkiYAwBJ9mmYkTQAAbPL7pPn4449rxowZmjZtmtLT0yVJZWVluvfee7Vs2TKVl5crOTlZzz//vKKionzbWHjM4U+f83UTUIvWl07xdRPghlFZ4esmnPP8Oml++umnevHFF9W7d2/T/OnTp2v16tVasWKFNm3apP3792v48OE+aiUAnLucDTQFCr9ta0lJiW677Ta99NJLat26tWt+UVGRFixYoLlz52ro0KFKSEjQokWL9NFHH2nLli0+bDEA4Fznt0Vz8uTJuuGGG5SUlGSan52drRMnTpjmx8fHKy4uTllZWZb7Ky8vV3FxsWkCANSu+kIgb0+Bwi/7NJctW6atW7fq008/rbGsoKBAwcHBioiIMM2PiopSQUGB5T7T0tI0Z84cTzcVANCI+F3SzM/P17Rp0/Taa68pNDTUY/udMWOGioqKXFN+fr7H9g0A5ypHA02Bwu+KZnZ2tg4cOKBLLrlETZo0UZMmTbRp0yY9++yzatKkiaKiolRRUaEjR46YtissLFR0dLTlfkNCQhQWFmaaAACoC787PfuLX/xC27dvN80bP3684uPjdd999yk2NlZNmzZVZmamRowYIUnKzc1VXl6eEhMTfdFkADhnORynJm8fI1D4XdFs1aqVevbsaZrXokULtW3b1jV/woQJSk1NVZs2bRQWFqapU6cqMTFRl112mS+aDABoJPyuaNrxzDPPyOl0asSIEabBDQAAnuWUQ04v9zp6e/+eFBBFc+PGjaavQ0NDlZGRoYyMDN80CADQKAVE0QQA+AZ9mmZ+d/UsAAD+iqQJALDk+Omft48RKEiaAADYRNIEAFiiT9OMpAkAgE0UTQAAbOL0LADAkqMBBjfgQiAAAM5BJE0AgCUuBDIjaQIAYBNJEwBgiaRpRtIEAMAmkiYAwBLD6JmRNAEAsImkCQCw5HScmrx9jEBB0gQAwCaSJgDAEn2aZiRNAABsImkCACxxn6YZSRMAAJtImgAASw55v88xgIImSRMAALsomgAA2MTpWQCAJQY3MCNpAgBgE0kTAGCJwQ3MSJoAANhE0gQAWGJwAzOSJgAANpE0AQCWHPL+4AMBFDRJmgAA2EXSBABYcsohp5c7HZ0BlDVJmgAA2ETSBABYok/TjKQJAIBNJE0AgDWipglJEwAAm0iaAABLjD1rRtIEAMAmkiYAwFoDjD0bQEGTpAkAgF0UTQAAbOL0LADAEnecmJE0AQCwiaQJALBG1DQhaQIAYBNJEwBgicENzEiaAADYRNIEAFhyNMDgBl4fPMGDSJoAANhE0gQAWOLiWTOSJgAANpE0AQDWiJomJE0AAGwiaQIALHGfphlJEwAAmyiaAABL1fdpenuqj4yMDHXu3FmhoaEaOHCgPvnkE8t1X3rpJV1xxRVq3bq1WrduraSkpFrXt0LRBAAEnOXLlys1NVWzZs3S1q1b1adPHyUnJ+vAgQNu19+4caNGjx6tDRs2KCsrS7Gxsbr22mu1b9++Oh2XogkACDhz587VxIkTNX78eF100UWaN2+emjdvroULF7pd/7XXXtPdd9+tvn37Kj4+Xn//+99VVVWlzMzMOh2XogkAsORooEmSiouLTVN5ebnbNlVUVCg7O1tJSUmueU6nU0lJScrKyrL1fR07dkwnTpxQmzZtbP4kfjpOndYGAMBLYmNjFR4e7prS0tLcrnfo0CFVVlYqKirKND8qKkoFBQW2jnXfffcpJibGVHjt4JYTAIC1BhzcID8/X2FhYa7ZISEhXjnc448/rmXLlmnjxo0KDQ2t07YUTQCAXwgLCzMVTSuRkZEKCgpSYWGhaX5hYaGio6Nr3fapp57S448/rvfff1+9e/eucxs5PQsAsORooH91ERwcrISEBNNFPNUX9SQmJlpu95e//EWPPPKI1q1bp/79+9fr50HSBAAEnNTUVI0dO1b9+/fXgAEDlJ6ertLSUo0fP16SNGbMGHXo0MHVL/rEE09o5syZWrp0qTp37uzq+2zZsqVatmxp+7gUTQCAJX99CPWoUaN08OBBzZw5UwUFBerbt6/WrVvnujgoLy9PTuf/Tqa+8MILqqio0MiRI037mTVrlmbPnm37uBRNAEBAmjJliqZMmeJ22caNG01ff/vttx45JkUTAGCJJ4OZcSEQAAA2kTQBANaImiYkTQAAbCJpAgAs8RBqM5ImAAA2kTQBAJb89T5NXyFpAgBgE0kTAGCJi2fNSJoAANhE0QQAwCa/LJqzZ8+Ww+EwTfHx8a7lZWVlmjx5stq2bauWLVtqxIgRNZ6rBgDwAEcDTQHCL4umJF188cX6/vvvXdOHH37oWjZ9+nStXr1aK1as0KZNm7R//34NHz7ch60FADQGfnshUJMmTdw+gbuoqEgLFizQ0qVLNXToUEnSokWL1KNHD23ZskWXXXZZQzcVAM5ZDG5g5rdJ8+uvv1ZMTIzOP/983XbbbcrLy5MkZWdn68SJE0pKSnKtGx8fr7i4OGVlZfmquQCARsAvk+bAgQO1ePFide/eXd9//73mzJmjK664Ql9++aUKCgoUHBysiIgI0zZRUVGuJ3G7U15ervLyctfXxcXF3mo+AJwzGNzAzC+L5rBhw1z/37t3bw0cOFCdOnXS66+/rmbNmtVrn2lpaZozZ46nmggAaIT89vTs6SIiInThhRdq165dio6OVkVFhY4cOWJap7Cw0G0faLUZM2aoqKjINeXn53u51QAQ+Lh41iwgimZJSYm++eYbtW/fXgkJCWratKkyMzNdy3Nzc5WXl6fExETLfYSEhCgsLMw0AQBQF355evYPf/iDUlJS1KlTJ+3fv1+zZs1SUFCQRo8erfDwcE2YMEGpqalq06aNwsLCNHXqVCUmJnLlLAB4GuPomfhl0dy7d69Gjx6tH374Qe3atdPll1+uLVu2qF27dpKkZ555Rk6nUyNGjFB5ebmSk5P1/PPP+7jVAIBznV8WzWXLltW6PDQ0VBkZGcrIyGigFgFA48R9mmYB0acJAIA/8MukCQDwEw1wn2YABU2SJgAAdpE0AQCWuHjWjKQJAIBNFE0AAGzi9CwAwBrnZ01ImgAA2ETSBABYYnADM5ImAAA2kTQBAJZ4CLUZSRMAAJtImgAAS1w8a0bSBADAJpImAMAaUdOEpAkAgE0kTQCAJe7TNCNpAgBgE0kTAGDJoQa4T9O7u/cokiYAADaRNAEAlrh41oykCQCATRRNAABs4vQsAMASA7abkTQBALCJpAkAqAWXAp2OpAkAgE0kTQCAJfo0zUiaAADYRNIEAFiiR9OMpAkAgE0kTQCAJfo0zUiaAADYRNIEAFjiIdRmJE0AAGwiaQIArHH5rAlJEwAAm0iaAABLBE0zkiYAADaRNAEAlrhP04ykCQCATRRNAABs4vQsAMASgxuYkTQBALCJpAkAsMY9JyYkTQAAbCJpAgAsETTNSJoAANhE0gQAWGJwAzOSJgAANpE0AQC18P59moHUq0nSBADAJpImAMASfZpmJE0AAGyiaAIAYBNFEwAAm+jTBABYok/TjKQJAIBNFE0AAGzi9CwAwBIPoTYjaQIAYBNJEwBgiQuBzEiaAADYRNIEAFjiIdRmJE0AAGwiaQIArBE1TUiaAADYRNIEAFjiPk0zkiYAADaRNAEAlrhP04ykCQCATSRNAIAlLp41I2kCAGATSRMAYI2oaULSBADAJoomAMCSo4H+1UdGRoY6d+6s0NBQDRw4UJ988kmt669YsULx8fEKDQ1Vr1699M4779T5mBRNAEDAWb58uVJTUzVr1ixt3bpVffr0UXJysg4cOOB2/Y8++kijR4/WhAkT9Pnnn+vmm2/WzTffrC+//LJOx3UYhmF44hsINMXFxQoPD1fhD0UKCwvzdXOAgNL60im+bgLcMCorVL79JRUVnf3nWvVnZMEh739GFhcXKzoyvE7tHjhwoC699FI999xzkqSqqirFxsZq6tSpuv/++2usP2rUKJWWlmrNmjWueZdddpn69u2refPm2W4rSRMAYKl6cANvT3VRUVGh7OxsJSUlueY5nU4lJSUpKyvL7TZZWVmm9SUpOTnZcn0rjfbq2eqAfbS42MctAQKPUVnh6ybAjerXxZMnEIsb4DOy+hg/P1ZISIhCQkJqrH/o0CFVVlYqKirKND8qKko7duxwe4yCggK36xcUFNSprY22aB49elSS1LVLrI9bAgCedfToUYWHh5/VPoKDgxUdHa1uDfQZ2bJlS8XGmo81a9YszZ49u0GOb1ejLZoxMTHKz89Xq1at5AikgQ8tFBcXKzY2Vvn5+fTR+hleG/91rr02hmHo6NGjiomJOet9hYaGas+ePaqoaJizCoZh1PgsdpcyJSkyMlJBQUEqLCw0zS8sLFR0dLTbbaKjo+u0vpVGWzSdTqc6duzo62Z4XFhY2Dnxy38u4rXxX+fSa3O2CfN0oaGhCg0N9dj+PCU4OFgJCQnKzMzUzTffLOnUhUCZmZmaMsX9RWqJiYnKzMzUPffc45q3fv16JSYm1unYjbZoAgACV2pqqsaOHav+/ftrwIABSk9PV2lpqcaPHy9JGjNmjDp06KC0tDRJ0rRp03TVVVfp6aef1g033KBly5bps88+0/z58+t0XIomACDgjBo1SgcPHtTMmTNVUFCgvn37at26da6LffLy8uR0/u8GkUGDBmnp0qV68MEH9ac//UndunXTW2+9pZ49e9bpuI32Ps1zTXl5udLS0jRjxgzLfgD4Bq+N/+K1QV1RNAEAsInBDQAAsImiCQCATRRNAABsomgCAGATRTOAlZSUaNasWbruuuvUpk0bORwOLV682NfNgqRPP/1UU6ZM0cUXX6wWLVooLi5Ov/rVr7Rz505fN63R+89//qNbb71V559/vpo3b67IyEhdeeWVWr16ta+bhgDAfZoB7NChQ3r44YcVFxenPn36aOPGjb5uEn7yxBNPaPPmzbr11lvVu3dvFRQU6LnnntMll1yiLVu21PneMHjOd999p6NHj2rs2LGKiYnRsWPH9M9//lM33nijXnzxRf3ud7/zdRPhx7jlJICVl5fr8OHDio6O1meffaZLL71UixYt0rhx43zdtEbvo48+Uv/+/RUcHOya9/XXX6tXr14aOXKkXn31VR+2Dj9XWVmphIQElZWVWT4lA5A4PRvQQkJC6jzYMBrGoEGDTAVTkrp166aLL75YX331lY9aBStBQUGKjY3VkSNHfN0U+DlOzwINxDAMFRYW6uKLL/Z1UyCptLRUx48fV1FRkVatWqV3331Xo0aN8nWz4OcomkADee2117Rv3z49/PDDvm4KJN1777168cUXJZ166tHw4cP13HPP+bhV8HcUTaAB7NixQ5MnT1ZiYqLGjh3r6+ZA0j333KORI0dq//79ev3111VZWdlgz45E4KJPE/CygoIC3XDDDQoPD9fKlSsVFBTk6yZBUnx8vJKSkjRmzBitWbNGJSUlSklJEddGojYUTcCLioqKNGzYMB05ckTr1q1TTEyMr5sECyNHjtSnn37KvbSoFadnAS8pKytTSkqKdu7cqffff18XXXSRr5uEWhw/flzSqT90ACskTcALKisrNWrUKGVlZWnFihVKTEz0dZPwkwMHDtSYd+LECS1ZskTNmjXjjxvUiqQZ4J577jkdOXJE+/fvlyStXr1ae/fulSRNnTpV4eHhvmxeo3Xvvfdq1apVSklJ0Y8//lhjMIPbb7/dRy3D//t//0/FxcW68sor1aFDBxUUFOi1117Tjh079PTTT6tly5a+biL8GCMCBbjOnTvru+++c7tsz5496ty5c8M2CJKkIUOGaNOmTZbL+bXznWXLlmnBggXavn27fvjhB7Vq1UoJCQmaOnWqbrzxRl83D36OogkAgE30aQIAYBNFEwAAmyiaAADYRNEEAMAmiiYAADZRNAEAsImiCQCATRRNAABsomgCHnbPPffI4XBo9uzZbpe//fbbatmypZxOpx555JGGbRyAs0LRBDxs27ZtkqS+ffvWWPbkk09q+PDhqqqq0vLly/XQQw81cOsAnA0GbAc8zF3RPHHihCZNmqSFCxcqOjpab7/9tgYMGOCjFgKoL8aeBTwoLy9PnTp1Unh4uI4cOSJJ+vHHHzV8+HBt2rRJvXv31po1axQbG+vbhgKoF07PAh6Uk5MjSerTp48kKTc3VwMHDtSmTZv0y1/+Ups3b6ZgAgGMogl4UPWp2T59+igzM1OJiYnatWuXpk+f7roACEDgomgCHlSdNLdv367rrrtOR48e1Ysvvqi5c+fK6eTXDQh09GkCHtS1a1d98803rq9XrlypESNG+LBFADyJP30BDzl69Kh2796tpk2batiwYZKkxYsXi79LgXMHRRPwkC+++EKGYSg+Pl5Lly7V+eefrzVr1mjOnDm+bhoAD6FoAh5y+pWzEREReuONN9SsWTM9/PDDWr16tW8bB8AjKJqAh/x8UIM+ffpo3rx5MgxDd9xxh77++msftg6AJ1A0AQ/5+T2akjRmzBhNmjRJRUVFuuWWW1RSUuKj1gHwBK6eBTygsrJSrVq10vHjx3XgwAG1a9fOtayiokJXXHGFPvnkE40cOVIrVqzwYUsBnA2SJuABO3fu1PHjx9W+fXtTwZSk4OBgrVy5Uu3atdPKlSv1xBNP+KiVAM4WRRPwgNNHAnInNjZW//jHPxQUFKQHHnhA69evb8jmAfAQTs8CAGATSRMAAJsomgAA2ETRBADAJoomAAA2UTQBALCJogkAgE0UTQAAbKJoAgBgE0UTAACbKJoAANhE0QQAwCaKJgAANlE0AQCw6f8DUfojtmFQxYUAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 500x700 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Community assignments\n",
        "u = np.zeros((60, 3))\n",
        "u[:20, 0] = 1\n",
        "u[20:40, 1] = 1\n",
        "u[40:, 2] = 1\n",
        "\n",
        "plt.matshow(u, aspect='auto', cmap='Blues')\n",
        "plt.gcf().set_size_inches(5, 7)\n",
        "plt.title(r'Community assignments $u$', fontsize=17)\n",
        "plt.xlabel(r'$K$', fontsize=15)\n",
        "plt.ylabel(r'$N$', fontsize=15)\n",
        "plt.xticks(ticks=[0,1,2], labels=[1,2,3], size=12)\n",
        "plt.tick_params(axis='x', bottom=True, top=False, labelbottom=True, labeltop=False)\n",
        "plt.colorbar()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "5b6d41df",
      "metadata": {},
      "outputs": [],
      "source": [
        "sampler = HyMMSBMSampler(\n",
        "    w=w, \n",
        "    u=u,\n",
        "    max_hye_size=10\n",
        ")\n",
        "sample_generator = sampler.sample()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "9a75a25b",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Getting sample number: 0\n",
            "Getting sample number: 1\n",
            "Getting sample number: 2\n",
            "Getting sample number: 3\n",
            "Getting sample number: 4\n",
            "Getting sample number: 5\n",
            "Getting sample number: 6\n",
            "Getting sample number: 7\n",
            "Getting sample number: 8\n",
            "Getting sample number: 9\n",
            "Getting another couple of samples...\n",
            "CPU times: user 7.26 s, sys: 202 ms, total: 7.47 s\n",
            "Wall time: 7.85 s\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "\n",
        "for i in range(10):\n",
        "    print(\"Getting sample number:\", i)\n",
        "    new_sample = next(sample_generator)\n",
        "\n",
        "# Get some more samples later in the code\n",
        "print(\"Getting another couple of samples...\")\n",
        "_ = next(sample_generator)\n",
        "_ = next(sample_generator)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "c023018b",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Extracted sample with N=60 nodes and |E|=851 hyperedges.\n"
          ]
        }
      ],
      "source": [
        "sample = next(sample_generator)\n",
        "print(f\"Extracted sample with N={sample.num_nodes()} nodes and |E|={sample.num_edges()} hyperedges.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1559418c",
      "metadata": {},
      "source": [
        "Notice, however, that all the samples generated from the same call of `HyMMSBMSampler.sample` will have the same degree and size sequence. To have a completely new sample, a new call to the method is needed."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "484387d3",
      "metadata": {},
      "source": [
        "## Conditioning the sampling with additional inputs"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0e3f1f05",
      "metadata": {},
      "source": [
        "#### 1. Providing input sequences"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "596ee699",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Does the sample have same degree sequence as the input one? False\n"
          ]
        }
      ],
      "source": [
        "w = np.eye(3)\n",
        "u = np.zeros((200, 3))\n",
        "u[:66, 0] = 1\n",
        "u[66:133, 1] = 1\n",
        "u[133:, 2] = 1\n",
        "\n",
        "deg_seq = np.random.randint(low=1, high=5, size=200)\n",
        "\n",
        "sampler = HyMMSBMSampler(\n",
        "    w=w, \n",
        "    u=u,\n",
        "    max_hye_size=4,\n",
        ")\n",
        "sample_generator = sampler.sample(deg_seq=deg_seq)\n",
        "sample = next(sample_generator)\n",
        "\n",
        "print(\n",
        "    \"Does the sample have same degree sequence as the input one?\", \n",
        "    np.all(binary_incidence_matrix(sample).sum(axis=1) == deg_seq)\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "692584ae",
      "metadata": {},
      "source": [
        "Notice that, due to the approximations in the MCMC procedure, sometimes the degree sequence in the samples could very slightly deviate from the input one. \n",
        "\n",
        "Similarly, to provide the dimension sequence:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "a7b847f7",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Does the sample have same dimension sequence as the input one? False\n"
          ]
        }
      ],
      "source": [
        "# 10 hyperedges of size 3, 7 hyperedges of size 4, etc...\n",
        "dim_seq = {\n",
        "    3: 10,\n",
        "    4: 7,\n",
        "    7: 5,\n",
        "    9: 5,\n",
        "    12: 3,\n",
        "}\n",
        "\n",
        "sampler = HyMMSBMSampler(\n",
        "    w=w, \n",
        "    u=u,\n",
        "    max_hye_size=4,\n",
        ")\n",
        "sample_generator = sampler.sample(dim_seq=dim_seq)\n",
        "sample = next(sample_generator)\n",
        "\n",
        "print(\n",
        "    \"Does the sample have same dimension sequence as the input one?\", \n",
        "    dict(Counter(len(hye) for hye in sample)) == dim_seq\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9ced670d",
      "metadata": {},
      "source": [
        "#### 2. Providing an input hypergraph"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b4c4f573",
      "metadata": {},
      "source": [
        "## TODO CHANGE DATASET LOADING HERE"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "dc9b1caa",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [],
      "source": [
        "def line_to_hyperedge(line):\n",
        "    hye = line\n",
        "    hye = line.strip(\"\\n\")\n",
        "    hye = [int(node) for node in line.split(\" \")]\n",
        "    return hye\n",
        "\n",
        "# Load Justice dataset.\n",
        "with open(\"./_example_data/justice_data/hyperedges.txt\", \"r\") as hye_file:\n",
        "    with open(\"./_example_data/justice_data/weights.txt\", \"r\") as weight_file:\n",
        "        justice = Hypergraph([\n",
        "            line_to_hyperedge(hye)\n",
        "            for hye, weight in zip(hye_file.readlines(), weight_file.readlines())\n",
        "        ])\n",
        "\n",
        "N = justice.num_nodes()\n",
        "K = 2  # arbitrarily chosen\n",
        "\n",
        "w = np.eye(K)\n",
        "# Random hard community assignments.\n",
        "u = np.zeros((N, K))\n",
        "u[np.arange(N), np.random.randint(0, K, size=N)] = 1  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "4c41945e",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [],
      "source": [
        "model = HyMMSBMSampler(\n",
        "    u=u, \n",
        "    w=w,\n",
        ")\n",
        "sample_generator = model.sample(initial_hyg=justice)\n",
        "_ = next(sample_generator)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "07ba1a2a",
      "metadata": {},
      "source": [
        "## Pre-adjusting the expected statistics"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "6396e020",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "53.75039682539682"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Diagonal affinity matrix\n",
        "w = np.eye(3)\n",
        "\n",
        "# Community assignments\n",
        "u = np.zeros((60, 3))\n",
        "u[:20, 0] = 1\n",
        "u[20:40, 1] = 1\n",
        "u[40:, 2] = 1\n",
        "\n",
        "max_hye_size=10\n",
        "\n",
        "model = HyMMSBM(\n",
        "    u=u, \n",
        "    w=w,\n",
        "    max_hye_size=max_hye_size\n",
        ")\n",
        "orig_deg = model.expected_degree()\n",
        "orig_deg"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6e231f2e",
      "metadata": {},
      "source": [
        "To obtain, for example, an expected degree of 100, one can simply rescale $w$ or $u$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "1d4c61ad",
      "metadata": {
        "pycharm": {
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected degree when rescaling w: 100.0\n",
            "Expected degree when rescaling u: 100.00000000000001\n"
          ]
        }
      ],
      "source": [
        "new_deg = 100\n",
        "\n",
        "rescaled_w = w / orig_deg * new_deg\n",
        "new_w_model = model = HyMMSBM(\n",
        "    u=u, \n",
        "    w=rescaled_w,\n",
        "    max_hye_size=max_hye_size\n",
        ")\n",
        "print(\"Expected degree when rescaling w:\", new_w_model.expected_degree())\n",
        "\n",
        "rescaled_u = u / np.sqrt(orig_deg) * np.sqrt(new_deg)\n",
        "new_u_model = HyMMSBM(\n",
        "    u=rescaled_u, \n",
        "    w=w,\n",
        "    max_hye_size=max_hye_size\n",
        ")\n",
        "print(\"Expected degree when rescaling u:\", new_u_model.expected_degree())"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "hgx-installation",
      "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.11.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}