{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "87d91a42",
   "metadata": {},
   "source": [
    "# Align all and Compute for Graphs\n",
    "\n",
    "$\\textbf{Lead Author: Anna Calissano}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e4098a0",
   "metadata": {},
   "source": [
    "Dear learner, \n",
    "\n",
    "the aim of the current notebook is to introduce the align all and compute as a learning method for graphs. The align all and compute allows to estimate the Frechet Mean, the Generalized Geodesic Principal Components and the Regression. In this notebook you will learn how use all the learning methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "319481d6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Using numpy backend\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "import networkx as nx\n",
    "\n",
    "import geomstats.backend as gs\n",
    "\n",
    "from geomstats.geometry.stratified.graph_space import GraphSpace\n",
    "from geomstats.learning.aac import AAC\n",
    "\n",
    "gs.random.seed(2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e476384d",
   "metadata": {},
   "source": [
    "Let's start by creating simulated data using `networkx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5d2c46cb",
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB8HUlEQVR4nO3dd0BWZf8/8DdDUNziSHOxkSEOFHGjKMuJgJVlpraXlZWWWpallaapaZZmZZYB4gQEEQQURXEgyBYU3IqLDTf3+f3xfTy/eJ5Skfvmusf79d/zIOe8s4A317nO5zKQJEkCEREREdFjMhQdgIiIiIi0GwslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1CAslERERETUICyURERERNQgLJRERERE1iLHoAES6oqxKgfPFZahWKGFibIie5s3R3JRfYkREpPv4046oAXKvlWBrciHisq+j8FY5pL99zABA93Zm8LDriGlu3WHTqaWomERERGplIEmS9PA/RkR/V3SrHB/tSENi3k0YGRqgVvnvX0b3Pz7Muj2+nOyMbu3MGjEpERGR+rFQEtXTtuOF+GT3WSiU0gOL5H8zMjSAsaEBFk9wxFMDuqsxIRERUeNioSSqh7VxuVgendPg68wda4s3PGxUkIiIiEg8vuVN9Ii2HS9USZkEgOXROfjreKFKrkVERCQaCyXRIyi6VY5Pdp9V6TUX7T6LolvlKr0mERGRCCyURI/gox1pUNRjv+SjUCglfLQjTaXXJCIiEoGFkughcq+VIDHvZr1ewHkUtUoJiXk3kXe9RKXXJSIiamwslEQPsTW5EEaGBmq5tpGhAX4/yr2URESk3VgoiR4iLvu6ylcn76tVSojLua6WaxMRETUWFkqiByitUqBQzS/OFBaXo6xKodZ7EBERqRMLJdEDXCgug7oHtUoAzheXqfkuRERE6sNCSfQA1QqlTt2HiIhIHVgoiR7AxLhxvkQa6z5ERETqwJ9iRA/Q07w51PN+999IEi5mnkJNTY2670RERKQWLJRED9Dc1Bjd25mp9R7Kkuvw9vRA+/btERgYiM2bN+Pq1atqvScREZEqsVASPYSHXUe1zqF8wWsgUlJS8N5776GoqAizZs1C586d4erqikWLFuHo0aOora1Vy/2JiIhUwUCSJHW/xEqk1XKvlWDMqgS1XT/mneGw7thS/t83btxAVFQUIiIisG/fPty+fRvm5ubw9vaGr68vvLy8YG5urrY8RERE9cVCSfQIntuUjKT8YpUOODcyNMBgS3NsmeX2r39GoVDg2LFjCA8PR0REBE6fPg1DQ0MMGjQIvr6+8PX1RZ8+fWBgoPadnkRERP+KhZLoERTdKofnynhUqXC8j6mxIWLeGYFu9dijeenSJezbtw8RERHYv38/SkpK0LlzZ/j4+MDX1xeenp5o3bq1yjISERE9ChZKoke07Xgh5oWlqex6X/k7Y+qA7o/9+dXV1Th06BAiIiIQERGBzMxMGBsbY+jQofDz84Ovry969erF1UsiIlI7Fkqielgbl4vl0TkNvs77Y+3wuoe1ChL9fwUFBYiMjERERARiY2NRUVGBHj16yI/GR40aBTMz9b6xTkRE+omFkqieth0vxKJd6aiqUcDA0OiRP09S1sLYEPjSv0+DViYfRUVFBQ4ePIiIiAiEh4ejoKAApqam8PDwkAumlZWVWjMQEZH+YKEkegzPvPgmEiu7wqirE4wMDR74ss79j3cxuItTP36A1KRY2NraNlpWSZKQk5MjPxqPj49HTU0NbG1t5XI5fPhwmJqaNlomIiLSLSyURPWUk5MDBwcHrFixAr5PzcTW5ELE5VxHYXE5/v7FZACgu7kZPGw74tlB3fFkS2M4ODjAzs4OkZGRwvY2lpSUIDY2Vn5z/NKlS2jevDk8PT3h6+sLHx8fdOvWTUg2IiLSTiyURPX0zDPPIDExEbm5uWjatKn8/5dVKXC+uAzVCiVMjA3R07w5mpsa1/ncXbt2YdKkSdi5cycmTpzY2NH/hyRJSEtLk1cvk5KSUFtbC2dnZ3n10t3dHU2aNBEdlYiINBgLJVE9pKWlwcXFBT/88ANeeumlen++JEnw9fVFVlYWMjIy0KxZMzWkfHy3b99GdHQ0IiIiEBkZiRs3bqB169bw8vKCr68vvL290alTJ9ExiYhIw7BQEtWDv78/UlNTkZWV9dirdjk5OXBycsKCBQuwaNEiFSdUHaVSiRMnTsirl8ePH4ckSXB1dZVXLwcMGABDQ57gSkSk71goiR7RiRMn4Orqil9//RXTp09v0LXmzZuH7777DpmZmejZs6dqAqrZ9evX5aHqUVFRuHPnDjp06CAfCTl27Fi0a9dOdEwiIhKAhZLoEfn4+OD8+fNIT0+HkdGjjwv6J6WlpbCzs4ObmxvCwsJUlLDxKBQKHD16VF69TE1NhaGhIdzd3eHr6ws/Pz/07t2bQ9WJiPQECyXRIzh06BCGDRuG4OBgBAYGquSa27Ztw9NPP419+/bBy8tLJdcU5eLFi/JQ9f3796OsrAxdunSRH417enqiZcuWomMSEZGasFASPYQkSfDw8MCdO3dw8uRJle0ZvH/dK1euIC0tDSYmJiq5rmhVVVV1joS8v9902LBhcsG0t7fn6iURkQ5hoSR6iAMHDsDT0xO7d+/G+PHjVXrttLQ09O3bF0uXLsX777+v0mtrivz8fLlcxsXFobKyEhYWFnK5HDlyJI+EJCLSciyURA8gSRLc3d0hSRKOHj2qllW1t99+Gz///DOysrLw5JNPqvz6mqS8vLzOkZDnz59H06ZN4eHhAT8/P/j6+sLCwkJ0TCIiqicWSqIH2Lt3L8aPH4/9+/fD09NTLfe4c+cObG1tMWbMGGzdulUt99BEkiQhKytLXr1MSEiAQqGAvb29vHo5bNgwndkKQESky1goif6FUqlE//790bp1a8TFxal1z9/mzZsxc+ZMxMfHY/jw4Wq7jya7d+8eDhw4IBfMy5cvo0WLFvD09ISfnx98fHx0fgWXiEhbsVAS/YvQ0FAEBgYiISEBw4YNU+u9lEolBg8ejPLycpw8eRLGxsYP/yQdJkkSUlNT5XJ55MgRKJVKuLi4yKuXgwYN0vu/JyIiTcFCSfQP7p9n3b17d+zbt69R7pmSkoKBAwfiu+++w5tvvtko99QWt27dqnMk5M2bN9GmTZs6R0J27NhRdEwiIr3FQkn0D7Zs2YLp06fj+PHjcHV1bbT7vvTSSwgJCUF2djYL0r+ora1FSkqKvHqZkpICAwMDDBgwQF697N+/P4+EJCJqRCyURP+lpqYG9vb26N27N3bs2NGo97558yZsbW3h7++PjRs3Nuq9tdW1a9ewb98+hIeHIzo6Gnfv3kXHjh3h4+MDX19fjBkzBm3bthUdk4hIp7FQEv2Xn376CS+//DJSU1Ph7Ozc6Pdft24dXn/9dSQnJ2PgwIGNfn9tVlNTgyNHjsirl2lpaTAyMsLgwYPl1UtnZ2cOVSciUjEWSqK/qayshI2NDYYOHYo///xTSIba2lq4urrC2NgYycnJfHTbAEVFRYiMjER4eDhiYmJQXl6Orl27yuVy9OjRaNGiheiYRERaj4WS6G9Wr16Nd955B5mZmbC1tRWW4/Dhwxg6dCh++uknzJ49W1gOXVJVVYWEhAR59TInJwcmJiYYPny4XDBtbW25eklE9BhYKIn+o6ysDFZWVvD19cXPP/8sOg6mT5+OyMhI5OTkcA+gGuTl5SEyMlI+ErKqqgqWlpbw9fWFn58fRowYgWbNmomOSUSkFVgoif7j66+/xoIFC5CTk4OePXuKjoMrV67Azs4Ozz//PNasWSM6jk4rKytDXFycfCRkYWEhmjVrhlGjRsmrl5rw3wQRkaZioSQCcPfuXVhaWmLq1KlYt26d6DiyFStW4IMPPsDJkyfh4uIiOo5ekCQJmZmZ8qPxxMREKBQKODg4yOVyyJAhPBKSiOhvWCiJACxevBjLli1DXl6eRh3vV1NTAxcXF5ibmyMhIYH7+wS4e/cuYmJi5IJ59epVtGzZEmPGjIGvry98fHzQpUsX0TGJiIRioSS9V1xcDEtLS8yePRsrVqwQHed/xMTEYMyYMfj9998xbdo00XH0mlKpRGpqKsLDwxEREYGjR49CkiT07dtXXr10c3ODkZGR6KhERI2KhZL03vz587FmzRrk5+dr7Ok0gYGBOHz4MLKzs9GyZUvRceg/bt68KR8JuW/fPhQXF6Ndu3bykZBeXl7o0KGD6JhERGrHQkl67erVq7CyssKcOXPwxRdfiI7zrwoLC2Fvb4833ngDX3/9teg49A9qa2tx/Phx+cWekydPwsDAAAMHDoSfnx98fX3Rt29fzhUlIp3EQkl6bc6cOfjll19QUFCg8aN5lixZgsWLFyMtLQ329vai49BDXLlyBfv27UNERASio6Nx7949dOrUST4ScuzYsWjdurXomEREKsFCSXqrqKgI1tbWWLhwIRYsWCA6zkNVVlbC0dERlpaWiI6O5gs6WqSmpgZJSUny6uXZs2dhZGSEoUOHynsvHR0d+e+UiLQWCyXprVdeeQWhoaEoKCjQmn2Je/fuxfjx47F9+3b4+/uLjkOP6cKFC/JQ9QMHDqC8vBzdunWrcyRk8+bNRcckInpkLJSkl/Lz82FnZ4elS5di7ty5ouPUy7hx45CWlobMzEyYmZmJjkMNVFlZiYSEBISHhyM8PBznzp2DiYkJRo4cKRdMGxsb0TGJiB6IhZL00vPPP4/o6GicO3dO60pZXl4eHB0dMW/ePCxevFh0HFKx3NxceeblwYMHUV1dDWtra7lcjhgxAk2bNhUdk4ioDhZK0juZmZlwcnLC6tWr8frrr4uO81g+/vhjrFixAhkZGbC0tBQdh9SktLQUcXFx8tzLoqIimJmZYfTo0fJQ9R49eoiOSUTEQkn6Z+rUqTh69ChycnJgamoqOs5jKSsrg729Pfr164ddu3aJjkONQJIknD17Vl69PHToEGpra+Ho6FjnSMgmTZqIjkpEeoiFkvTK6dOn0bdvX2zcuBGzZs0SHadBQkJCEBQUhPDwcPj6+oqOQ43szp07iImJQXh4OCIjI3Ht2jW0atUKY8eOha+vL7y9vdG5c2fRMYlIT7BQkl6ZMGECMjMzkZGRofUrOZIkwdPTE4WFhUhPT9fa1VZqOKVSiVOnTsmrl8nJyZAkCf369YOvry/8/PwwYMAAHglJRGrDQkl6Izk5GYMGDcLWrVvxzDPPiI6jEhkZGXBxccHnn3+OefPmiY5DGuLmzZuIiopCeHg49u3bh9u3b8Pc3Bze3t7ykZDm5uaiYxKRDmGhJL0xduxYXL58GampqTq1UvPuu+9iw4YNyM7ORteuXUXHIQ1TW1uL5ORkefXy1KlTMDQ0hJubm7z3sm/fvhyqTkQNwkJJeiE+Ph4jR47UyYHgd+/ehZ2dHUaOHIlt27aJjkMa7vLly/JQ9f3796OkpASdO3eWj4T09PTkkZBEVG8slKTzJEnC8OHDUV5ejpSUFJ1cifntt9/w/PPPIzY2Fh4eHqLjkJaorq7G4cOH5dXLjIwMGBsbY+jQofDz84Ovry969eqlk18zRKRaLJSk86KiouDt7Y2IiAj4+PiIjqMWSqUSw4YNw927d3Hq1Cmtf+GIxDh//rxcLmNjY1FRUYEePXrIj8Y9PDx4JCQR/SMWStJpkiRh4MCBMDExwaFDh3R6peXUqVPo378/Vq5cibffflt0HNJyFRUViI+PR0REBMLDw5Gfnw9TU1OMHDlSXr20srISHZOINAQLJem0nTt3YvLkyXrzKPi1117D1q1bkZOTg06dOomOQzpCkiTk5ubKJ/bEx8ejpqYGtra28url8OHDObqKSI+xUJLOUiqVcHFxQceOHXHgwAHRcRpFcXExbG1tMWHCBGzevFl0HNJRpaWlOHDggPx4/OLFi2jevLl8JKSvry+6desmOiYRNSIWStJZ27Ztw9NPP42kpCS4u7uLjtNoNmzYgFdeeUXv/rlJDEmSkJ6eLq9eJiUloba2Fs7OznK5dHd3575eIh3HQkk6SaFQwNHRETY2Nti7d6/oOI2qtrYWAwcOBAAcO3ZMp2Zukua7ffs29u/fj4iICERGRuL69eto3bo1vLy85CMhuR2DSPewUJJO2rx5M2bOnIkTJ06gX79+ouM0uiNHjmDw4MHYsGEDXnrpJdFxSE8plUqcOHFCfjR+/PhxSJIEV1dXefXS1dWVv/QQ6QAWStI51dXVsLW1haurK0JDQ0XHEeaFF17A7t27kZOTw2P2SCNcv34dUVFRiIiIwL59+3Dnzh20b99eHqo+duxYtGvXTnRMInoMLJSkc9avX4/XX38d6enpcHBwEB1HmGvXrsHW1hbTpk3DunXrRMchqkOhUODo0aPy6mVqaioMDQ3h7u4ur166uLjo9KgvIl3CQkk6paKiAtbW1hg1ahS2bNkiOo5wq1atwrvvvouUlBS9fPRP2uPSpUt1joQsLS1Fly5d5HLp6emJli1bio5JRP+ChZJ0ysqVK/H+++8jKysL1tbWouMIV1NTg759+6JVq1Y4dOgQDA0NRUcieqjq6mocOnRIfnM8KysLTZo0wbBhw+SCaW9vz9VLIg3CQkk6o7S0FJaWlpg4cSJ++ukn0XE0RlxcHEaNGoVff/0V06dPFx2HqN7y8/Pl1cvY2FhUVlbCwsJCLpcjR46EmZmZ6JhEeo2FknTGl19+icWLFyM3Nxfdu3cXHUejPPXUUzh48CBycnLQqlUr0XGIHlt5eTkOHjwoHwl5/vx5NG3aFB4eHnLBtLS0FB2TSO+wUJJOuHPnDiwsLPDcc89h9erVouNonIsXL8LOzg6vvPIKVqxYIToOkUpIkoTs7Gz5xZ6EhATU1NTA3t5eLpfDhg2DiYmJ6KhEOo+FknTCokWLsHz5cpw7dw6dO3cWHUcjLV26FAsXLkRqaiocHR1FxyFSuXv37tU5EvLy5cto0aIFPD094evrCx8fH3Tt2lV0TCKdxEJJWu/GjRuwtLTEq6++iq+//lp0HI1VVVUFJycndO/eHTExMXyhgXSaJEk4c+aMXC6TkpKgVCrh4uIir14OGjQIxsbGoqMS6QQWStJ677//PjZs2ICCggIO8H6IyMhI+Pr6Ijg4GIGBgaLjEDWaW7duITo6Wj4S8ubNm2jTpk2dIyE7duwoOiaR1mKhJK125coVWFpa4oMPPsDixYtFx9EKEydOxKlTp5CZmYnmzZuLjkPU6JRKJVJSUuocCWlgYIABAwbIq5f9+/fnmC2iemChJK325ptvYuvWrcjPz0ebNm1Ex9EK+fn5cHBwwNy5c7FkyRLRcYiEu3btGvbt24eIiAhERUXh7t276NChQ50jIdu2bSs6JpFGY6EkrXXhwgXY2Nhg8eLFmD9/vug4WmXRokX46quvcPbsWQ6AJ/obhUKBI0eOyGOJ0tLSYGRkhMGDB8url87OztyDTPRfWChJa82ePRt79uxBfn4+H93WU3l5OXr16gVnZ2fs3btXdBwijVVUVCQPVY+JiUFZWRmefPLJOkdCtmjRQnRMIuFYKEkr5ebmolevXli+fDnmzJkjOo5WCgsLw5QpU7Bnzx6MGzdOdBwijVdVVYXExER572V2djaaNGmCESNGyAXT1taWq5ekl1goSSs9++yzOHjwIPLy8tC0aVPRcbSSJEnw8vJCfn4+0tPT+fdIVE95eXny6mVcXByqqqpgaWkJX19f+Pn5YcSIEWjWrJnomESNgoWStE56ejp69+6NdevW4ZVXXhEdR6tlZWXB2dkZn376KT7++GPRcYi0Vnl5OeLi4uS9lxcuXECzZs0watQoefWyZ8+eomMSqQ0LJWmdKVOm4NSpU8jKyuKRairwwQcfYO3atcjKyuIZ6EQqIEkSMjMz5UfjiYmJUCgU6NWrl7x6OWTIEH7/Ip3CQkla5cSJE3B1dcUvv/yC559/XnQcnVBSUgI7OzsMGTIEISEhouMQ6Zx79+4hJiZGLphXrlxBy5YtMWbMGPlIyC5duoiOSdQgLJSkVfz8/HDu3Dmkp6fzyDQV2rp1K5599lns378fnp6eouMQ6SxJknD69Gm5XB49ehRKpRJ9+vSRVy/d3NxgZGQkOipRvbBQktZISkrCkCFDsG3bNkydOlV0HJ0iSRKGDx+O4uJipKamokmTJqIjEemF4uLiOkdCFhcXo127dvKRkF5eXujQoYPomEQPxUJJWmPUqFEoLi7GqVOneCSaGqSmpqJfv3745ptv8O6774qOQ6R3amtrcfz4cXn18sSJEzAwMMDAgQPh5+cHX19f9O3bl9//SCOxUJJWiI2NxejRo7Fr1y5MmDBBdByd9eabb+LXX39FdnY2OnfuLDoOkV67evVqnSMh7927h06dOslHQo4ZM4ZHzpLGYKEkjSdJEoYMGQKFQoHk5GQODVaj27dvw9bWFj4+Pvjtt99ExyGi/6ipqUFSUpK8epmeng4jIyMMGTJEXr10dHTk90cShoWSNF54eDjGjRuHqKgojB07VnQcnbdx40a8+OKLSExMxNChQ0XHIaJ/UFhYWOdIyPLycnTr1k2eeTlq1CgeCUmNioWSNJpSqYSrqytatmyJgwcP8rfvRqBUKjFo0CDU1NQgJSWFb5sSabjKykokJCTIq5e5ubkwMTHBiBEj5NVLGxsb0TFJx7FQkkbbvn07AgICEB8fj+HDh4uOozeOHTsGNzc3rFu3Dq+++qroOERUD7m5uYiMjER4eDgOHjyI6upqWFtby6uXI0aM4FGrpHIslKSxamtr0bt3b3Tt2hVRUVGi4+id2bNnIywsDDk5OWjfvr3oOET0GMrKyhAbGysfCVlUVAQzMzOMHj1aHqreo0cP0TFJB7BQksb6/fff8dxzzyE5ORkDBw4UHUfv3LhxA7a2tggKCsKGDRtExyGiBpIkCRkZGXK5PHToEGpra+Ho6CivXg4ZMoRzaOmxsFCSRqqpqUGvXr3g5OSEnTt3io6jt9asWYO3334bx44dg6urq+g4RKRCd+/exf79++W9l9euXUOrVq0wduxY+Pr6wtvbm+PD6JGxUJJGuv+mcWpqKnr37i06jt5SKBTo168fzMzMkJSUxIHKRDpKqVTKR0KGh4cjOTkZkiShX79+8urlwIED+ZIe/SsWStI4VVVVsLGxweDBg7Ft2zbRcfReQkICRowYgc2bN2PGjBmi4xBRI7h58yaioqIQERGBffv24datWzA3N4e3t7d8JKS5ubnomKRBWChJ46xZswZz5szB2bNnYW9vLzoOAZg2bRpiYmKQnZ3NkzmI9ExtbS2OHTuG8PBwREREyMffurm5yauXffv25Vg3PcdCSRqlvLwclpaW8PHxwebNm0XHof+4fPky7OzsMGvWLKxatUp0HCIS6PLly/KRkNHR0SgpKUHnzp3lIyE9PT3RunVr0TGpkbFQkkb55ptv8NFHHyEnJwcWFhai49DffPPNN5g/fz5OnToFZ2dn0XGISANUV1cjKSlJXr3MyMiAsbExhg4dKq9eOjg4cPVSD7BQksa4d+8eLC0tERgYiPXr14uOQ/+luroavXv3xhNPPIG4uDj+gCCi/3H+/Hn5SMgDBw6goqICPXr0kMulh4cHmjdvLjomqQELJWmMzz77DF9++SXy8vLQtWtX0XHoH0RHR8PLywt//vknnnrqKdFxiEiDVVZWIj4+HuHh4QgPD0d+fj5MTU0xcuRIuWBaW1uLjkkqwkJJGuHWrVuwsLDArFmz8O2334qOQw/g7++PY8eOISsrCy1atBAdh4i0gCRJyM3NlWdexsfHo7q6Gra2tnK5HD58OExNTUVHpcfEQkka4aOPPsJ3332HgoICdOzYUXQceoDz58+jV69emDNnDpYuXSo6DhFpodLS0jpHQl68eBHNmzeXj4T09fVFt27dRMekemChJOGuX78OCwsLvP322/jyyy9Fx6FH8Nlnn2HJkiVIT0+Hra2t6DhEpMUkSUJ6erq8enn48GHU1tbC2dlZLpfu7u4adyRkWZUC54vLUK1QwsTYED3Nm6O5qbHoWMKwUJJw77zzDn7++WcUFBSgXbt2ouPQI6ioqICDgwPs7OwQGRnJF3SISGXu3LmD6OhoREREIDIyEtevX0fr1q3lIyF9fHzQqVMnIdlyr5Vga3Ih4rKvo/BWOf5eoAwAdG9nBg+7jpjm1h02nVoKySgKCyUJdfHiRVhbW+Ojjz7CokWLRMeheti1axcmTZqEnTt3YuLEiaLjEJEOUiqVOHnypLx6eezYMUiSBFdXV3n10tXVVe1HQhbdKsdHO9KQmHcTRoYGqFX+e3W6//Fh1u3x5WRndGtnptZsmoKFkoR69dVXERISgvz8fLRq1Up0HKoHSZLg6+uL7OxsnD17Fs2aNRMdiYh03I0bN+Sh6lFRUbh9+zbat29f50hIVT/p2na8EJ/sPguFUnpgkfxvRoYGMDY0wOIJjnhqQHeVZtJELJQkTEFBAWxtbfHll1/i/fffFx2HHkNOTg6cnJywYMECrjATUaNSKBRITk6WVy9Pnz4NQ0NDuLu7y6uXLi4uDdqSszYuF8ujcxqcde5YW7zhYdPg62gyFkoSZsaMGdi3bx/y8/NhZqYfjwR00fz587Fq1SpkZmaiZ8+eouMQkZ66dOkS9u3bh/DwcOzfvx+lpaXo0qWLXC49PT3RsuWj72vcdrwQ88LSVJbvK39nTNXhlUoWShIiKysLjo6OWLVqFd58803RcagBSktLYW9vj4EDByIsLEx0HCIiVFdX49ChQ/LqZWZmJpo0aYJhw4bJBdPe3v5fVy+LbpXDc2U8qhRKlWUyNTZEzDsjdHZPJQslCfHUU08hKSkJubm5HGSrA7Zt24ann34aUVFRGDt2rOg4RER1FBQUyOUyNjYWlZWV6NmzJ3x9feHn54eRI0fWeVL23KZkJOUX12vP5MMYGRpgsKU5tsxyU9k1NQkLJTW6M2fOwMXFBT/99BNmz54tOg6pgCRJ8PDwwNWrV3HmzBmYmJiIjkRE9I8qKipw8OBBeah6QUEBmjZtCg8PD/j6+sLBfTRmhuar7f4x7wyHdUfdGynEQkmNbtKkSUhPT5cfQZBuSEtLQ9++fbF06VK+ZEVEWkGSJGRnZ8urlwkJCWgx4gW06ucHGKp+FJGRoQGec+uBTyc4qvzaorFQUqM6duwY3NzcsGXLFjz77LOi45CKzZkzB5s2bUJWVhaefPJJ0XGIiOqlpKQEI1fEo7hKfYc19DA3Q/xcD7VdXxQWSmpUXl5euHjxIs6cOaP2QbTU+O7cuQM7Ozt4enpi69atouMQEdVLaZUCzp9GQZ3FyABA+qdeOndMo6HoAKQ/EhISEB0djc8++4xlUke1adMGy5Ytwx9//IGEhATRcYiI6uVCcZlayyQASADOF5ep+S6NjyuU1CgkScKIESNQWlqKlJQUGBrydxldpVQqMXjwYJSXl+PkyZMwNtat38KJSHedKryNyeuT1H6fHa8ORt/ubdV+n8bEn+rUKPbv34/ExER8/vnnLJM6ztDQEGvXrkV6ejp++OEH0XGIiB6ZiXHj/HxqrPs0Jq5QktpJkgQ3NzcYGxvj8OHDDToGi7THyy+/jODgYGRnZ6Njx46i4xARPVRZlQJO3EP5WHSvIpPG2bNnD44fP44lS5awTOqRL774AgYGBvjoo49ERyEieqjS0lLs3bkdJtX31Hqf7uZmOlcmARZKUjOlUomFCxfCw8MDo0aNEh2HGlH79u3xxRdfYNOmTTh27JjoOERE/6OsrAwhISEIDAxEx44d8dRTT8HwaiYM1LRGaWRoAA9b3Xxiw0fepFZ//fUXnnrqKRw+fBiDBw8WHYcaWW1tLVxdXdGkSRMcPXqU+2eJSLjy8nJEREQgODgY4eHhKC8vR//+/REUFISAgADUNu+AMavUN6WCJ+UQ1ZNCoYCTkxOsrKwQHh4uOg4JcvjwYQwdOhQbN27ErFmzRMchIj1UXl6OyMhIBAcHY+/evSgvL0e/fv3kEmllZVXnz/Ms7/pjoSS1+fXXXzFjxgykpKSgf//+ouOQQNOnT0dkZCRycnLQtq1ujcogIs1UUVFRp0SWlZWhb9++CAwMRGBgIKytrf/1c4tulcNzZTyqFEqV5TE1NkTMOyPQrZ2Zyq6pSVgoSS2qq6thb2+Pvn37Yvv27aLjkGBXrlyBnZ0dnn/+eaxZs0Z0HCLSURUVFdi3bx9CQkKwe/dulJWVoU+fPnKJtLGxeeRrbTteiHlhaSrL9pW/M6YO6K6y62kaFkpSix9++AGvvfYa0tLS4OjoKDoOaYBvv/0W77//Pk6ePAkXFxfRcYhIR1RWVtYpkaWlpejduzeCgoIQGBgIW1vbx7722rhcLI/OaXDG98fa4XWPf18R1QUslKRyFRUVsLa2hoeHB37//XfRcUhD1NTUwMXFBe3bt0d8fDxHSBHRY6usrERUVJRcIktKSuDs7CyXSDs7O5Xda9vxQnyy+yyqahSAwaO/WGhkaABjQwN8NsFRp1cm72OhJJVbtWoV5s6di8zMzHo9XiDdFxMTgzFjxmDr1q145plnRMchIi1SVVWF6OhoBAcHY9euXSgpKYGTk5NcIu3t7dV277DoBLz6czyaWfSDkaHBA1/Wuf/xYdbt8eVkZ53dM/nfWChJpUpLS2FlZYXx48dj48aNouOQBgoMDMThw4eRnZ2Nli11b3QGEalOVVUV9u/fL5fIe/fuwdHRUS6RvXr1UnsGSZLg4eGBu3fvYlvkQfx57CLicq6jsLi8zrRKA/zf0HIP2454dlB3nRwN9CAslKRSS5cuxSeffILc3Fz06NFDdBzSQIWFhbC3t8cbb7yBr7/+WnQcItIw1dXVdUrk3bt30atXL0ydOhWBgYFwcHBo1Dz3n6zs2bMH48aNk///sioFzheXoVqhhImxIXqaN9fJE3AeFQslqcydO3dgYWGBZ599lm/y0gN98cUX+PTTT5GWlqbWx1REpB2qq6sRExMjl8g7d+7A3t4eQUFBCAoKEvZypyRJcHd3BwAcOXKEe78fgIWSVOaTTz7B119/jfz8fHTu3Fl0HNJglZWVcHJygqWlJaKiovhNmkgPVVdX48CBAwgJCcGOHTtw584d2NnZ1SmRor837N27F+PHj0dMTAxGjx4tNIumY6Eklbh58yYsLS3x8ssv45tvvhEdh7TA/W/UYWFhmDx5sug4RNQIampq6pTI27dvw9bWVi6RTk5OwkvkfUqlEv3790ebNm0QGxurMbk0FQslqcQHH3yA9evXo6CgAO3btxcdh7TEuHHjkJ6ejoyMDJiZ6cebkET6pqamBnFxcQgODsaOHTtw69Yt2NjYyC/W9O7dWyPLWkhICIKCgpCYmIihQ4eKjqPxWCipwa5cuQIrKyvMnTsXn332meg4pEXy8vLg6OiIefPmYfHixaLjEJGKKBQKuUSGhYXh1q1bsLa2lkuki4uLRpbI+2pra+Hk5ISePXsiMjJSdBytwEJJDfbWW29hy5YtKCgoQJs2bUTHIS2zYMECLF++HBkZGbC0tBQdh4gek0KhwMGDB+USWVxcDCsrKwQGBiIoKAh9+vTR6BL5d1u2bMH06dORkpKC/v37i46jFVgoqUEKCwthY2ODTz75BB999JHoOKSFysrK0KtXL/Tt2xe7du0SHYeI6kGhUCA+Pl4ukff3098vkX379tWaEnlfTU0N7O3t4eLigrCwMNFxtIb+Dkwilfj888/RunVrvPXWW6KjkJZq3rw5VqxYgaCgIERGRsLHx0d0JCJ6AIVCgYSEBLlE3rhxAxYWFpg5cyaCgoLQr18/rSuRf7d582YUFBRg586doqNoFa5Q0mPLy8uDvb09vvnmG7zzzjui45AWkyQJnp6eKCoqQlpaGkxNTUVHIqK/qa2tRUJCAkJCQrB9+3Zcv34dPXr0kN/O7t+/v1aXyPsqKythY2ODYcOG4Y8//hAdR6uwUNJje+655xAbG4u8vDw0a9ZMdBzSchkZGXBxccHnn3+OefPmiY5DpPdqa2uRmJgol8hr166he/fucol0dXXViRL5d6tXr8a7776LjIwM2Nraio6jVVgo6bFkZGTAyckJ33//PV599VXRcUhHvPfee/jhhx+QnZ2Nrl27io5DpHdqa2tx+PBhBAcHIzQ0FNeuXUO3bt3kt7MHDhyocyXyvrKyMlhZWcHPzw+bNm0SHUfrsFDSYwkICMCJEyeQnZ0NExMT0XFIR9y7dw+2trYYOXIktm3bJjoOkV5QKpV1SuTVq1fRtWtXuUS6ubnpbIn8u6+++goLFy5ETk4OevbsKTqO1mGhpHo7efIk+vfvj82bN2PGjBmi45COuT+uIzY2Fh4eHqLjEOkkpVKJpKQkuUReuXIFXbt2RUBAAIKCguDm5gZDQ0PRMRvN3bt3YWFhgaeffhrff/+96DhaiYWS6m3cuHHIzc3F2bNnYWzMQQGkWpIkYejQobh37x5OnjyJJk2aiI5EpBOUSiWOHDkil8jLly/jySeflEvkoEGD9KpE/t3ixYuxbNkynDt3Dl26dBEdRyuxDVC9HDlyBOHh4fjzzz9ZJkktDAwMsHbtWvTv3x/r1q3D22+/LToSkdZSKpU4evQoQkJCEBISgkuXLqFz587ynEh3d3e9LZH3FRcX49tvv8Xrr7/OMtkAXKGkehk9ejRu3LiB06dP6/03IVKv1157DVu3bkVOTg46deokOg6R1lAqlUhOTpZL5MWLF9G5c2d5JXLw4MH8/v038+bNw9q1a1FQUIAOHTqIjqO1WCjpkcXGxmL06NHYuXMnJk6cKDoO6bhbt27B1tYW48ePx+bNm0XHIdJokiTVKZFFRUV44okn5BI5ZMgQlsh/cPXqVVhaWuLdd9/FkiVLRMfRaiyU9Eju72urrq7GsWPH9OKNPxLvxx9/xMsvv4ykpCS4u7uLjkOkUSRJwvHjxxEcHIyQkBAUFhaiU6dOCAgIQGBgIIYOHQojIyPRMTXanDlz8Msvv6CgoABt27YVHUersVDSI4mMjISvry/27dsHLy8v0XFIT9TW1sLNzQ0AkJyczB+OpPckSUJKSopcIi9cuICOHTvKJXLYsGH8OnlERUVFsLa2xqJFi/Dxxx+LjqP1WCjpoSRJgqurK8zMzJCQkMDVSWpUR48ehbu7OzZs2ICXXnpJdByiRidJEk6cOCGXyPPnz6NDhw6YMmUKgoKCMHz4cJbIx/Dyyy8jLCwM+fn5aNmypeg4Wo+Fkh4qLCwMU6ZMQXx8PIYPHy46DumhF154Abt370ZOTg7Mzc1FxyFSO0mScPLkSblEFhQUoH379nVKJCdtPL5z587B3t4ey5Ytw3vvvSc6jk5goaQHqq2thYuLC7p06YLo6GjRcUhPXbt2Dba2tpg2bRrWrVsnOg6RWkiShFOnTiEkJATBwcHIz8+Hubm5XCJHjBjBEqkizz//PPbv349z586hWbNmouPoBP6XSQ/0119/4ezZszzXlITq1KkTPvvsM7zzzjuYPXs2+vXrJzoSkUpIkoTTp0/LJfLcuXMwNzeHv78/NmzYgJEjR7JEqlhmZiZ+//13rFmzhmVShbhCSf+qpqYGDg4O6NWrF3bv3i06Duk5hUKBvn37olWrVjh06BD38pLWkiQJZ86cQXBwMIKDg5GXl4d27drB398fgYGB8PDw4AlRahQUFIRjx44hOzsbpqamouPoDP7aQ//qt99+Q15eHkJDQ0VHIYKxsTHWrFkDDw8P/P7773juuedERyJ6ZJIkIS0tTS6Rubm5aNu2Lfz9/bF27VqMGjWKJbIR3F8N3rRpE8ukinGFkv5RVVUVbG1tMWjQIPz111+i4xDJnnrqKRw8eBA5OTlo1aqV6DhE/0qSJKSnp8slMicnB23btsXkyZMRGBiI0aNHs0Q2sgkTJiArKwsZGRncSqBiLJT0j77//nu89dZbSE9PR69evUTHIZJdvHgRdnZ2eOWVV7BixQrRcYjqkCQJZ8+eld/OzsrKQps2bTBp0iQEBQVh9OjRMDExER1TLyUnJ2PQoEH4448/8PTTT4uOo3NYKOl/lJeXw8rKCmPHjsWvv/4qOg7R/1i2bBkWLFiA1NRUODo6io5DhLNnz8ov1mRmZqJ169ZyifT09GSJ1ABjxozB1atXkZqaymMo1YCFkv7H8uXLMX/+fGRnZ8PS0lJ0HKL/UVVVBWdnZ3Tr1g0xMTF8QYeEyMjIkEtkRkYGWrVqVadEco+e5jh48CA8PDwQFhaGyZMni46jk1goqY6SkhJYWFggICAAP/zwg+g4RP9q37598PHxQUhICAICAkTHIT2RmZkpl8izZ8+iVatWmDhxIoKCgjBmzBiWSA0kSRKGDx+OiooKHD9+nL+AqgkLJdWxZMkSLFmyBHl5eejatavoOEQPNHHiRJw6dQqZmZlo3ry56Diko7Kzs+UXa9LT09GyZUtMnDgRgYGBGDt2LJo2bSo6Ij1AVFQUvL29ERkZCW9vb9FxdBYLJclu374NCwsLzJgxA6tWrRIdh+ih8vPz4eDggLlz52LJkiWi45AOycnJkV+sOXPmDFq0aCGXSC8vL5ZILSFJEgYOHAhTU1MkJiZydVKN+M48yZYvX46amhrMnz9fdBSiR2JpaYkPP/wQy5Ytw4wZM2BtbS06Emmx3NxcuUSmpqaiRYsWGD9+PBYvXgwvLy+eqqKFdu3ahZSUFMTFxbFMqhlXKAkAcP36dVhaWuLNN9/E0qVLRcchemTl5eVwcHCAk5MT9u7dKzoOaZm8vDx5T+Tp06fRvHlzjB8/HkFBQfD29maJ1GJKpRIuLi7o1KkTYmJiRMfReVyhJADAV199BSMjI7z//vuioxDVi5mZGVauXAl/f3/s3bsX48aNEx2JNNy5c+fkEnnq1CmYmZlh/PjxWLhwIXx8fFgidcT9Pa8//fST6Ch6gSuUhEuXLsHa2hrz5s3DJ598IjoOUb1JkgQvLy/k5+cjPT2d+9vof+Tn58sl8uTJkzAzM8O4ceMQFBQEHx8fmJmZiY5IKqRQKODo6AhbW1vs2bNHdBy9wEJJeO211/DXX3+hoKCAR9mR1srKyoKzszM+/fRTfPzxx6LjkAYoKCiQS+SJEyfQrFmzOiWSkwF01+bNmzFz5kycPHkSffv2FR1HL7BQ6rmCggLY2dlhyZIl+OCDD0THIWqQDz74AGvXrkVWVha6d+8uOg4JcP78eblEpqSkoFmzZvDz80NgYCD8/PxYIvVAdXU1bG1tMWDAAISEhIiOozdYKPXczJkzERERgXPnzvEbLWm9kpIS2NnZYciQIfxBokcuXLiAkJAQhISE4NixY2jatCl8fX0RFBQEPz8/tGjRQnREakTr16/HG2+8gfT0dPTq1Ut0HL3BQqnHcnJy0KtXL6xcuRJvvfWW6DhEKvHHH39g2rRpiImJwejRo0XHITUpLCyUS2RycjKaNm0KHx8fBAUFYdy4cSyReqqiogLW1tYYPXo0fvvtN9Fx9AoLpR57+umncejQIeTm5vIlBtIZkiRhxIgRuHnzJlJTU9GkSRPRkUhFioqKEBoaiuDgYBw9ehSmpqZ1SmTLli1FRyTBvv32W3zwwQfIzs6GlZWV6Dh6hYVST505cwYuLi748ccf8eKLL4qOQ6RSqamp6NevH7755hu8++67ouNQA1y8eFEukUeOHIGJiUmdEskXCem+0tJSWFhYYPLkyfjxxx9Fx9E7LJR6avLkyThz5gyysrK4gkM66c0338Svv/6K7OxsdO7cWXQcqodLly7JJTIpKQkmJibw9vZGUFAQxo8fzxJJ/+jLL7/E4sWLkZeXh27duomOo3dYKPVQSkoKBgwYgN9++w3PPfec6DhEanH79m3Y2trCx8eHe6m0wKVLl7B9+3YEBwfj8OHDMDExgZeXFwIDAzFhwgS0bt1adETSYHfu3IGFhQWmT5+O7777TnQcvcRCqYe8vb1RWFiItLQ0GBkZiY5DpDabNm3C7NmzcejQIQwZMkR0HPovly9frlMijY2N65TINm3aiI5IWmLhwoVYsWIF8vPz8cQTT4iOo5dYKPVMYmIihg8fjpCQEAQEBIiOQ6RWSqUSgwYNQk1NDVJSUvgLlAa4cuUKtm/fjpCQECQmJsLY2BhjxoxBUFAQJk6cyBJJ9Xbjxg1YWlritddew1dffSU6jt5iodQjkiRh5MiRuHfvHk6cOAFDQ0PRkYjU7tixY3Bzc8O6devw6quvio6jl65evYqwsDAEBwcjISEBRkZGdUpk27ZtRUckLfb+++9jw4YNKCgogLm5ueg4eouFUo/ExMRgzJgx2LNnD8aNGyc6DlGjmT17NsLCwpCTk4P27duLjqMXrl27JpfI+Ph4GBkZwdPTUy6R7dq1Ex2RdMCVK1dgaWmJDz/8EJ9++qnoOHqNhVJPSJKEQYMGwcDAAEeOHIGBgYHoSESN5saNG7C1tUVQUBA2bNggOo7Oun79ep0SaWBgIJfISZMmsUSSyr3xxhv4448/UFBQwBe3BDMWHYAax969e3Hs2DHExMSwTJLe6dChAz7//HO89dZbePHFF+Hq6io6ks64ceOGXCIPHjwIAwMDjB49Gj/++CMmTZrER5CkNhcuXMCPP/6Izz77jGVSA3CFUg8olUr069cPbdu2RWxsLAsl6SWFQoH+/fvDzMwMhw8f5h7iBrhx4wZ27NiB4OBgxMXFwcDAAKNGjUJgYCAmT57MbQXUKGbPno09e/YgPz8fzZs3Fx1H73GFUg9s374dqampSExMZJkkvWVsbIy1a9di+PDh+O233zBjxgzRkbTKzZs3sWPHDoSEhCA2NhaSJMHDwwPr16/H5MmT0aFDB9ERSY/k5ubil19+wYoVK1gmNQRXKHVcbW0tnJyc0LNnT0RGRoqOQyTctGnTEBMTg+zsbI6oeYji4mLs3LkTwcHBOHDggDwpIigoCJMnT0bHjh1FRyQ9NW3aNMTHxyMvLw9NmzYVHYfAQqnzfvvtNzz//PNISUlB//79RcchEu7y5cuws7PDrFmzsGrVKtFxNM6tW7fqlEilUokRI0YgKCgI/v7+LJEkXHp6Onr37o3169fj5ZdfFh2H/oOFUofV1NTAzs4Offr0QVhYmOg4RBrjm2++wfz583H69Gk4OTmJjiPc7du35RIZExOD2traOiWyU6dOoiMSyaZMmYJTp04hKysLJiYmouPQf7BQ6rAff/wRr7zyClJTU+Hs7Cw6DpHGqK6uRu/evdG5c2e9fVHt9u3b2LVrF4KDg7F//37U1tZi+PDhconk8XWkiU6cOAFXV1f8+uuvmD59uug49DcslDqqsrISNjY2GDZsGP744w/RcYg0TnR0NLy8vLBt2zZMnTpVdJxGcefOnTolUqFQYNiwYQgMDMSUKVPQuXNn0RGJHsjX1xcFBQVIT0/nUaoahoVSR3333Xd47733kJGRAVtbW9FxiDSSv78/jh07hqysLLRo0UJ0HLW4e/euXCKjo6OhUCgwdOhQuUR26dJFdESiR3L48GEMHToUf/31F4KCgkTHof/CQqmDysrKYGlpiXHjxmHTpk2i4xBprPPnz6NXr16YM2cOli5dKjqOyty9exe7d+9GSEgIoqKiUF1djSFDhiAoKAhTpkzBk08+KToiUb2NGjUKt27dwsmTJzlHVgNxDqUOWrt2LW7fvo2FCxeKjkKk0Xr27In58+djyZIleOGFF7R6Nf/evXvYs2cPgoODsW/fPlRXV2Pw4MH4+uuvMWXKFHTt2lV0RKLHFhsbi7i4OOzevZtlUkNxhVLH3L17FxYWFnj66afx/fffi45DpPEqKirg6OgIOzs7REREaNULOiUlJXVKZFVVFdzd3eWVyG7duomOSNRgkiRh8ODBUCqVOHr0qFZ9jeoTrlDqmJUrV6KiogIff/yx6ChEWqFZs2ZYtWoVJk6ciD179mDChAl1Pl5WpcD54jJUK5QwMTZET/PmaG4q7ltnSUkJ9u7di+DgYERGRqKqqgqDBg3C0qVLERAQwBJJOiciIgJHjx5FdHQ0y6QG4wqlDikuLoaFhQVeeuklLF++XHQcIq0hSRL8/PyQlZWFs2fP4uI9BbYmFyIu+zoKb5Xj798kDQB0b2cGD7uOmObWHTadWqo9X2lpaZ0SWVlZCTc3NwQFBSEgIADdu3dXewYiEZRKJVxdXdGyZUscPHiQhVKDsVDqkHnz5mHt2rUoKCjgubpE9ZSTkwOXwaPQ9+WvcVnZGkaGBqhV/vu3x/sfH2bdHl9Odka3dmYqzVNaWorw8HAEBwcjIiIClZWVGDhwIAIDAxEQEICePXuq9H5Emmj79u0ICAhAQkIChg0bJjoOPQALpY64evUqLC0t8e6772LJkiWi4xBpnW3HC/Fx2GkolICB4aPPtzMyNICxoQEWT3DEUwMatlJYVlaG8PBwhISEIDw8HBUVFRgwYIBcIi0sLBp0fSJtUltbC2dnZ3Tv3h379u0THYcegnsodcTSpUthamqKuXPnio5CpHXWxuVieXQOACMY1PMF0lqlhFqlhHlhabhZWoU3PGzq9fllZWWIiIhASEgI9u7di4qKCvTv3x+ffvopAgMDWSJJb/3555/IzMzEr7/+KjoKPQKuUOqAoqIiWFtbY9GiRXwZh6ieth0vxLywNJVd7yt/Z0x9yEpleXk5IiMjERwcjL1796K8vBz9+vVDUFAQAgMDYWlpqbI8RNqopqYGvXr1gpOTE3bu3Ck6Dj0CFkod8PLLLyMsLAz5+flo2VL9LwgQ6YqiW+XwXBmPKoVSZdc0NTZEzDsj/mdPZUVFRZ0SWVZWhr59+8ol0srKSmUZiLTdxo0b8dJLLyE1NRXOzs6i49AjYKHUcufOnYO9vT2WLVuG9957T3QcIq3y3KZkJOUXP/Dlm/oyMjTAYEtzbJnlhoqKCuzbtw/BwcHYs2cPysrK0KdPH7lEWltbq+y+RLqiqqoKNjY2GDJkCP7880/RcegRsVBquenTpyMmJgbnzp1Ds2bNRMch0hq510owZlWC2q7f59IeHNixFaWlpXBxcUFgYCACAwO1+jQeosawZs0azJkzBxkZGbCzsxMdhx4RX8rRYhkZGfj999+xdu1alkmietqaXPjQ0UCPS1LWIqPGHB9++CECAwP5Q5HoEZWXl+OLL77A9OnT+XWjZbhCqcWCgoJw7NgxZGdnw9TUVHQcIq0y4ps4XLhVrrbr9zA3Q/xcD7Vdn0gXffPNN/j444+Rk5PDWatahiuUWur06dMICQnBpk2bWCaJ6qm0SoFCNZZJACgsLkdZlULoMY1E2uTevXtYtmwZZs+ezTKpheo5cY00xcKFC2FjY4Pp06eLjkKkdS4Ul0Hdj2YkAOeLy9R8FyLdsWrVKpSVlXH8nZbir85a6OjRo9i7dy/++OMPGBvzXyFRfVWrcEyQJtyHSNvdunULK1aswGuvvYYnn3xSdBx6DFyh1EILFy6Ek5MTpk6dKjoKkVYyMW6cb32NdR8ibbd8+XLU1tZi3rx5oqPQY+LylpY5ePAgYmJiEBYWBkND/rAiehw9zZvDAFDrY2+D/9yHiB7s2rVr+O677zBnzhx07NhRdBx6TCyUWkSSJCxYsAD9+/fHpEmTRMch0kolJSXYu3cvTKqBKpNWartPd3MzvpBD9AiWLVuGJk2aYO7cuaKjUANwiUuLREVF4fDhw1iyZAkMDAxExyHSGnfv3sXvv/+OiRMnokOHDnjmmWdgeC0TBmpaozQyNICHLVdaiB7m4sWLWL9+Pd577z20bdtWdBxqAM6h1BKSJGHAgAFo2rQpEhMTWSiJHuL27dvYvXs3QkJCsH//flRXV8Pd3R0BAQGYMmUKqpu2U+tJOTHvDId1x5Zquz6RLnj11VcREhKCgoICtGzJrxdtxucxWmLXrl04ceIE4uLiWCaJ/kVxcTF27tyJ0NBQxMTEoLa2FkOGDMHXX38Nf39/dOvWrc6fH2bdXm1nebNMEj1Yfn4+Nm7ciKVLl7JM6gCuUGqB2tpa9OnTB0888QT2798vOg6RRrlx4wZ27NiB0NBQxMbGQqlUYvjw4QgICIC/vz+6dOnyr59bdKscnivjUaXC8T6mxoaIeWcEurUzU9k1iXTRjBkzEBUVhXPnzsHMjF8v2o4rlFogODgY6enp+Omnn0RHIdII165dQ1hYGEJDQ3Hw4EEAwMiRI7FmzRpMnjwZTzzxxCNdp1s7Myye4Ih5YWkqy/bZBEeWSaKHyMrKwpYtW/Ddd9+xTOoIrlBqOIVCAQcHB9jZ2WHPnj2i4xAJc/nyZblEJiQkwNDQEKNGjUJgYCAmTZqEDh06PPa118blYnl0ToMzvj/WDq97WDf4OkS67qmnnsKRI0eQk5PD44N1BFcoNdyWLVuQm5uL4OBg0VGIGt3Fixexfft2hIaG4vDhwzAyMoKnpyc2btyIiRMnwtzcXCX3ecPDBu1bmOKT3WehUEr121OprIWpSRN8NsERUwd0V0keIl2WmpqKv/76Cxs3bmSZ1CFcodRgVVVVsLOzw8CBA1koSW9cuHBBLpFHjhxBkyZNMHbsWAQGBmLChAlqHS1SdKscH+1IQ2LeTRgZGjywWN7/eEXBKbzpZo4F77yqtlxEumTixInIyMhARkYGmjRpIjoOqQgLpQZbt24d3nzzTaSnp6NXr16i4xCpTX5+PrZv346QkBAcP34cpqam8Pb2RkBAAMaPH4/WrVs3ap7cayXYmlyIuJzrKCwurzOt0gD/N7Tcw7Yjnh3UHV8veA8hISHIzc1F+/btGzUnkbY5duwY3Nzc8Pvvv2PatGmi45AKsVBqqIqKClhZWcHT0xO//fab6DhEKpeXl4eQkBCEhobi5MmTaNq0KXx9fREQEAA/Pz+0aqW+U2zqo6xKgfPFZahWKGFibIie5s3rnIBz48YN2NjY4KmnnsIPP/wgMCmR5vPy8sKlS5eQmpoKIyMj0XFIhVgoNdS3336LDz/8EFlZWbCyshIdh0glsrOz5RKZmpoKMzMz+Pn5ISAgAL6+vmjRooXoiI/lu+++wzvvvIOTJ0+iT58+ouMQaaSEhASMGDEC27dvh7+/v+g4pGIslBqopKQElpaW8Pf3x4YNG0THIWqQs2fPIjQ0FKGhoUhPT0eLFi0wbtw4BAQEwMfHRydGhtTU1KBPnz4wNzdHfHw8Dx8g+i+SJGHEiBEoLS3FiRMn+DWig/iWtwZavXo17t27hwULFoiOQlRvkiQhLS1NLpGZmZlo2bIlJkyYgCVLlmDs2LFo1qyZ6Jgq1aRJE6xatQpjx45FcHAwpk6dKjoSkUbZv38/EhMTER4ezjKpo7hCqWFu374NCwsLPP/88/juu+9ExyF6JJIk4fTp03KJzMnJQevWrTFx4kQEBgZizJgxejEeZPLkyUhJSUFWVhaaN28uOg6RRpAkCW5ubmjSpAkOHTrEQqmjuEKpYb799ltUV1dj/vz5oqMQPZAkSThx4oRcIs+dO4e2bdti8uTJWLVqFUaPHg0TExPRMRvVihUr4ODggK+++gqfffaZ6DhEGmH37t04fvw4YmNjWSZ1GFcoNciNGzdgaWmJ119/HcuWLRMdh+h/SJKEY8eOySXy/PnzaN++PSZPnoyAgAB4eHjo/Vy5jz/+GN9++y0yMzPRs2dP0XGIhFIqlejTpw86dOiAAwcOiI5DasRCqUHmzp2LH3/8EQUFBSo7AYSooZRKJY4ePSqXyKKiInTs2BH+/v4ICAjAiBEjYGzMhx33lZaWws7ODu7u7ggNDRUdh0iov/76C0899RSSkpLg7u4uOg6pEQulhrh8+TKsrKzw4Ycf4tNPPxUdh/RcbW0tkpKSEBISgu3bt+Py5ct44oknMGXKFAQEBGDYsGGcIfcAf/zxB6ZNm4YDBw5g1KhRouMQCaFQKODk5ARra2vs3btXdBxSMxZKDfHGG2/gjz/+QEFBQaOfCkIE/F+JTEhIQGhoKMLCwnD16lU8+eSTcokcPHgwS+QjkiQJw4YNw927d3Hq1Cmu4JJe+uWXX/DCCy/gxIkT6Nevn+g4pGYslBrg/PnzsLW1xeeff44PP/xQdBzSIwqFAgcPHkRoaCh27NiB69evo1u3bggICEBAQAAGDRoEQ0ND0TG10smTJ+Hq6orVq1fjjTfeEB2HqFFVV1fDzs4O/fv359YPPcFCqQFmzZqFvXv3Ij8/n6NGSO1qamoQGxsrl8ji4mL07NkTgYGBCAgIwIABA/gmpoq8+OKLCA0N5TnfpHd++OEHvPbaa0hLS4Ojo6PoONQIWCgFy8nJgYODA1asWIG3335bdBzSUdXV1YiJiUFoaCh27tyJ27dvw8rKSi6R/fr1Y4lUg+vXr8PW1hZPP/001q9fLzoOUaOoqKiAtbU1PDw88Pvvv4uOQ42EhVKwadOmISEhAbm5uWjatKnoOKRDKisrsX//foSGhmLXrl24e/cubG1t5RLp4uLCEtkIVq1ahffeew8nT56Ei4uL6DhEardy5Uq8//77yMrKgrW1teg41EhYKAVKT09H79698cMPP+Cll14SHYd0QEVFBaKiohAaGordu3ejpKQEvXr1kkukk5MTS2Qjq6mpgYuLCzp06ICDBw/y7590WmlpKSwtLTFx4kT89NNPouNQI2KhFMjf3x+pqanIysrS+2HQ9PjKy8sRGRmJkJAQ7N27F2VlZXB2dpZfrHFwcBAdUe9FR0fDy8sL27Zt4znfpNOWLl2KTz/9FLm5uejevbvoONSIWCgFOXHiBFxdXfHrr79i+vTpouOQliktLUVERARCQkIQERGB8vJy9OnTRy6RdnZ2oiPSf5k0aRJOnjyJrKwsmJmZiY5DpHJ37tyBhYUFnn32WaxZs0Z0HGpkLJSC+Pr6oqCgAOnp6ZztR4/k3r172Lt3L0JDQxEZGYnKykr0799fLpHcq6TZzp07BwcHB8ybNw+LFy8WHYdI5RYtWoTly5fj3Llz6Ny5s+g41MhYKAU4fPgwhg4dir/++gtBQUGi45AGu3PnDvbs2YPQ0FBERUWhqqoKAwcOlEukhYWF6IhUDx999BFWrlzJc75J59y8eRMWFhZ49dVX8fXXX4uOQwKwUDYySZLg4eGBO3fu4OTJkxwaTf/j1q1b2L17N0JDQxEdHY2amhq4u7sjMDAQ/v7+6NGjh+iI9Jh4zjfpqg8++ADr169HQUEBZ67qKZ4H1shiY2MRHx+P3bt3s0yS7ObNm9i1axdCQkJw4MAB1NbWYujQoVi+fDn8/f3RtWtX0RFJBVq0aIGvv/4azz77LOLi4uDh4SE6ElGDXblyBWvXrsXcuXNZJvUYVygbkSRJGDx4MJRKJY4ePcrxIXru+vXr2LlzJ0JCQhAXFwdJkjB8+HAEBARg8uTJ6NKli+iIpAaSJGHo0KG4d+8ez/kmnfDmm29i69atyM/PR5s2bUTHIUH4nawRhYeH4+jRo4iOjmaZ1FNXr17Fjh07EBISgvj4eACAh4cH1q5di8mTJ6NTp06CE5K6GRgYYPXq1RgwYAA2bNiA119/XXQkosd24cIFbNiwAYsXL2aZ1HNcoWwkSqUS/fv3R+vWrREXF8dCqUcuX76M7du3IzQ0FImJiTA0NMTo0aMREBCASZMmoUOHDqIjkgAvvvgitm/fjtzcXJibm4uOQ/RYXnzxRezatQv5+flo0aKF6DgkEAtlIwkNDUVgYCASEhIwbNgw0XFIzYqKiuQSefjwYTRp0gSenp4ICAjAxIkTWSBIPuf7mWeewbp160THIaq3vLw82Nvb45tvvsE777wjOg4JxkLZCGpra+Hs7Izu3btj3759ouOQmpw/f14ukUePHoWJiQnGjh2LgIAATJgwAW3bthUdkTTMypUrMXfuXJ7zTVrp2WefxcGDB5GXl4emTZuKjkOCsVA2gi1btmD69Ok4duwYBgwYIDoOqVB+fj5CQ0MRGhqK48ePw9TUFN7e3ggMDMS4cePQunVr0RFJg9XU1KB3797o2LEjz/kmrXL27Fk4Oztj3bp1eOWVV0THIQ3AQqlmNTU1sLe3R+/evbFjxw7RcUgFcnNz5RJ58uRJNG3aFL6+vggMDISfnx9atmwpOiJpkfvnfPOgA9ImAQEBOHHiBLKzs2FiYiI6DmkAFko1++mnn/Dyyy8jNTUVzs7OouPQY8rKykJoaChCQkJw5swZmJmZYdy4cQgICICPjw83o1ODTJw4EadOneI536QVTp48if79+2Pz5s2YMWOG6DikIVgo1aiyshI2NjYYOnQo/vzzT9FxqB4kSUJGRgZCQkIQGhqKs2fPokWLFhg/fjwCAgLg7e3NH/ykMvfP+Z4/fz4+/fRT0XGIHsjPzw/nzp1Deno656iSjIVSjdasWYM5c+YgIyMDdnZ2ouPQQ0iShLS0NLlEZmVloVWrVpgwYQICAgLg5eXFjeekNvPnz8eqVauQlZXF4zVJYyUlJWHIkCHYtm0bpk6dKjoOaRAWSjUpKyuDlZUVfH198fPPP4uOQ/9CkiScPn1aLpG5ublo06YNJk6ciICAAIwZMwampqaiY5IeuH/O9+DBgxESEiI6DtE/Gj16NG7cuIHTp0/z+GCqg2vVavL999/j1q1bWLRokego9F8kSUJKSor8Yk1+fj7atWuHSZMmYfXq1Rg1ahQ3mVOja9GiBb766is899xzPOebNFJsbCxiY2Oxc+dOlkn6H1yhVIN79+7BwsICU6dO5cBiDSFJEpKTk+USeeHCBbRv3x6TJ09GYGAgRo4ciSZNmoiOSXpOkiQMGTIEpaWlOHnyJPenkca4/9+mQqFAcnIyR1zR/+B3KzVYtWoVysvL8fHHH4uOoteUSiWOHDmC0NBQbN++HUVFRejYsSP8/f0RGBiI4cOH8wc2aZT753wPHDiQ53yTRomMjMSRI0cQFRXFMkn/iCuUKnbr1i1YWFhg9uzZWLFiheg4eqe2thaHDx+WS+Tly5fRuXNnTJkyBQEBARg6dCiMjIxExyR6oNmzZyMsLIznfJNGkCQJ/fv3R4sWLRAfH89CSf+IhVLF5s+fjzVr1iA/Px8dO3YUHUcvKBQKJCYmIjQ0FGFhYbh69SqefPJJBAQEICAgAIMHD+Z+H9Iq169fh42NDZ599ll8//33ouOQngsLC8OUKVMQHx+P4cOHi45DGoqFUoWuXbsGS0tLzJkzB1988YXoODpNoVDg4MGDCAkJwY4dO3Djxg10795dLpFubm4skaTVvv32W7z//vs4deoUevfuLToO6ana2lr07t0bXbt2RVRUlOg4pMFYKFXonXfewebNm1FQUIC2bduKjqNzampqcODAAYSGhmLnzp0oLi5Gz549ERgYiICAAAwYMICPYkhnVFdXw8XFBZ06dUJcXBz/2yYhtm7dimeffRbJyckYOHCg6DikwVgoVeTixYuwtrbGggULsGDBAtFxdEZVVRViYmIQGhqKXbt24fbt27CyspJLZL9+/fiDlnRWVFQUvL29ERwcjMDAQNFxSM/U1NTAwcEBDg4O2LVrl+g4pOFYKFXklVdeQWhoKAoKCtCyZUvRcbRaZWUloqOjERoait27d+Pu3buwtbVFYGAgAgMD0bt3b5ZI0hsTJkxAamoqMjMzedwnNapNmzZh9uzZOH36NFxcXETHIQ3HQqkC+fn5sLOzw9KlSzF37lzRcbRSRUUF9u3bh9DQUOzZswclJSVwcHCQVyIdHR1ZIkkv5eXlwdHRked8U6OqqqqCra0t3N3dsW3bNtFxSAuwUKrAjBkzEBUVhXPnznEFoR7KysoQGRmJ0NBQ7N27F2VlZXB2dkZgYCCmTJkCBwcH0RGJNALP+abGtnbtWrz99ts4e/Ys7O3tRcchLcBC2UCZmZlwcnLC6tWrOYT4EZSWliI8PByhoaGIiIhAeXk5+vTpI5dIOzs70RGJNE5JSQns7OwwdOhQBAcHi45DOq68vBxWVlbw8vLCL7/8IjoOaQkWygaaOnUqjh49ipycHJiamoqOo5Hu3buHvXv3IiQkBPv27UNlZSVcXV0REBCAKVOmwNraWnREIo23ZcsWTJ8+HXFxcRg5cqToOKTDli9fjvnz5yM7OxuWlpai45CWYKFsgNTUVPTp0wcbN27ErFmzRMfRKHfu3MGePXsQEhKCqKgoVFdXw83NTS6RFhYWoiMSaRWlUomhQ4fynG9Sq5KSElhYWCAwMBDr168XHYe0CAtlA0ycOBEZGRnIyMhAkyZNRMcR7tatW9i1axdCQ0Oxf/9+1NTUYPDgwXKJ7N69u+iIRFotJSUFAwcOxNq1a/Haa6+JjkM66PPPP8cXX3yBvLw8dO3aVXQc0iIslI8pOTkZgwYNwtatW/HMM8+IjiPMzZs3sXPnToSGhuLAgQOora3F0KFD5RL55JNPio5IpFNmzZqFnTt3Ijc3F+3atRMdh3TIrVu3YGFhgZkzZ2LlypWi45CWYaF8TGPHjsXly5eRmpoKIyMj0XEa1fXr17Fjxw6EhoYiLi4OkiRh+PDhCAgIgL+/Pzp37iw6IpHOunbtGmxtbfHcc89h7dq1ouOQDvn444+xatUq5Ofno1OnTqLjkJZhoXwM8fHxGDlyJLZv3w5/f3/RcRrF1atXERYWhtDQUMTHx8PAwAAjR45EYGAgJk2axG8+RI2I53yTql2/fh2WlpZ466238OWXX4qOQ1qIhbKeJEnCiBEjUFZWhpSUFJ0etn3p0iW5RCYmJsLQ0BCjR49GYGAgJk6ciA4dOoiOSKSXqqur0bt3b3Tu3BmxsbE6/X2IGse7776LTZs2oaCggFsp6LHwNcF6io6ORmJiIiIiInTym3hRURG2b9+OkJAQJCUloUmTJhgzZgw2bdqECRMmwNzcXHREIr1nYmKCVatWwcfHB9u3b0dAQIDoSKTFLl26hHXr1uGjjz5imaTHxhXKepAkCQMHDoSJiQkOHTqkM4Xy/PnzcolMTk6GiYkJvLy8EBAQgPHjx6Nt27aiIxLRPxg/fjzOnDnDc76pQV577TX89ddfKCgoQKtWrUTHIS3FFcp62L17N1JSUnTiEdO5c+cQGhqK0NBQpKSkwNTUFD4+Pvj9998xbtw4tG7dWnREInqIlStXwtHREd988w0++eQT0XFICxUUFGDjxo1YsmQJyyQ1CFcoH5FSqUSfPn3QoUMHHDhwQHScx5KTkyOXyFOnTqFZs2bw9fVFQEAA/Pz80LJlS9ERiaie5s2bh9WrVyMrK4uzXqneXnjhBURGRiI/P5+r3NQgLJSPaNu2bXj66aeRlJQEd3d30XEeWWZmplwiz5w5AzMzM4wbNw4BAQHw9fVF8+bNRUckoga4f873sGHD8Ndff4mOQ1okOzsbDg4OWLVqFd58803RcUjLsVA+AoVCAUdHR9jY2GDv3r2i4zyQJEk4e/asXCLPnj2LFi1aYPz48QgMDISXlxd/CyXSMb/99huef/55nvNN9fL000/j8OHDyM3Nhampqeg4pOVYKB/BL7/8ghdeeAEnTpxAv379RMf5H5Ik4cyZM3KJzMrKQqtWrTBx4kQEBARg7NixaNq0qeiYRKQmSqUSQ4YMQXl5OU6cOMFzvumhzpw5AxcXF/z444948cUXRcchHcBC+RDV1dWws7ND//79ERoaKjqOTJIknDp1CqGhoQgJCUFeXh7atGmDSZMmISAgAJ6envyNk0iPHD9+HAMHDsS6devw6quvio5DGm7SpElIT09HZmYmmjRpIjoO6QAWyodYv349Xn/9daSnp8PBwUFoFkmSkJKSgpCQEISGhsoDaCdPnoyAgACMGjUKJiYmQjMSkTgzZ87Erl27eM43PdD9Xz62bNmCZ599VnQc0hEslA9QUVEBa2trjBo1Clu2bBGSQalU4tixY3KJLCwsRPv27eHv74+AgACMHDmSv10SEQCe802PxtvbG4WFhUhLS4ORkZHoOKQjuNHmAX744Qdcu3at0ee7KZVKJCUlITQ0FNu3b8fFixfRqVMnuUQOHz6ce6SI6H906tQJixYtwgcffICXX34Zzs7OoiORhklMTERUVBRCQkJYJkml9H6FsqxKgfPFZahWKGFibIie5s3R3NQYpaWlsLS0xMSJE/HTTz+pPUdtbS0OHTokl8grV66gc+fOmDJlCgIDAzFkyBB+8RPRQ90/57tLly44cOCA1h/CQKojSRJGjhyJkpISpKSkwNDQUHQk0iF6ucyVe60EW5MLEZd9HYW3yvH3Rm0AoHs7M7QouYAyo5ZYuHCh2nIoFAokJCQgNDQUYWFhuHbtGp588kkEBQUhMDAQ7u7u/IInonoxMTHBypUr4evri7CwMEyZMkV0JNIQMTExSEhIwN69e/mzhVROr1Yoi26V46MdaUjMuwkjQwPUKv/9H11S1sLA0AjDrNvjy8nO6NZONbMba2pqcPDgQYSGhmLHjh24ceMGevTogYCAAAQEBGDgwIH8QieiBhs/fjzS0tKQmZmJZs2aiY5DgkmShEGDBsHQ0BBJSUlcuSaV05tCue14IT7ZfRYKpfTAIvnfjAwNYGxogMUTHPHUgMc71qy6uhqxsbEIDQ3Fzp07UVxcDAsLCwQGBiIgIACurq784iYilcrNzYWjoyMWLFiARYsWiY5Dgu3ZswcTJkxATEwMRo8eLToO6SC9KJRr43KxPDqnwdeZO9YWb3jYPNKfraqqQkxMDEJCQrBr1y7cuXMH1tbWcons27cvSyQRqdWHH36INWvW8JxvPadUKtG3b1+Ym5sjNjZWdBzSUTpfKLcdL8S8sDSVXe8rf2dM/ZeVysrKSkRHRyMkJAS7d+/GvXv3YGdnJ5fI3r17s0QSUaMpKSmBra0tRowYgW3btomOQ4IEBwdj6tSpOHToEIYMGSI6DukonS6URbfK4bkyHlUKpcquaWpsiJh3Rsh7KisqKhAZGYnQ0FDs2bMHpaWlcHR0lPdEOjo6skQSkTD3z/k+ePAgRowYIToONTKFQgEnJydYWloiIiJCdBzSYTpdKJ/blIyk/OJ67Zl8GCNDA7j1aIPJbS8hNDQU4eHhKCsrQ+/eveUS2atXL5Xdj4ioIZRKJQYPHozKykqcOHGC48f0zP1fKFJSUtC/f3/RcUiH6WyhzL1WgjGrEtR2/Us/vQLn7h3kEmlra6u2exERNcSxY8fg5uaG9evX45VXXhEdhxpJdXU17O3t0bdvX2zfvl10HNJxOlsoP919FluSL6h0dfI+A0iY5NAOK58brPJrExGpwwsvvIDdu3fznG89smHDBrz66qs4c+YMnJycRMchHaezAw/jsq+rpUwCgAQDnLxWpZZrExGpw9KlS1FTU9PoR8mSGJWVlfj888/x9NNPs0xSo9DJQllapUDhrXK13qOwuBxlVQq13oOISFWeeOIJLFq0COvXr0d6erroOKRmGzZswNWrV/Hpp5+KjkJ6QicL5YXiMqj7Ob4E4HxxmZrvQkSkOm+99RasrKzw9ttvQ0d3OxGAsrIyfPnll5gxYwZsbB5tdjJRQ+lkoaxW4ZggTbgPEZEqmJiYYNWqVYiNjUVYWJjoOKQma9aswe3bt7Fw4ULRUUiP6GShNDFunH+sxroPEZGq+Pj4wM/PD++99x4qKipExyEVu3v3Lr7++mu89NJL6NGjh+g4pEd0shH1NG+Oxhgl3qm5Tv71EZGOW7lyJS5fvozly5eLjkIqtnLlSlRUVOCjjz4SHYX0jE42ouamxuj+n5Ns1KXm1mV0e6Ij/Pz88P3336OgoECt9yMiUhUbGxvMmTMHS5cuRVFRkeg4pCLFxcX49ttv8cYbb6BLly6i45Ce4RzKx2BkaAA/25bocfMYwsPDkZiYCIVCgV69esHX1xd+fn4YMmQITExMVH5vIiJVuHfvHuzs7HjOtw758MMPsW7dOuTn56NDhw6i45Ce0dlCqe6TcmLeGQ7rji0B/N835v379yMiIgIRERG4evUqWrZsiTFjxsDPzw8+Pj7o3Lmz2rIQET2OX3/9FTNmzEB8fDyGDx8uOg41wNWrV2FpaYn33nsPn3/+ueg4pId0tlAC6jvLe7ClObbMcvvHjyuVSpw+fVoul0ePHoUkSejbty/8/Pzg6+uLgQMH8jxdIhJOqVTC3d0dVVVVPOdby7399tv47bffUFBQgDZt2oiOQ3pIpwtl0a1yeK6MR5UKx/uYGhsi5p0R6PaIezRv3ryJqKgoREREYN++fbh16xbMzc3h7e0NX19feHl5wdzcXGX5iIjqIzk5GYMGDeI531qssLAQNjY2+OSTT/gyDgmj04USALYdL8S8sDSVXe8rf2dMHdD9sT63trYWycnJiIiIQHh4OE6fPg1DQ0MMGjQIvr6+8PX1RZ8+fWBg0BjvqBMR/Z8XXngBe/bsQW5uLtq2bSs6DtXTSy+9hB07dqCgoAAtWrQQHYf0lM4XSgBYG5eL5dE5Db7O+2Pt8LqHtQoS/Z9Lly4hMjISERER2L9/P0pLS9G5c2f5xR5PT0+0bNlSZfcjIvonV69eha2tLWbMmIHVq1eLjkP1kJeXB3t7e3z99dd49913RcchPaYXhRL4v5XKT3afhUIp1WtPpZGhAYwNDfDZBMfHXpl8FNXV1UhMTJT3XmZlZaFJkyYYNmyYvPfSzs6Oq5dEpBbffPMN5s+fj9OnT8PJyUl0HHpE06dPx4EDB5CXl4dmzZqJjkN6TG8KJfB/eyo/2pGGxLybMDI0eGCxvP/xYdbt8eVk50feM6kq+fn5crmMjY1FVVUVLC0t5UfjI0eO5DcPIlKZ6upqODk5oVu3boiJieEvr1ogIyMDTk5O+P777/Hqq6+KjkN6Tq8K5X2510qwNbkQcTnXUVhcjr//BRgA6G5uBg/bjnh2UHd5NJBI5eXliIuLk/deXrhwAc2aNcOoUaPkgtmzZ0/RMYlIy0VERMDPzw/bt2+Hv7+/6Dj0EIGBgTh+/DhycnI495iE08tC+XdlVQqcLy5DtUIJE2ND9DRvjuamxqJj/StJkpCZmYnw8HBERETg0KFDUCgUcHBwqDNUvUmTJqKjEpEW8vPzQ0ZGBjIyMvgURIOdOnUK/fr1w88//4wXXnhBdBwiFkptd/fu3TpD1a9du4ZWrVrJQ9W9vb05VJ2IHllOTg6cnJywaNEiLFiwQHQc+hfjx49HTk4Ozp49C2NjzV0EIf3BQqlDlEolTp06JZfL5ORkSJKEfv36yS/2DBgwgMOLieiBPvjgA6xduxbZ2dno1q2b6Dj0X44ePQp3d3f8+eefeOqpp0THIQLAQqnTbty4UWeo+u3bt9G+ffs6Q9XbtWsnOiYRaZh79+7B1tYWHh4e+PPPP0XHof/i6emJa9euITU1FYaGhqLjEAFgodQbCoUCycnJ8t7L+9+I3N3d5b2XvXv35pudRAQA+OWXX/DCCy8gISEBw4YNEx2H/iMuLg6jRo3Cjh07MGnSJNFxiGQslHrq4sWLdYaql5WVoUuXLnK5HD16NIeqE+mx++d8V1dXIyUlhVtlNIAkSRg2bBgqKytx/PhxLgCQRmGhJFRVVdUZqp6dnY0mTZpg+PDh8t5LW1tbfvMi0jP3z/n+4Ycf8PLLL4uOo/ciIyPh6+uLyMhIeHt7i45DVAcLJf2PvLw8REZGIjw8HAcPHkRVVRWsrKzqDFVv2rSp6JhE1AhmzJiBvXv38pxvwSRJgqurK8zMzJCQkMBf8EnjsFDSA5WVlSEuLk7ee1lYWIhmzZph9OjRcsHs0aOH6JhEpCZXrlyBra0tZs6cie+++050HL21Y8cO+Pv74+DBgxgxYoToOET/g4WSHpkkScjIyKgzVL22thaOjo7y3svBgwdzqDqRjvn666/x0Ucf8ZxvQWpra+Hi4oLOnTtj//79ouMQ/SMWSnpsd+7cqTNU/fr162jdujXGjh0LX19feHt744knnhAdk4gaqKqqCs7OzujevTv279/Px62N7I8//sC0adNw9OhRuLm5iY5D9I9YKEkllEolTp48KZ83fvz4cUiShP79+9cZqs6ZaUTaKTw8HOPGjUNYWBgmT54sOo7eUCgU6NWrF3r16oXdu3eLjkP0r1goSS2uX7+OqKgohIeHIyoqCnfu3EGHDh3qDFXnBn8i7eLr64vMzEye892Ifv75Z8yaNQunTp1Cnz59RMch+lcslKR2CoUCR48elfdenjlzBoaGhhg8eLC899LZ2ZmP0Yg0XHZ2NpycnPDJJ5/wnO9GUFVVBVtbW7i5uSE4OFh0HKIHYqGkRldUVCQPVY+JiUFZWRmefPLJOkPVW7RoITomEf2D999/H+vWrUN2dja6du0qOo5O+/777/HWW28hPT0dvXr1Eh2H6IFYKEmoqqoqJCQkyHsvc3NzYWJighEjRshjiWxtbUXHJKL/uH/O96hRo/DHH3+IjqOzysvLYWVlhbFjx+LXX38VHYfooVgoSaPk5ubWGapeXV0Na2truVyOGDGCQ9WJBNu8eTNmzpyJxMREDB06VHQcnbRixQrMmzcPWVlZsLKyEh2H6KFYKEljlZWVITY2Vt57WVRUBDMzszpD1bt37y46JpHeUSqVGDRoEGpqanjOtxqUlJTA0tIS/v7+2LBhg+g4RI+EhZK0giRJOHv2rFwuDx8+jNraWjg5Ocljidzd3TlUnaiR3D/ne8OGDXjppZdEx9EpS5Ysweeff468vDx069ZNdByiR8JCSVrp9u3b8lD1yMjIOkPV/fz84O3tjU6dOomOSaTTnn/+eURERCAnJ4djwFTk9u3bsLCwwIwZM7Bq1SrRcYgeGQslaT2lUokTJ07UGaoOAAMGDJAfjbu6unKoOpGK3T/ne9asWSw/KrJgwQKsXLkS586d40ljpFVYKEnnXLt2rc5Q9bt376JDhw7w8fGBr68vxo4dy9UUIhW5f853amoqHB0dRcfRajdu3ICFhQXeeOMNLFu2THQconphoSSdplAocOTIEXnvZVpaGoyMjOoMVXdycuJQdaLHVFVVBScnJ/Ts2RPR0dH8WmqA9957Dxs3bkR+fj7Mzc1FxyGqFxZK0itFRUWIiIiQh6qXl5ejW7du8qPxUaNGcag6UT3t3bsX48ePx44dOzBp0iTRcbTSpUuXYG1tjXnz5uGTTz4RHYeo3lgoSW9VVlbWGaqel5cHExMTjBw5Ui6YNjY2omMSaTxJkuDr64vs7GxkZGRwVuxjeP3117Ft2zbk5+ejdevWouMQ1RsLJdF/5ObmyuUyPj4e1dXVsLGxqTNU3dTUVHRMIo10/5zvTz/9FB9//LHoOFrl/PnzsLW1xeeff44PP/xQdByix8JCSfQPSktL6wxVv3jxIpo3b47Ro0fDz88PPj4+nA9H9F/mzp2L9evX85zvepo5cyYiIiJw7tw5NG/eXHQcosfCQkn0EJIkIS0tTd57mZSUhNraWjg7O9cZqm5sbCw6KpFQd+/eha2tLTw9PbF161bRcbRCTk4OevXqhZUrV+Ktt94SHYfosbFQEtXT7du3ER0dLQ9Vv3HjBtq0aQMvLy/4+vrC29sbHTt2FB2TSAie810/zzzzDBITE5Gbm8u9p6TVWCiJGkCpVCIlJUXee5mSkgIDAwN5qLqfnx/69evHoeqkN+6f861QKHD8+HGe8/0AaWlpcHFxwQ8//MDjK0nrsVASqdC1a9cQGRmJiIgIREVF4d69e+jYsSN8fHzg5+eHMWPGoE2bNqJjEqnV0aNH4e7ujh9//BEvvvii6Dgaa/LkyThz5gyysrLQpEkT0XGIGoSFkkhNampqkJSUJO+9TE9Ph5GREYYMGSLvvXR0dOQgaNJJ06dPR2RkJHJzc/lL1D9ISUnBgAED8Ntvv+G5554THYeowVgoiRrJhQsX5NXLAwcOoLy8HN27d68zVJ1veJKuuHz5Muzs7HjO97/w8fHB+fPn5V80ibQdCyWRAJWVlYiPj0d4eDjCw8ORn58PU1NTeai6n58frKysRMckapCvvvoKH3/8Mc6cOQMHBwfRcTTGoUOHMGzYMAQHByMwMFB0HCKVYKEkEkySJOTm5sozL+Pj41FTUwNbW1u5XA4bNoxD1Unr3D/n28LCAlFRUdzegf/7evfw8MCdO3dw8uRJvrBHOoOFkkjDlJSU4MCBA/Ley0uXLqF58+bw9PSUh6pzaDRpiz179mDChAnYuXMnJk6cKDqOcDExMRgzZgz27NmDcePGiY5DpDIslEQaTJIknDlzRh5LdOTIESiVSvTu3Vt+sWfQoEEcqk4ai+d8/3+SJMHd3R0AcOTIEa7Ykk5hoSTSIrdu3aozVP3mzZto27ZtnaHqHTp0EB2TqI6srCw4Oztj8eLF+Oijj0THEWbv3r0YP3489u/fD09PT9FxiFSKhZJIS9XW1iIlJUXee3nixAkYGBhg4MCB8t7Lvn37co8WaYT33nsPP/zwg96e861UKtGvXz+0adMGcXFxXJ0kncNCSaQjrl69Ko8lio6Oxr1799CpU6c6Q9Vbt24tOibpKX0/5zskJARBQUE8kpJ0FgslkQ6qqanB4cOH5Rd7zp49C2Nj4zpD1R0cHLhKQo3q559/xqxZs3Do0CEMGTJEdJxGU1tbCycnJ/Ts2RORkZGi4xCpBQslkR44f/48IiMjER4ejtjYWFRUVKBHjx51hqqbmZmJjkk6TqlUws3NDUqlEseOHdObgd5btmzB9OnTcfz4cbi6uoqOQ6QWLJREeqaioqLOUPWCggKYmprCw8ND3ntpaWkpOibpqCNHjmDw4MH46aefMHv2bNFx1K6mpgb29vbo3bs3duzYIToOkdqwUBLpMUmSkJOTI7/Yk5CQgJqaGtjZ2dUZqm5iYiI6KukQfTrn+8cff8Qrr7yC1NRUODs7i45DpDYslEQkKykpQUxMjLz38vLly2jRokWdoepPPvmk6Jik5S5fvgxbW1u8+OKLWLlypeg4alNZWQkbGxsMGzYMf/zxh+g4RGrFQklE/0iSJKSmpsrl8v5QdRcXlzpD1fVlHxyp1rJly7Bw4UKcOXMGvXr1Eh1HLVavXo133nkHmZmZsLW1FR2HSK1YKInokRQXFyM6Ohrh4eHYt28fiouL0a5duzpD1du3by86JmmJqqoqODo6wsrKCvv27dO5iQNlZWWwtLSEn58ffv75Z9FxiNSOhZKI6q22thbHjx+X916ePHkSBgYGcHNzk/de9unTh0PV6YF0+Zzvr776CgsXLkROTg569uwpOg6R2rFQElGDXblypc5Q9ZKSEjzxxBPyUHVPT08OVaf/IUkSfHx8kJubi7Nnz+rMOd93796FhYUFnn76aXz//fei4xA1ChZKIlKp6upqeah6eHg4MjMzYWxsjGHDhslzL3v16qVzjzjp8dw/5/uzzz7D/PnzRcdRicWLF2PZsmU4d+4cunTpIjoOUaNgoSQitSooKKgzVL2yshI9e/aUy6WHhweHquu5d999Fz/++COys7O1fopAcXExLCws8OKLL2LFihWi4xA1GhZKImo0FRUVOHjwoDxU/fz582jatGmdoeoWFhaiY1Iju3/O95gxY/D777+LjtMg8+bNw9q1a5Gfn4+OHTuKjkPUaFgoiUgISZKQnZ1dZ6i6QqGAvb29PJZo6NChHKquJzZt2oTZs2fj8OHDGDx4sOg4j+Xq1auwtLTEu+++iyVLloiOQ9SoWCiJSCPcu3evzlD1K1euoGXLlhgzZgx8fX3h4+PD/Wg6TKlUYuDAgQCA5ORkrZxvOmfOHPzyyy8oKChA27ZtRcchalQslESkcSRJwunTp+UXe44ePQpJktC3b19576Wbm5tWlg76d0lJSRgyZIhWnvNdVFQEa2trLFy4EAsWLBAdh6jRsVASkca7efNmnaHqt27dQrt27eDt7Q0/Pz94eXnB3NxcdExSgeeeew5RUVHIycnRqnO+X375ZYSFhSE/Px8tW7YUHYeo0bFQEpFWqa2tRXJysvxo/NSpUzA0NISbm5u897JPnz4cS6Sl7p/z/dJLL+Hbb78VHeeRnDt3Dvb29li2bBnee+890XGIhGChJCKtdvnyZXmo+v79+1FSUoLOnTvLj8Y9PT3RqlUr0TGpHpYuXYpFixZpzTnfzz//PPbv34+8vDyOwCK9xUJJRDqjuroahw4dkvdeZmVloUmTJnWGqtvb23P1UsNVVlbCyclJK875zszMhJOTE1avXo3XX39ddBwiYVgoiUhn5efny0PV4+LiUFlZCQsLizpD1Zs1ayY6Jv2D3bt3Y+LEidi1axcmTJggOs6/CgoKQnJyMnJycmBqaio6DpEwLJREpBfKy8sRFxcnr15euHABTZs2xahRo+S9lz179hQdk/5DkiR4e3vj3LlzOHv2rEaWtdOnT6Nv377YtGkTZs6cKToOkVAslESkdyRJQmZmpvxiT2JiIhQKBXr16iWXyyFDhnCoumCZmZno3bs3Pv/8c8ybN090nP8xYcIEZGZmyufVE+kzFkoi0nt3796tM1T96tWraNmyJcaOHSsPVe/cubPomHpJU8/5Tk5OxqBBg7B161Y888wzouMQCcdCSUT0N0qlss5Q9eTkZEiShH79+snnjQ8YMIBD1RvJnTt3YGtrCy8vL2zZskV0HNmYMWNw5coVpKam8r8FIrBQEhE90M2bN7Fv3z5ERERg3759uH37NszNzeWh6mPHjuVQdTXbuHEjXnzxRY055/vgwYPw8PBAWFgYJk+eLDoOkUZgoSQiekQKhaLOUPXTp0/D0NAQgwYNkvdeuri4aPSYG21UW1uLgQMHwsDAAMeOHYOhoaGwLJIkYfjw4SgvL0dKSgr/XRP9BwslEdFjunTpUp2h6qWlpejSpUudoeo8hk81Dh8+jKFDh2Ljxo2YNWuWsBxRUVHw9vZGREQEfHx8hOUg0jQslEREKlBVVVVnqHp2djaaNGmC4cOHywXTzs6OK1oN8Oyzz2L//v3IyclB69atG/3+kiRh4MCBMDExwaFDh/jvkuhvWCiJiNTg3LlzdYaqV1VVwdLSUn6xZ8SIERyqXk+XLl2CnZ0dXn75ZaxYsaLR779z505MnjwZcXFxGDlyZKPfn0iTsVASEalZeXk5YmNj5dXLwsJCNGvWrM5Q9R49eoiOqRW+/PJLfPLJJ41+zrdSqYSLiws6deqEmJiYRrsvkbZgoSQiakSSJCEjI0N+sefQoUNQKBRwcHCoM1S9SZMmoqNqpMrKSjg6OsLGxgaRkZGN9th527ZtePrpp5GUlAR3d/dGuSeRNmGhJCIS6O7du9i/f79cMK9du4ZWrVrVGar+xBNPiI6pUXbt2oVJkyZh9+7dGD9+vNrvp1Ao5BK7d+9etd+PSBuxUBIRaQilUolTp04hPDwcEREROHbsGCRJQv/+/eW9l66urno/SFuSJHh5eSE/P79RzvnevHkzZs6ciZMnT6Jv375qvReRtmKhJCLSUDdu3KgzVP3OnTto3759naHq7dq1Ex1TiMY657u6uhq2trYYMGAAQkJC1HYfIm3HQklEpAUUCgWOHj0qPxpPTU2FoaEh3N3d5b2XvXv31qtRNu+88w5++ukn5OTkoEuXLmq5x/r16/H6668jPT0dDg4OarkHkS5goSQi0kIXL16sM1S9rKwMTz75ZJ2h6i1atBAdU63un/Pt7e2N3377TeXXr6iogLW1NUaNGqVR54gTaSIWSiIiLVdVVYXExER572VOTg5MTEzkoep+fn6wsbHRydXLn376CS+99JJa3r7+9ttv8cEHHyA7OxtWVlYqvTaRrmGhJCLSMXl5efKj8YMHD6KqqgpWVlZ1hqo3bdpUdEyVuH/Ot6GhIZKTk1V2zndpaSksLCwwefJk/Pjjjyq5JpEuY6EkItJhZWVldYaqFxUVwczMrM5Q9e7du4uO2SD3z/netGkTZs6cqZJrfvnll1i8eDFyc3O1/u+HqDGwUBIR6QlJknD27Fm5XB4+fBi1tbVwdHSUy+XgwYO1cqj6tGnTEBMTo5Jzvu/cuQMLCws899xzWL16tYoSEuk2FkoiIj11584d7N+/H+Hh4YiMjMT169fRunXrOkPVO3XqJDrmI7l48SLs7OzwyiuvNPic74ULF2LFihXIz8/nUHmiR8RCSUREUCqVOHnypPxiz/HjxyFJElxdXesMVVfVHkV1uH/Od1paGuzt7R/rGjdu3IClpSVee+01fPXVVypOSKS7WCiJiOh/XL9+XR6qHhUVhTt37qBDhw51hqq3bdtWdMw67p/zbWtri4iIiMd6q/3999/Hhg0bUFBQAHNzczWkJNJNLJRERPRACoUCR44ckfdepqWlwcjIqM5QdWdnZ40YS7Rz505MnjwZe/bswbhx4+r1uZcvX4aVlRU++OADLF68WE0JiXQTCyUREdVLUVERIiMjER4ejpiYGJSXl6Nr167yUPXRo0cLG6ouSRLGjh2LgoKCfzznu6xKgfPFZahWKGFibIie5s3R3NQYAPDGG2/gjz/+QEFBQYNf7CHSNyyURET02KqqqpCQkCDvvczNzYWJiQlGjBhRZ6h6Y8rIyEDv3r3xxRdf4MMPP0TutRJsTS5EXPZ1FN4qx99/6BkA6N7ODP27NMUP703DojkvqfVscCJdxUJJREQqk5ubW2eoenV1NaytreVyOXz48EYZqj5nzhxsDt4NrwW/4FhRCYwMDVCr/PcfdwaSEpKBIQZbtMVXAX3QrZ2Z2jMS6RIWSiIiUovS0tI6Q9UvXrwIMzMzeHp6ymOJ1DU0/Of4LCzemwkDI2PA4NHfTDcyNICxoQEWT3DEUwM40JzoUbFQEhGR2kmShPT0dLlcJiUloba2Fs7OzvLey8GDB8PY2LjB91obl4vl0TkAJPzfQ+3HM3esLd7waNzH9UTaioWSiIga3e3bt+sMVb9x4wZat24NLy8vefWyY8eO9b7utuOFmBeWprKcX/k7YypXKokeioWSiIiEUiqVOHHiRJ2h6gYGBnB1dZXHEvXv3/+hQ9WLbpXDc2U8qhRKlWUzNTZEzDsjuKeS6CFYKImISKNcu3atzlD1u3fvokOHDvDx8ZGHqrdp0+Z/Pu+5TclIyi9+4Ms39WVkaIDBlubYMstNZdck0kUslEREpLFqamrqDFVPT0+HkZERhgwZIu+9dHJyQt71UoxZlaC2HDHvDId1x5Zquz6RtmOhJCIirVFYWCgPVT9w4ADKy8vRrVs39JzyPi42tYCyAS/h/BsjQwM859YDn05wVPm1iXQFCyUREWmlyspKeaj6LmVfoEUHtd2rh7kZ4ud6qO36RNqOhZKIiLRaaZUCzp9GQZ0/zAwApH/qJR/TSER1Pfq0VyIiIg10obhMrWUS+L+JlueLy9R8FyLtxUJJRERarVqFY4I04T5E2oiFkoiItJqJceP8KGus+xBpI351EBGRVutp3lwN73bXZfCf+xDRP2OhJCIirdbc1Bjd1XySTXdzM76QQ/QALJRERKT1POw6wshQPeuURoYG8LCt/7niRPqEhZKIiLTeNLfuKj1y8e9qlRKeHdRdLdcm0hUslEREpPVsOrXEMOv2Kl+lNDI0wDDr9jx2keghWCiJiEgnfDnZGcYqLpTGhgb4crKzSq9JpItYKImISCd0a2eGxSo+b/uzCY7opuYXfoh0AQslERHpjKcGdMfcsbYqudb7Y+0wdQD3ThI9Cp7lTUREOmfb8UJ8svssFEqpXi/rGBkawNjQAJ9NcGSZJKoHFkoiItJJRbfK8dGONCTm3YSRocEDi+X9jw+zbo8vJzvzMTdRPbFQEhGRTsu9VoKtyYWIy7mOwuJy/P2HngH+b2i5h21HPDuoO9/mJnpMLJRERKQ3yqoUOF9chmqFEibGhuhp3pwn4BCpAAslERERETUI3/ImIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogZhoSQiIiKiBmGhJCIiIqIGYaEkIiIiogb5f8O4IKqSOPjrAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "graphset_1 = gs.array(\n",
    "    [\n",
    "        nx.to_numpy_array(nx.erdos_renyi_graph(n=5, p=0.6, directed=True))\n",
    "        for i in range(10)\n",
    "    ]\n",
    ")\n",
    "graphset_2 = gs.array(\n",
    "    [\n",
    "        nx.to_numpy_array(nx.erdos_renyi_graph(n=5, p=0.6, directed=True))\n",
    "        for i in range(100)\n",
    "    ]\n",
    ")\n",
    "graphset_3 = gs.array(\n",
    "    [\n",
    "        nx.to_numpy_array(nx.erdos_renyi_graph(n=3, p=0.6, directed=True))\n",
    "        for i in range(1000)\n",
    "    ]\n",
    ")\n",
    "\n",
    "nx.draw(nx.from_numpy_array(graphset_1[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d766dd07",
   "metadata": {},
   "source": [
    "### A primer in space, metric and aligners"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0e65cec",
   "metadata": {},
   "source": [
    "The first step is to create the total space and then add quotient structure to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bade9271",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_space = GraphSpace(n_nodes=5)\n",
    "total_space.equip_with_group_action()  # permutations by default\n",
    "\n",
    "graph_space = total_space.equip_with_quotient_structure()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36336c04",
   "metadata": {},
   "source": [
    "By default, the total space comes equipped with the Frobenius metric (`MatricesMetric`) and graph space with a quotient metric."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8656a77",
   "metadata": {},
   "source": [
    "With the FAQ alignment and the default Frobenius norm on the total space, we match two graphs and a set of graphs to a base graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4074d7e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "permutated_graph = total_space.aligner.align(graphset_1[1], graphset_1[0])\n",
    "\n",
    "permuted_graphs = total_space.aligner.align(graphset_1[1:3], graphset_1[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d7bfb1e",
   "metadata": {},
   "source": [
    "To compute the distance we can either call the distance function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "abe70991",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.6457513110645907"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "graph_space.metric.dist(graphset_1[0], graphset_1[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "440de4b1",
   "metadata": {},
   "source": [
    "Or, if matching has been already done, we can use the total space distance, to avoid computing the matching twice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8fb68953",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.6457513110645907"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_space.metric.dist(graphset_1[0], permutated_graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aefc53f6",
   "metadata": {},
   "source": [
    "We can also align points to geodesics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "488ff262",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init_point, end_point = graph_space.random_point(2)\n",
    "\n",
    "geodesic_func = graph_space.metric.geodesic(init_point, end_point)\n",
    "\n",
    "aligned_init_point = total_space.aligner.align_point_to_geodesic(\n",
    "    geodesic_func, init_point\n",
    ")\n",
    "\n",
    "total_space.metric.dist(init_point, aligned_init_point)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18fe3623",
   "metadata": {},
   "source": [
    "This short introduction should be enough to set you up for experimenting with the learning algorithms on graphs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dc0c7a5",
   "metadata": {},
   "source": [
    "### Frechet Mean\n",
    "Reference: Calissano, A., Feragen, A., & Vantini, S. (2020). Populations of unlabeled networks: Graph space geometry and geodesic principal components. MOX Report.\n",
    "\n",
    "Given $\\{[X_1], \\dots, [X_k]\\}, [x_i] \\in X/T$, we estimate the Frechet Mean using AAC consisting on two steps:\n",
    "1. Compute $\\hat{X}$ as arithmetic mean of $\\{X_1, \\dots, X_k\\}, X_i \\in X$ \n",
    "2. Using graph to graph alignment to find $\\{X_1, \\dots, X_k\\}, X_i \\in X$ optimally aligned with $\\hat{X}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfbe01a3",
   "metadata": {},
   "source": [
    "Let's instantiate the graph space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e1b4b3d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_space = GraphSpace(n_nodes=5)\n",
    "total_space.equip_with_group_action()\n",
    "total_space.equip_with_quotient_structure();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa8ede70",
   "metadata": {},
   "source": [
    "And now create the estimator, and fit the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5073938a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Maximum number of iterations 20 reached. The estimate may be inaccurate\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0.  , 0.9 , 0.37, 0.9 , 0.59],\n",
       "       [0.61, 0.  , 0.25, 0.83, 0.66],\n",
       "       [0.29, 0.79, 0.  , 0.44, 0.49],\n",
       "       [0.19, 0.61, 0.23, 0.  , 0.29],\n",
       "       [0.87, 0.92, 0.73, 0.86, 0.  ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aac_fm = AAC(space=total_space, estimate=\"frechet_mean\", max_iter=20)\n",
    "\n",
    "fm = aac_fm.fit(graphset_2)\n",
    "\n",
    "fm.estimate_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83e636a5",
   "metadata": {},
   "source": [
    "### Principal Components\n",
    "Reference: Calissano, A., Feragen, A., & Vantini, S. (2020). Populations of unlabeled networks: Graph space geometry and geodesic principal components. MOX Report.\n",
    "\n",
    "We estimate the Generalized Geodesics Principal Components Analysis (GGPCA) using AAC. Given $\\{[X_1], \\dots, [X_k]\\}, (s_i,[X_i]) \\in X/T $ we are searching for:\n",
    "$\\gamma: \\mathbb{R}\\rightarrow X/T$ generalized geodesic principal component capturing the majority of the variability of the dataset. The AAC for ggpca works in two steps: \n",
    "\n",
    "1. finding $\\delta: \\mathbb{R}\\rightarrow X$ principal component in the set of adjecency matrices $\\{X_1, \\dots, X_k\\}, X_i \\in X$ \n",
    "2. finding $\\{X_1, \\dots, X_k\\}, X_i \\in X$ as optimally aligned with respect to $\\gamma$. The estimation required a point to geodesic aligment defined in the metric."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb886a18",
   "metadata": {},
   "source": [
    "As before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c1eae258",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_space = GraphSpace(n_nodes=3)\n",
    "total_space.equip_with_group_action()\n",
    "total_space.equip_with_quotient_structure();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d2c0141",
   "metadata": {},
   "source": [
    "For GGPCA, we also need the point to geodesic aligner."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9754a5aa",
   "metadata": {},
   "source": [
    "Again, create the estimator and fit the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "434fef57",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "aac_ggpca = AAC(space=total_space, estimate=\"ggpca\", n_components=2)\n",
    "\n",
    "aac_ggpca.fit(graphset_3);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c85724c",
   "metadata": {},
   "source": [
    "## Regression\n",
    "Reference: Calissano, A., Feragen, A., & Vantini, S. (2022). Graph-valued regression: Prediction of unlabelled networks in a non-Euclidean graph space. Journal of Multivariate Analysis, 190, 104950.\n",
    "\n",
    "We estimate a graph-to-value regression model to predict graph from scalar or vectors. Given $\\{(s_1,[X_1]), \\dots, (s_k, [X_k])\\}, (s_i,[X_i]) \\in \\mathbb{R}^p\\times X/T $ we are searching for:\n",
    "$$f: \\mathbb{R}^p\\rightarrow X/T$$\n",
    "where $f\\in \\mathcal{F}(X/T)$ is a generalized geodesic regression model, i.e., the canonical projection onto Graph Space of a regression line $h_\\beta : \\mathbb{R}^p\\rightarrow X$ of the form $$h_\\beta(s) = \\sum_{j=1}^{p} \\beta_i s_i$$\n",
    "The AAC algorithm for regression combines the estimation of $h_\\beta$ given $\\{X_1, \\dots, X_k\\}, X_i \\in X$\n",
    "$$\\sum_{i=0}^{k} d_X(h_\\beta(s_i), X_i)$$\n",
    "and the searching for $\\{X_1, \\dots, X_k\\}, X_i \\in X$ optimally aligned with respect to the prediction along the current regression model:\n",
    "$$\\min_{t\\in T}d_X(h_\\beta(s_i),t^TX_it)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "893da39b",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_space = GraphSpace(n_nodes=5)\n",
    "total_space.equip_with_group_action()\n",
    "total_space.equip_with_quotient_structure();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6f33d152",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = gs.array([random.randint(0, 10) for i in range(10)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a4cb1e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "aac_reg = AAC(space=total_space, estimate=\"regression\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0a2152ae",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Maximum number of iterations 20 reached. The estimate may be inaccurate\n"
     ]
    }
   ],
   "source": [
    "aac_reg.fit(s, graphset_1);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46835c03",
   "metadata": {},
   "source": [
    "The coefficients are saved in the following attributes and they can be changed into a graph shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9b71203c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.        ],\n",
       "       [ 0.02115813],\n",
       "       [ 0.00890869],\n",
       "       [-0.0233853 ],\n",
       "       [-0.00668151],\n",
       "       [-0.02561247],\n",
       "       [-0.        ],\n",
       "       [ 0.04899777],\n",
       "       [-0.04120267],\n",
       "       [ 0.02004454],\n",
       "       [ 0.0233853 ],\n",
       "       [ 0.06458797],\n",
       "       [-0.        ],\n",
       "       [-0.06904232],\n",
       "       [-0.02115813],\n",
       "       [-0.02115813],\n",
       "       [ 0.07238307],\n",
       "       [ 0.04008909],\n",
       "       [-0.        ],\n",
       "       [ 0.07572383],\n",
       "       [-0.        ],\n",
       "       [-0.04788419],\n",
       "       [ 0.0701559 ],\n",
       "       [ 0.00445434],\n",
       "       [-0.        ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aac_reg.total_space_estimator.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "874fefa2",
   "metadata": {},
   "source": [
    "A graph can be predicted using the fit model and the corresponding prediction error can be computed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b03476b8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.93063035236635"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "graph_pred = aac_reg.total_space_estimator.predict(s)\n",
    "\n",
    "gs.sum(graph_space.metric.dist(graphset_1, graph_pred))"
   ]
  }
 ],
 "metadata": {
  "backends": [
   "numpy"
  ],
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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
}