{
 "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": [],
   "source": [
    "import random\n",
    "\n",
    "import networkx as nx\n",
    "\n",
    "import geomstats.backend as gs\n",
    "from geomstats.learning.frechet_mean import FrechetMean\n",
    "from geomstats.learning.pca import GGPCA\n",
    "from geomstats.learning.regression import GeneralizedGeodesicRegression\n",
    "from geomstats.metric_geometry.graph_space import GraphSpace\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": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAemlJREFUeJzt3WdYVFf7NfA1gKAidrFSVEAUEQwg2HtL7Dj22EvsIhiNRk1MHmMMTcXYFRsWpNi7iF1EI2KnCINdwIJ0mHk/PH9940OKygx7yvpdV658QM9ZJAqLfe99jkShUChARERERPSZ9EQHICIiIiLNxkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFYiA6ABGpn8zcAiSlZSKvQA5DAz1YVjGGsRG/XBAR0V/jdwgiAgDEPcvA9ssyRNx7Dll6FhR/+pgEgHnlsmjfwBRDXc1hXd1EVEwiIlJDEoVCofj3X0ZE2iolPQtzw2JxNj4V+noSFMr//kvCu4+3tqqKxX3tYVa5bAkmJSIidcVCSaTDdl6RYeG+WyiQK/6xSP4vfT0JDPQk+LGXHQa5mKswIRERaQIWSiIdFRARB+9j94t9Ha8uNpjS3loJiYiISFPxlDeRDtp5RaaUMgkA3sfuY9cVmVKuRUREmomFkkjHpKRnYeG+W0q95oJ9t5CSnqXUaxIRkeZgoSTSMXPDYlHwCfslP0aBXIG5YbFKvSYREWkOFkoiHRL3LANn41M/6QDOxyiUK3A2PhXxzzOUel0iItIMLJREOmT7ZRn09SQquba+ngTbLnEvJRGRLmKhJNIhEfeeK3118p1CuQIR95+r5NpERKTeWCiJdMTb3ALIVHxwRpaWhczcApXeg4iI1A8LJZGOSE7LhKofOqsAkJSWqeK7EBGRumGhJNIReQVyrboPERGpDxZKIh1haFAyf91L6j5ERKQ++JWfSEdYVjGGas53/3+S/7sPERHpFhZKIh1hbGQA88plVXoP8yplYWxkoNJ7EBGR+mGhJNIh7RuYquw5lAp5IQxT4/Do0SOVXJ+IiNQXCyWRDhnqaq6y51BK9PQRvdMPlpaWGD58OGJiYlRyHyIiUj8slEQ6xLq6CVpbVVX6KqW+ngStraoi+cZlLF26FJGRkXB0dETnzp1x5MgRKBSqfmARERGJxEJJpGMW97WHgZILpYGeBIv72qN8+fLw8PBAQkICduzYgVevXqF79+6wt7fHpk2bkJubq9T7EhGRemChJNIxZpXL4sdedkq95qJedjD704EfAwMDDBo0CFFRUYiMjES9evUwevRoWFpaYvHixUhPT1fq/YmISCyJgrMoIp2jUCjQZuJipFRyLPa1ZnVpgMntrf711929exd+fn7YsmUL9PT0MHr0aMyYMQP169cvdgYiIhKLK5REOsjHxwfn1nwPqUUejAz0PnlPpb6eBEYGevi1n/1HlUkAsLW1xZo1ayCTyfDtt99i586dsLa2hru7Oy5evPg5nwYREakJrlAS6Zjjx4+jW7dumD17NhYvXoyU9CzMDYvF2fhU6OtJ/vkUuLwQ0NNHa6uqWNzX/oMx96fKzs7Gli1b4Ovri/v376N58+bw8vJC7969oa+v/9nXJSKiksdCSaRDEhMT4ezsDDc3N+zfv/+D4hb3LAPbL8sQcf85ZGlZ+PMXBgmAyoZyPLhwACFLZqC9U0OlZZLL5Th48CC8vb1x5swZ1K9fHzNmzMCoUaNgbMy37hARaQIWSiIdkZmZiebNmyM7OxtRUVGoVKnS3//a3AIkpWUir0AOQwM9WFYxhiI/B9WqVcOPP/6Ib7/9ViUZo6Oj4ePjg+DgYJQvXx4TJ07ElClTULNmTZXcj4iIlIOFkkgHKBQKDBo0CIcOHcKlS5dgZ/d5p7zd3d0hk8lw5coVJSf8UHJyMpYtW4Z169YhLy8PQ4cOhaen52fnJiIi1eKhHCIdsHTpUuzevRubN28uVimTSqWIjo7GgwcPlJiuKAsLC/j6+iIlJQU///wzjh07hsaNG6N79+44ceIEH5RORKRmWCiJtNzRo0fx3XffYd68eejXr1+xrtWjRw+ULl0ae/bsUVK6f1axYkXMmjULiYmJ2Lp1K548eYLOnTujadOm2Lp1K/Ly8kokBxER/TOOvIm0WHx8PFxcXNCyZUvs3btXKaen+/Xrh4cPHyIqKkoJCT+NQqHAqVOn4OPjg8OHD6N27dqYNm0axo8fj4oVK5Z4HiIi+i8WSiIt9fbtW7i5uSEvLw9RUVFKK1w7duzAkCFD8ODBA1haWirlmp/j1q1b8PX1xbZt22BoaIixY8di+vTpQjMREekqjryJtJBCocCoUaOQnJyM8PBwpa7e9ejRA0ZGRiU29v47dnZ22LBhA5KTkzF9+nRs2bIFVlZWGDRokMoPDRER0YdYKIm00JIlS7Bnzx5s3boVjRo1Uuq1TUxM0L17dwQHByv1up+rRo0a+PnnnyGTybBs2TJER0ejWbNmaNu2Lfbt2we5XC46IhGR1mOhJNIyhw8fxrx587BgwQL06dNHJfeQSqWIiopCcnKySq7/OYyNjTF58mTcu3cPISEhKCgoQO/evdGwYUOsWbMG2dnZoiMSEWkt7qEk0iJxcXFwcXFBmzZtEB4eDj091fzM+ObNG5iamuI///kPPD09VXIPZbh48SJ8fHwQFhaGypUrY/LkyZg0aRJMTU1FRyMi0ioslERaIiMjA25ubigsLMTly5dRoUIFld6vd+/eePbsGS5duqTS+yhDQkIC/P39sXHjRhQWFmLEiBGYOXMmGjRoIDoaEZFW4MibSAsoFAqMHDkSKSkpCA8PV3mZBP479r58+TJkMpnK71Vc9evXx4oVK5CSkoIFCxZg3759sLW1Ra9evRAZGckHpRMRFRMLJZEWWLx4MUJDQ7Ft2zbY2tqWyD179uwJQ0NDhISElMj9lKFy5cqYO3cukpKSsGnTJiQmJqJdu3ZwcXHBzp07UVBQIDoiEZFGYqEk0nAHDx7E/Pnz8cMPP6BXr14ldt8KFSqga9euanPa+1MYGRlh5MiRiI2NxZEjR1CpUiUMHjwY9evXh5+fH968eSM6IhGRRuEeSiINdv/+fbi4uKB9+/YIDQ1V2SGcv7N161YMHz4cMpkMZmZmJXpvZYuJiYGPjw927NiBsmXLYvz48Zg2bZrGf15ERCWBhZJIQ7158wZubm5QKBS4fPkyypcvX+IZXr9+DVNTU/z666+YMWNGid9fFR49eoQVK1Zg9erVyMzMxMCBA+Hp6YmmTZuKjkZEpLZYKIk0kFwuh7u7O06dOoWoqCihp5V79uyJ9PR0nD9/XlgGVcjIyMDGjRvh5+eH5ORkdOjQAV5eXujWrRskEonoeEREaoV7KIk00M8//4y9e/di+/btwh99I5VKceHCBTx8+FBoDmUzMTHB9OnTER8fj127diEjIwNffvklGjdujA0bNiAnJ0d0RCIitcFCSaRh9u/fj4ULF+LHH39Ejx49RMdBr169UKpUKY067f0pDAwMMGDAAFy+fBlnzpyBtbU1xo0bB0tLS/z8889IS0sTHZGISDiOvIk0yN27d+Hq6oqOHTtiz549JX4I5+/06NEDr169wrlz50RHKRH379+Hn58fAgMDIZFIMGrUKHh4eMDKykp0NCIiIVgoiTTE69ev4erqCn19fVy6dAkmJiaiI723efNmjBw5Eg8fPkTt2rVFxykxL168wKpVqxAQEIDU1FT06dMHnp6eaNGiBfdZEpFOUY/lDSL6R3K5HF9//TWePn2K8PBwtSqTwH9fw6jNY++/U61aNSxYsADJyclYs2YN7ty5g1atWqFFixbYs2cPCgsLRUckIioRLJREGmDRokU4cOAAgoKCYG1tLTpOERUrVkTnzp018iHnylCmTBmMGzcOt27dwv79+1G6dGlIpVLY2NhgxYoVePv2reiIREQqxUJJpOb27t2LH3/8ET///DO+/PJL0XH+llQqxfnz5/H48WPRUYTR09NDjx49EBERgejoaLi6usLDwwPm5uaYO3cunjx5IjoiEZFKcA8lkRq7c+cOmjVr9v4Vh+q8L+/ly5eoXr06fHx8MHXqVNFx1IZMJsOyZcuwbt065OTkYOjQoZg5cybs7e1FRyMiUhoWSiI19fr1azRr1gylSpXCpUuXUK5cOdGR/tWXX36Jt2/f4syZM6KjqJ3Xr19j3bp1WLZsGR4+fIguXbrAy8sLnTp1UusfFIiIPgZH3kRqSC6XY9iwYXj+/DnCw8M1okwC/x17nzt3jqPdv1ChQgV4eXkhMTER27dvx4sXL9ClSxc4Ojpiy5YtyMvLEx2RiOizsVASqaEffvgBBw8exI4dOzTq2Ya9e/eGvr6+zp32/hSlSpXCkCFDcPXqVZw8eRJ16tTBiBEjULduXfz66694+fKl6IhERJ+MI28iNRMaGgp3d3f88ssvmDNnjug4n6x79+7IyspCZGSk6Cga4/bt2/D19cXWrVtRqlQpjBkzBjNmzEDdunVFRyMi+igslERq5Pbt23B1dUX37t2xa9cujdxbt3HjRowdOxaPHj1CzZo1RcfRKE+fPsXKlSvx+++/49WrV3B3d4enpydcXV1FRyMi+kcslERq4tWrV2jWrBlKly6NCxcuaMy+yf+Vnp6O6tWrw9/fH5MnTxYdRyNlZWVh8+bN8PX1RXx8PFq1agVPT0/07NkT+vr6ouMRERXBPZREaqCwsBBDhgxBamoqwsLCNLZMAkDlypXRsWNHnX3IuTKULVsWEydOxN27dxEWFgaFQoG+ffuiYcOGWLVqFbKyskRHJCL6AAslkRpYsGABjh49ip07d6J+/fqi4xSbVCrFmTNn8OzZM9FRNJq+vj769OmDc+fO4eLFi3BwcMCUKVNgbm6OBQsW8L8vEakNFkoiwUJCQrB48WL88ssv6NKli+g4StGnTx/o6ekhNDRUdBSt4ebmhuDgYMTFxWHo0KHw9fWFhYUFxo0bhzt37oiOR0Q6jnsoiQS6efMm3Nzc0KNHD+zYsUMjD+H8na5duyI/Px+nTp0SHUUrvXz5EmvWrMHy5cvx5MkTfPXVV/D09ES7du206s8REWkGrlASCfLy5Uv06dMH9evXx4YNG7SuBPTv3x+RkZF4/vy56ChaqVKlSpgzZw6SkpIQGBgImUyGDh06wNnZGUFBQcjPzxcdkYh0CAslkQCFhYUYPHgwXr58ifDwcBgbG4uOpHR9+/aFRCLh2FvFDA0NMWLECMTExODo0aOoWrUqhg4divr168PHxwdv3rwRHZGIdABH3kQCfPfdd1i6dCmOHj2KTp06iY6jMp07d4ZcLsfJkydFR9EpN27cgK+vL4KCglCmTBmMGzcO06dPh5mZmehoRKSluEJJVMJ2796NJUuWYOnSpVpdJoH/nvY+ffo0x94lrEmTJggMDERSUhImTZqEDRs2oG7duhg6dCiuXbsmOh4RaSGuUBKVoBs3bqB58+bo06cPtm3bpnX7Jv/XixcvULNmTaxcuRITJkwQHUdnvX37Fhs3boS/vz8ePHiAdu3awcvLC927d4eeHtcViKj4WCiJSkh6ejqcnZ1RoUIFnD9/HmXLlhUdqUS8W4U9ceKE4CRUUFCAsLAw+Pj44PLly7C1tYWnpyeGDRuG0qVLi45HRBqMP5oSlYB3h3DevHmDsLAwnSmTwH/H3hEREXjx4oXoKDrPwMAAUqkUFy9exLlz52Bra4vx48fDwsICP/30E1JTU0VHJCINxUJJVALmzp2LkydPYteuXbC0tBQdp0T17dsXABAWFiY4Cb0jkUjQsmVLhIWF4d69e3B3d8cvv/wCc3NzTJo0CXFxcaIjEpGG4cibSMV27tyJwYMHw9fXFx4eHqLjCNGxY0fo6enh+PHjoqPQ30hNTcWqVasQEBCAFy9eoFevXvDy8kLLli21fq8vERUfVyiJVCgmJgajR4/G0KFDMWPGDNFxhHk39uZIVX1VrVoV8+fPR3JyMtatW4f79++jdevWcHNzw+7du1FQUCA6IhGpMRZKIhVJS0tDnz59YGtri7Vr1+r0Kk+/fv2gUCg49tYApUuXxpgxY3Dz5k0cPHgQxsbGGDhwIKytrbFs2TJkZGSIjkhEaogjbyIVKCgoQPfu3XH9+nVER0fDwsJCdCThOnToAAMDAxw7dkx0FPpE165dg4+PD3bt2gUTExNMmDABU6dORe3atUVHIyI1wRVKIhWYM2cOIiIiEBwczDL5f6RSKU6dOsWxtwb64osvsH37diQmJmLMmDH4/fffUbduXYwYMQI3btwQHY+I1AALJZGSBQUFwcfHBz4+PmjXrp3oOGrj3dg7PDxcdBT6TObm5vD29sbDhw+xZMkSREREwMHBAV26dMHRo0fBgReR7uLIm0iJrl+/jhYtWkAqlSIwMFCn903+lfbt28PQ0BBHjx4VHYWUID8/H3v27IG3tzeuXbsGe3t7zJw5E4MHD4aRkZHoeERUgrhCSaQkqamp6NOnDxo1aoTVq1ezTP4FqVSKkydPIi0tTXQUUoJSpUph8ODBiI6ORkREBCwsLDBq1CjUrVsXv/zyC9LT00VHJKISwkJJpAQFBQUYOHAgsrKyEBYWhjJlyoiOpJb69esHuVzOsbeWkUgkaNeuHfbv34/bt2+jR48e+PHHH2FmZoZp06YhMTFRdEQiUjGOvImUYObMmVixYgVOnDiBtm3bio6j1tq1a4fSpUvjyJEjoqOQCj179gy///47Vq5ciZcvX6Jfv37w9PSEm5ub6GhEpAJcoSQqpm3btsHPzw++vr4skx/h3dib41DtVr16dfz444+QyWRYuXIlYmJi0Lx5c7Rq1QphYWEoLCwUHZGIlIiFkqgYrl27hnHjxmHkyJGYMmWK6Dgawd3dHYWFhdi7d6/oKFQCypYti2+++QZ3795FeHg49PT00K9fP9ja2uL3339HVlaW6IhEpAQceRN9phcvXsDZ2RnVq1fHmTNnULp0adGRNEbbtm1hbGyMQ4cOiY5CAkRFRcHHxwd79uxBxYoVMXHiREyZMgU1atQQHY2IPhNXKIk+Q35+PgYMGICcnByEhoayTH4iqVSKEydO4OXLl6KjkADNmjXDrl27EB8fj6+//hr+/v6wsLDA2LFjcfv2bdHxiOgzsFASfYZZs2bh3Llz2LNnD+rUqSM6jsZxd3dHQUEBx946rm7duvD390dKSgoWLVqEw4cPw87ODl9++SVOnTrFB6UTaRAWSqJPtGXLFixbtgzLli1D69atRcfRSDVr1kSrVq0QHBwsOgqpgUqVKmH27Nl48OABNm/ejEePHqFjx47vX/mYn58vOiIR/QsWSqJPEB0djfHjx2P06NGYOHGi6DgarX///jh+/DhevXolOgqpCUNDQwwfPhzXr1/H8ePHUb16dQwbNgz16tXDb7/9htevX4uOSER/g4dyiD7S8+fP4eTkhFq1aiEyMpL7Jovp0aNHqFOnDgIDAzFixAjRcUhN3bx5E76+vti2bRtKly6NsWPHYvr06bCwsBAdjYj+hIWS6CPk5+ejU6dOuHfvHq5evYratWuLjqQVWrVqhYoVK+LAgQOio5Cae/LkCQICArBq1Sq8efMGUqkUnp6ecHZ2Fh2NiMCRN9FHmTlzJi5evIiQkBCWSSWSSqU4duwYx970r2rWrIn//Oc/kMlk8Pf3x+XLl+Hi4vL+lY9yuVx0RCKdxkJJ9C8CAwMREBCA5cuXo2XLlqLjaBV3d3fk5+dj3759oqOQhihXrhymTJmCuLg47NmzB7m5uejVqxcaNWqEtWvXIjs7W3REIp3EkTfRP7hy5Qpat26Nr7/+GmvXroVEIhEdSeu0bNkSlStXxv79+0VHIQ114cIF+Pj4ICwsDFWrVsXkyZMxadIkVKtWTXQ0Ip3BQkn0N549ewYnJyeYmZnh9OnTMDIyEh1JK/n7+2P27Nl4/vw5KlSoIDoOabD4+Hj4+/tj06ZNkMvlGDFiBDw8PNCgQQPR0Yi0HkfeRH8hLy8P/fv3R2FhIUJCQlgmVah///7Iy8vj2JuKzcrKCgEBAZDJZPj+++8RHh4OW1tb9O7dG2fOnOGD0olUiIWS6C94eHjg8uXLCAkJQa1atUTH0Wp16tRB8+bN+ZBzUpoqVapg3rx5SEpKwoYNGxAfH4+2bdvC1dUVu3btQkFBgeiIRFqHhZLof2zcuBG///47AgIC0KJFC9FxdIJUKsXRo0f54GpSqtKlS2P06NG4efMmDh06hPLly2PQoEGwsrKCv78/MjIyREck0hrcQ0n0J5cvX0abNm0watQorF69WnQcnZGSkgJzc3Ns3boVw4YNEx2HtNj169fh4+ODnTt3wtjYGBMmTMDUqVNRp04d0dGINBoLJdH/efr0KZycnGBpaYmIiAgYGhqKjqRTmjdvDlNTU+zdu1d0FNIBDx8+xPLly7FmzRpkZWVh0KBB8PT0hKOjo+hoRBqJI28i/P9DOACwZ88elkkB3o2937x5IzoK6YA6depg6dKlSElJwdKlS3HmzBk0bdoUnTp1wpEjR3iAh+gTsVASAZg+fTquXLmCkJAQ1KxZU3QcndS/f3/k5ubyeZRUosqXLw8PDw8kJCRgx44deP36Nbp37w57e3ts2rQJubm5oiMSaQQWStJ569atw+rVq7Fy5Uq4ubmJjqOzzM3N4erqytPeJISBgQEGDRqEqKgoREZGol69ehg9ejQsLS2xePFipKeni45IpNa4h5J02sWLF9G2bVuMHTsWv//+u+g4Os/Hxwfz5s3D8+fPUb58edFxSMfdvXsXfn5+2LJlC/T09DBq1Ch4eHigfv36oqMRqR0WStJZjx8/hrOzM+rXr4+TJ09y36QaSE5OhqWlJYKCgjB48GDRcYgAAC9evHj/KLG0tDT07dsXnp6efKwY0Z+wUJJOys3NRfv27SGTyRAdHY0aNWqIjkT/x9XVFbVr10ZoaKjoKEQfyM7OxpYtW+Dr64v79++jefPm8PT0RJ8+faCvry86HpFQ3ENJOmnq1Km4du0aQkNDWSbVjFQqxeHDh/H27VvRUYg+UKZMGUyYMAF37tzBvn37UKpUKfTv3x82NjYICAhAZmam6IhEwrBQks5Zs2YN1q1bh1WrVqFZs2ai49D/6N+/P3JycnDgwAHRUYj+kp6eHnr27InIyEhcuXIFzZo1w4wZM2BmZoZ58+bhyZMnoiMSlTiOvEmnnD9/Hu3bt8f48eMREBAgOg79jWbNmsHMzAwhISGioxB9lOTkZCxbtgzr1q1DXl4ehg4dipkzZ6Jx48aioxGVCBZK0hmPHj2Cs7MzbGxscOLECZQqVUp0JPobv/32GxYsWIAXL16gXLlyouMQfbRXr15h3bp1WLZsGR49eoRu3brB09MTHTt2hEQiER2PSGU48iadkJubC3d3dxgYGCA4OJhlUs29G3sfPHhQdBSiT1KxYkXMmjULiYmJ2Lp1K548eYLOnTujadOm2Lp1K/Ly8kRHJFIJFkrSegqFApMnT8b169cRFhYGU1NT0ZHoX9StWxfOzs58yDlpLENDQwwbNgx//PEHTpw4gVq1amH48OGoW7culi5dilevXomOSKRULJSk9VavXo0NGzZgzZo1cHZ2Fh2HPlL//v1x6NAhnpwljSaRSNCxY0ccOnQIN2/eRLdu3TB//nyYmZnBw8MDSUlJoiMSKQX3UJJWO3v2LDp06IBJkyZh2bJlouPQJ0hMTET9+vWxa9cuDBgwQHQcIqV5+vQpAgICsGrVKrx69QpSqRSenp5wcXERHY3os7FQktZ6+PAhnJyc0LBhQxw/fpz7JjWQk5MT6tWrx9E3aaXMzEwEBgbCz88PCQkJaN26Nby8vNCjRw/o6XGASJqFf2JJK+Xk5KBfv34wMjLC7t27WSY1lFQqxcGDBzn2Jq1kbGyMyZMn4969ewgJCUFhYSF69+6Nhg0bYvXq1cjOzhYdkeijsVCS1lEoFJg0aRJiY2N5CEfDSaVSZGdn49ChQ6KjEKmMvr4++vXrh/Pnz+PChQuwt7fH5MmTYW5ujh9++AHPnz8XHZHoX3HkTVpn5cqVmDJlCrZs2YKvv/5adBwqpi+++AJWVlbYvXu36ChEJSYhIQH+/v7YuHEjCgsLMXz4cMycORO2traioxH9Ja5QklY5c+YMZsyYgRkzZrBMaol3Y++srCzRUYhKTP369bFixQqkpKRgwYIF2L9/Pxo2bPj+lY9cCyJ1wxVK0hopKSlwcnJC48aNcezYMRgYGIiOREoQHx8Pa2trBAcHo3///qLjEAmRm5uLHTt2wNvbG7du3YKTkxM8PT3Rv39/7hEntcAVStIK2dnZ6NevH8qUKYNdu3axTGoRKysrODo68qQ36TQjIyOMHDkSsbGxOHLkCCpVqoQhQ4bAysoKvr6+ePPmjeiIpONYKEnjKRQKTJw4Ebdu3UJ4eDiqVasmOhIpmVQqxYEDBzj2Jp0nkUjQtWtXHD9+HNevX0fbtm0xe/ZsmJmZYdasWUhJSREdkXQUCyVpvBUrVmDz5s1Yv349mjZtKjoOqYBUKkVWVhYOHz4sOgqR2nBwcMCWLVuQlJSEiRMnYt26dahXr977Vz4SlSTuoSSNdvr0aXTq1AnTp0+Hj4+P6DikQo6OjrC1tcXOnTtFRyFSSxkZGdi4cSP8/PyQnJyMDh06wNPTE926deOD0knl+CeMNJZMJoNUKkW7du3w66+/io5DKvZu7M2HPRP9NRMTE0yfPh3x8fHYtWsXMjIy8NVXX8He3h4bNmxATk6O6IikxVgoSSNlZ2ejb9++KFeuHHbu3MlDODpAKpUiMzMTR44cER2FSK0ZGBhgwIABuHz5Ms6cOQNra2uMGzcOFhYW+Pnnn5GWliY6ImkhjrxJ4ygUCgwfPhwhISG4ePEiHBwcREeiEuLg4AA7OzsEBQWJjkKkUe7fvw8/Pz8EBgZCIpFg1KhR8PDwgJWVlehopCW4QkkaZ9myZdi2bRs2btzIMqljpFIp9u/fz7E30SeysbHBqlWrIJPJMGfOHAQHB8PGxub9Kx+5tkTFxUJJGuXUqVPw8vLCrFmzMGjQINFxqIRJpVK8ffsWR48eFR2FSCNVq1YNCxYsQHJyMtasWYM7d+6gVatWaN68Ofbs2YPCwkLREUlDceRNGiMpKQnOzs744osvcPjwYejr64uORAI0adIE9vb22L59u+goRBpPLpfj0KFD8PHxwenTp1G3bl14eHhg1KhRKFeunOh4pEG4QkkaISsrC3379kX58uWxc+dOlkkd9m7szROrRMWnp6eHHj16ICIiAtHR0XBzc4OHhwfMzMzw3Xff4fHjx6IjkoZgoSS1p1AoMG7cONy/fx/h4eGoXLmy6EgkkFQqRUZGBsfeRErm5OSEoKAgJCYmYvTo0Vi5ciUsLS3fv/KR6J+wUJLa8/PzQ1BQEDZt2oQmTZqIjkOC2draonHjxny3N5GKmJubw8fHBykpKVi8eDFOnjyJJk2avH/lI3fK0V9hoSS1duLECcyaNQuzZ8/GgAEDRMchNSGVSrFv3z6OvYlUqEKFCvDy8kJiYiK2b9+OFy9eoEuXLnB0dMTmzZuRl5cnOiKpER7KIbX14MEDODs7w8XFBQcPHuS+SXrvzp07aNSoEfbu3YtevXqJjkOkExQKBU6fPg1vb28cOnQItWrVwtSpUzFhwgRUqlRJdDwSjIWS1FJWVhZatGiBjIwMXLlyhfsmqYjGjRujadOm2Lp1q+goRDrn9u3b8PX1xdatW1GqVCmMGTMGM2bMQN26dUVHI0E48ia1o1AoMGbMGMTHx/MQDv2t/v37Y9++fcjNzRUdhUjnNGrUCOvXr0dycjI8PDywbds2WFlZvX/lI+keFkpSO97e3ti5cycCAwNhb28vOg6pKalUijdv3uDYsWOioxDprBo1auCnn35CSkoKAgIC8Mcff8DNzQ2tW7dGeHg4H5SuQ1goSa0cP34cc+bMwXfffYf+/fuLjkNqzM7ODg0bNuRpbyI1ULZsWUycOBF3795FWFgYAKBv375o2LAhVq1ahaysLMEJSdW4h5LURmJiIpydneHm5ob9+/fzEA79q4ULF8Lf3x/Pnz+HkZGR6DhE9CeXL1+Gj48PQkJCUKlSJUyaNAmTJ09G9erVRUcjFWChJLWQmZmJ5s2bIzs7G1FRUTwxSB/l5s2bsLe3x/79+9GjRw/RcYjoLyQmJmLZsmXYsGEDCgoK8PXXX2PmzJlo2LCh6GikRBx5k3AKhQKjR4/GgwcPEB4ezjJJH83Ozg62trYcexOpsXr16mHZsmVISUnBDz/8gIMHD6JRo0bvX/nIdS3twEJJwi1duhS7d+/G5s2bYWdnJzoOaRCJRAKpVIq9e/fytDeRmqtUqRLmzJmDpKQkBAYGQiaToUOHDnB2dkZQUBDy8/NFR6RiYKEkoY4ePYrvvvsO8+bNQ79+/UTHIQ0klUrx+vVrnDhxQnQUIvoIhoaGGDFiBGJiYnD06FFUrVoVQ4cORf369eHj44PXr1+LjkifgXsoSZj4+Hi4uLigZcuW2Lt3Lw/h0GdRKBRo2LAh3NzcEBgYKDoOEX2GGzduwNfXF0FBQShdujTGjRuH6dOnw9zcXHQ0+kgslCTE27dv4ebmhry8PERFRaFixYqiI5EGmz9/PgICAvDs2TMYGhqKjkNEn+nx48dYsWIFVq9ejYyMDAwYMACenp5wcnISHY3+BUfeVOIUCgVGjRqF5ORkhIeHs0xSsUmlUrx69YpjbyINV6tWLfzyyy9ISUmBr68vLl26BGdnZ7Rv3x4HDx6EXC4XHZH+BgsllbglS5Zgz5492Lp1Kxo1aiQ6DmkBe3t72NjYYM+ePaKjEJESlCtXDtOmTcP9+/exe/duZGdno0ePHrCzs8P69euRk5MjOiL9DxZKKlGHDh3CvHnzsGDBAvTp00d0HNIS7057h4eH86QokRYxMDCAVCrFxYsXce7cOdja2mL8+PGwsLDAokWLkJqaKjoi/R/uoaQSExcXBxcXF7Rp0wbh4eHQ0+PPM6Q8MTExcHR0xOHDh9GtWzfRcYhIReLi4uDn5/f+EN6IESPg4eEBGxsbscF0HAsllYiMjAy4ubmhsLAQly9fRoUKFURHIi2jUCjQoEEDtG7dGhs2bBAdh4hULDU1FatWrUJAQABevHiBXr16wdPTE61atYJEIhEdT+dwiYhUTi6XY8SIEUhJSUF4eDjLJKkEx95EuqVq1aqYP38+kpOTsW7dOty/fx9t2rSBm5sbdu/ejYKCAtERdQoLJanc4sWLERYWhm3btsHW1lZ0HNJiUqkU6enpOHXqlOgoRFRCSpcujTFjxuDmzZs4ePAgjI2NMXDgQFhbW2PZsmXIyMhQyX0zcwtw6/Fr/CF7iVuPXyMzV7cLLEfepFIHDx5Ez549sXDhQixcuFB0HNJyCoUCNjY2aNu2LdavXy86DhEJcu3aNfj6+mLnzp0wMTHBhAkTMHXqVNSuXbtY1417loHtl2WIuPccsvQs/LlASQCYVy6L9g1MMdTVHNbVTYp1L03DQkkqc//+fbi4uKB9+/YIDQ3lIRwqEXPnzsWaNWvw9OlTlCpVSnQcIhIoJSUFy5cvx5o1a5CTk4PBgwfD09MTTZo0+bTrpGdhblgszsanQl9PgkL531endx9vbVUVi/vaw6xy2eJ+GhqBhZJU4s2bN3Bzc4NCocDly5dRvnx50ZFIR/zxxx/44osvcPToUXTp0kV0HCJSA2/evMH69evh7++PlJQUdO7cGZ6enujSpcu/HuDZeUWGhftuoUCu+Mci+b/09SQw0JPgx152GOSi/a+Q5JIRKZ1cLsfw4cPx6NEjhIeHs0xSiXJ0dET9+vURHBwsOgoRqYny5ctj5syZSEhIQFBQENLS0tCtWzc0adIEgYGByM3N/cvfFxARhzmhscgtkH9SmQSAQrkCuQVyzAmNRUBEnDI+DbXGQklK9/PPP2Pfvn3Yvn07GjRoIDoO6Zh3p73DwsJ42puIPlCqVCkMHjwY0dHRiIiIgKWlJUaNGoW6devil19+QXp6+vtfu/OKDN7H7ivlvt7H7mPXFZlSrqWuOPImpdq/fz969eqFRYsWYf78+aLjkI66du0anJyccOzYMXTu3Fl0HCJSY3fu3IGfnx+2bNkCfX19jBkzBgNHT8LokAfILVDeu8ONDPRwwqOt1u6pZKEkpbl79y5cXV3RsWNH7Nmzh4dwSBiFQgErKyt07NgRa9euFR2HiDTA8+fPsXLlSqxcuRIGnWagtKUDIFHe9zF9PQla1KuCrWNclXZNdcJCSUrx+vVruLq6Ql9fH5cuXYKJiW49LoHUz+zZs7Fx40Y8efIEBgYGouMQkYa4kfwCvVZHqez6JzzawMpU+75HcgmJik0ul+Prr7/G06dPER4ezjJJakEqlSI1NRWnT58WHYWINEhozHPo66nm1Y36ehJsu6SdeylZKKnYFi1ahAMHDiAoKAjW1tai4xABAJycnGBpacnT3kT0SSLuPf/kE90fq1CuQMT95yq5tmgslFQse/fuxY8//oiff/4ZX375peg4RO+9O+0dGhrKd/oS0Ud5m1sAWXqWSu8hS8vSytc0slDSZ7tz5w6GDRsGd3d3fPfdd6LjEBXxbuwdGRkpOgoRaYDktEyo+mCJAkBSWqaK71LyWCjps7x+/Rp9+vSBhYUFAgMD//VNA0QiODs7c+xNRB8tT4mPCVKH+5QkHn2kTyaXyzF06FA8f/4cV65cQbly5URHIvpLEokE/fv3x+bNmxEQEMDT3kT0AblcjoSEBFy/fh3Xr1/HpXsPAesBKr+voYH2refxqyt9soULF+LQoUM4dOgQrKysRMch+kdSqRTe3t44c+YMOnToIDoOEQmSlZWFmzdvvi+PMTExiImJQWbmf8fPNWvWRJMvXPDfobTqpm4SAJZVjFV2fVFYKOmThIaG4ueff8Yvv/yCbt26iY5D9K9cXFxgYWGBPXv2sFAS6Yhnz569L47vyuO9e/cgl8uhp6cHW1tbODo6ok+fPnB0dISDgwNMTU0BAG1/i0CyCg/mmFcpC2Mj7atffLA5fbTbt2/D1dUV3bt3x65du7hvkjSGl5cXtm3bhkePHkFfX190HCJSksLCQsTFxRUpj0+fPgUAlCtXDg4ODnB0dHz/j52dHcqUKfO31/xh3y1svZyskkcH6etJ8LWrBX7oZaf0a4vGQkkf5dWrV3BxcUGZMmVw4cIF7pskjXL58mW4ubkhIiIC7dq1Ex2HiD7D27dvERsb+0F5jI2NRXZ2NgDAzMysSHmsW7fuJ78GOO5ZBjr7n1HFpwBAe9+Uo31rrqR0hYWFGDJkCNLS0hAdHc0ySRqnWbNmMDc3R3BwMAslkZpTKBR4/PgxYmJiPiiP8fHxUCgUMDAwQKNGjeDo6IiBAwe+H1lXqVJFKfe3rm6C1lZVcSExTamrlO/e5a2NZRLgCiV9hHnz5mHJkiU4fPgwunTpIjoO0Wfx9PTE9u3bOfYmUiP5+fm4d+9ekfKYmpoKAKhQocL71cZ3q4+NGjWCkZGRSnOlpGehk18kcpX4eB8jAz2c8GgLs8pllXZNdcJCSf8oJCQE/fv3x6+//opvv/1WdByiz3bp0iU0b94cp0+fRtu2bUXHIdI5r1+/xo0bN97vc7x+/Tpu3ryJ3NxcAIClpWWR8mhhYSFsv/7OKzLMCY1V2vV+7WePgS7mSrueumGhpL918+ZNuLm5oUePHtixYwcP4ZBGUygUsLCwQK9evRAQECA6DpHWUigUSElJ+eCQzPXr15GYmAgAMDQ0hJ2d3Qd7HZs0aYKKFSuKDf4XRnnvQkRaORT3UUKzujTA5Pba/Zg9Fkr6Sy9fvoSLiwuMjY1x4cIFGBtr3zOzSPfMnDkTO3bswMOHDzn2JlKCvLw83Llzp0h5fPnyJQCgcuXKHxRHR0dH2NraolSpUoKT/7uoqCi0bt0aHcfPR0JFJxTIFZ+0p1JfTwIDPQkW9bLT6pXJd1goqYjCwkJ89dVXuHLlCqKjo1G3bl3RkYiU4uLFi2jRogUiIyPRpk0b0XGINMrLly8/2OsYExODW7duIT8/HwBgZWVV5JR17dq1NXK69fTpUzg7O8Pc3BwRERF4nlmIuWGxOBufCn09yT8Wy3cfb21VFYv72mvtnsn/xUJJRXz33XdYunQpjh49ik6dOomOQ6Q0crkcFhYW6NOnD1asWCE6DpFaUigUSEpK+uCQzPXr1yGTyQAApUuXhr29/QflsUmTJjAx0Y7Ty3l5eejQoQMSExNx9epV1KxZ8/3H4p5lYPtlGSLuP4csLQt/LlAS/Peh5e1tTDHMzVxrT3P/HRZK+sDu3bsxcOBAeHt7w9PTU3QcIqXz8PDArl278PDhw09+Ph2RtsnJycHt27eLPBj8zZs3AIBq1aqhadOmH5RHGxsbGBho71MHJ02ahPXr1yMyMhLNmzf/21+XmVuApLRM5BXIYWigB8sqxlr5BpyPxUJJ7924cQPNmzdHnz59sG3bNo0cUxD9mwsXLqBly5Y4c+YMWrduLToOUYlJTU0t8nieO3fuoLCwEBKJBDY2NkVOWdeoUUOnvhds2LABY8eOxdq1azFu3DjRcTQKCyUBANLT0+Hs7IwKFSrg/PnzKFtWN/Z8kO6Ry+UwNzdHv379sHz5ctFxiJROLpcjISHhg0My169fx6NHjwAAZcuWRZMmTT4oj/b29jp/+PLSpUto27YtRo8ejVWrVomOo3FYKAmFhYX48ssvcfXqVURHR8PS0lJ0JCKVmjFjBoKDg5GSksKxN2m0rKws3Lx584PyGBMTg8zMTABAzZo1i5yyrl+/Pp9y8D+ePHkCJycn1KtXD6dOnYKhoaHoSBqHhZIwe/Zs+Pj44OjRo+jYsaPoOEQqd/78ebRq1Qpnz55Fq1atRMch+ijPnj0r8niee/fuQS6XQ19fHw0aNPigODo4OMDU1FR0bLWXl5eH9u3bIykpCVevXkWNGjVER9JIurt7lAAAO3fuxNKlS+Hr68sySTqjefPmqFWrFoKDg1koSe0UFhYiLi6uyEGZp0+fAgBMTEzg4OCADh06YObMmXB0dISdnR3KlCkjOLlmmjZtGqKjo3HmzBmWyWLgCqUOi4mJQfPmzdGvXz9s3bpVpzZeE02bNg0hISEce5NQb9++RWxs7AflMTY2FtnZ2QAAMzOzDw7JODo6om7duvwzqyRr167FhAkTsGHDBowePVp0HI3GQqmj0tLS4OzsjMqVK+PcuXP8yZZ0ztmzZ9GmTRucO3cOLVu2FB2HtJxCocDjx4+LnLKOj4+HQqGAgYEBGjVq9EF5dHBwQJUqVURH11oXLlxAu3btMG7cOKxcuVJ0HI3HQqmDCgoK0K1bN8TExODq1aswN9f+V0IR/S+5XI46depgwIAB8Pf3Fx2HtEhBQQHu3r1bpDympqYCACpUqFDk8TyNGjWCkZGR4OS64/Hjx3BycoKVlRVOnjzJQzhKwEKpg7y8vODv748TJ06gXbt2ouMQCTN16lSEhYVBJpNxhEif5c2bN4iJifmgPN68eRO5ubkAAEtLyyIHZSwsLLjFSKDc3Fy0a9cOKSkpuHr1KqpXry46klZgodQxQUFBGDp0KPz9/TF9+nTRcYiEOnPmDNq2bYvz58+jRYsWouOQGlMoFEhJSSnybMfExEQAgKGhIezs7D4oj02aNEHFihXFBqcPKBQKjB8/Hlu3bsXZs2fh4uIiOpLWYKHUIX/88QdatmwJqVSKwMBA/oRMOq+wsBB16tTBoEGD4OfnJzoOqYm8vDzcuXOnyCN6Xr58CQCoXLlykWc72traolSpUoKT079ZvXo1Jk6ciE2bNmHkyJGi42gVFkodkZqaCmdnZ1StWhVnz57lIRyi/zNlyhTs3bsXycnJHHvroJcvX36w4hgTE4Nbt24hPz8fAGBlZfXBCWtHR0fUrl2bP5BroPPnz6N9+/aYMGECVqxYITqO1mGh1AEFBQXo2rUrYmNjcfXqVZiZmYmORKQ2IiMj0a5dO1y8eBFubm6i45CKKBQKJCUlFXm2Y3JyMgCgdOnSsLe3/+CgTJMmTWBiYiI4OSnDo0eP4OTkhAYNGuDEiRNcTVYBPthcB3z77bc4c+YMTpw4wTJJ9D9atWqF6tWrIzg4mIVSS+Tk5OD27dtFyuObN28AANWqVUPTpk0xcODA9+XRxsYGBgb8lqiNcnJy0K9fPxgaGiI4OJhlUkW4Qqnltm3bhq+//horVqzAlClTRMchUkuTJ0/GgQMHkJSUxFGmhklNTS3yeJ67d++ioKAAEokENjY2RR7RU6NGDf5/1hEKhQJjx47F9u3bce7cOTg7O4uOpLVYKLXYtWvX0LJlSwwaNAgbN27kF1Civ3H69Gm0b98ely5dgqurq+g49BfkcjkSEhKKlMdHjx4BAMqWLYsmTZp8sNexcePGMDY2FpycRPr9998xefJkbN68GcOHDxcdR6uxUGqpFy9ewNnZGdWrV8eZM2dQunRp0ZGI1FZhYSFq166NYcOGwdvbW3QcnZeVlYWbN29+UB5v3LiBt2/fAgBq1apV5KBM/fr1oa+vLzg5qZOzZ8+iQ4cOmDx5Ml9eUAJYKLVQfn4+unTpgtu3b+Pq1auoU6eO6EhEam/SpEk4dOgQHjx4wNX8EvTs2bMij+e5d+8e5HI59PX1YWtr+0F5dHBwgKmpqejYpOZSUlLg7OyMRo0a4dixY9w3WQJYKLXQjBkzsHLlSpw6dQqtW7cWHYdII0RERKBDhw64fPkymjVrJjqO1iksLERcXFyR8vj06VMAgImJCRwcHD4oj3Z2dnzEGX2ynJwctGnTBs+ePUN0dDSqVasmOpJO4JE2LbNlyxYsW7YMK1euZJkk+gRt2rSBqakpgoODWSiL6e3bt4iNjf2gPN64cQPZ2dkAADMzMzg6OmLcuHHvC2TdunX5HFAqNoVCgYkTJyI2Nhbnz59nmSxBXKHUItHR0WjVqhWGDh2K9evXc2xH9IkmTpyIw4cPc+z9kRQKBZ48efLBIZnr168jPj4eCoUCBgYGaNSo0QfjagcHB1SpUkV0dNJSAQEBmDp1KrZu3Yphw4aJjqNTWCi1xPPnz+Hk5IRatWohMjKSh3CIPsOpU6fQsWNHREVF8R2//6OgoAD37t0rUh5TU1MBABUqVPjgkIyDgwMaNWoEIyMjwclJV0RGRqJjx46YNm0afH19RcfROSyUWiA/Px+dOnXCvXv3cPXqVdSuXVt0JCKNVFBQgFq1amHkyJFYunQpMnMLkJSWibwCOQwN9GBZxRjGRtq/U+jNmzeIiYn54JT1zZs3kZubCwCwtLQsUh4tLCy4qkvCyGQyODs7w97eHkePHuVD6gVgodQCU6dOxZo1axAREYGWLVuKjkOk0YZN8sK5ZxKYu32JlPQs/PkLpASAeeWyaN/AFENdzWFdXbNfy6dQKJCSkvLBIZnr168jMTERAGBoaAg7O7sPymOTJk1QsWJFscGJ/iQ7OxutW7dGamoqoqOjUbVqVdGRdBILpYYLDAzEqFGjsGrVKnzzzTei4xBprJT0LMwNi8XZ+FQo5IWQ6P39Mw319SQolCvQ2qoqFve1h1nlsiWY9PPk5eXhzp07Rcrjy5cvAQCVK1dG06ZNPzhlbWtry8etkFpTKBQYOXIkgoODcf78eTRt2lR0JJ3FQqnBoqKi0KZNGwwfPhxr1qzhuInoM+28IsPCfbdQIFegUP7xXxL19SQw0JPgx152GORirsKEn+bly5fvS+O7f9+6dQv5+fkAACsrqw9eRejo6IjatWvzawhpnGXLlmHGjBkICgrC4MGDRcfRaSyUGurp06dwdnaGubk5IiIiuPGd6DMFRMTB+9j9Yl/Hq4sNprS3VkKij6dQKJCUlFTk2Y7JyckAgNKlS8Pe3v6D8tikSROYmGj2qJ4I+O+zYzt37owZM2bwDVdqgIVSA+Xl5aFjx46Ij4/H1atXUatWLdGRiDTSzisyzAmNVdr1fu1nj4EqWqnMycnB7du3PzhhHRMTgzdv3gAAqlWrhqZNm35wUMbGxoaHE0grJScnw9nZGQ4ODjhy5Aj/nKsBFkoNNHnyZKxbtw6nT59GixYtRMch0kgp6Vno5BeJ3AK50q5pZKCHEx5ti72nMjU19YN9jtevX8fdu3dRUFAAiUQCGxubIqesa9SowZE16YSsrCy0atUKL1++RHR0NJ9rqiZY6TXMhg0b8Pvvv2Pt2rUsk0TFMDcsFgWfsF/yYxTIFZgbFoutY1w/6tfL5XIkJiYWebbjo0ePAABly5ZFkyZN0KpVK0yZMgWOjo5o3LgxjI2NlZqbSFMoFAqMHz8ed+/excWLF1km1QgLpQa5dOkSJk2ahAkTJmDcuHGi4xBprLhnGTgbn6r06xbKFTgbn4r45xmwMv1wn2JWVhZu3rz5wcrjjRs38PbtWwBArVq14ODggOHDh79feaxfvz709f/+tDmRrvH398f27duxY8cOODg4iI5Df8KRt4Z4+vQpnJycYGlpiYiICBgaGoqORKSxfth3C1svJ3/Sie6Ppa8ngXuTauhYMf2D8njv3j3I5XLo6+vD1tb2g4MyDg4OMDU1VXoWIm1y8uRJdO3aFZ6envj1119Fx6H/wUKpAfLy8tChQwc8ePAA0dHRqFmzpuhIRBqt7W8RSE7PUtn189Mf4/Ha8TAxMfmgNDo6OsLOzg5lypRR2b2JtFFSUhKcnZ3h5OSEQ4cOceVeDXHkrQGmT5+OK1euIDIykmWSqJje5hZApsIyCQClKtdE7J37aGRTH3p6eiq9F5G2y8rKQt++fVGhQgXs2LGDZVJNsVCquXXr1mH16tVYv3493NzcRMch0njJaZlQ/VhGAkl5U5ZJomJSKBQYO3Ys7t+/j0uXLqFy5cqiI9HfYKFUYxcvXsTkyZMxceJEjBkzRnQcIq2Qp8THBKnDfYi0ma+vL3bs2IHdu3fD3t5edBz6ByyUaurx48dwd3eHq6sr/P39Rcch0nhv3rzBjRs3cCzqNoDaKr+foQFXJ4mK4/jx4/j2228xZ84cSKVS0XHoX/BQjhrKzc1F+/btIZPJEB0djRo1aoiORKQxFAoFHj58WOSNMgkJCQAAI+PyqDFlO6DCh4BLANz8oSuMjfgzO9HnePDgAZydndGsWTMcOHCA+yY1AL/aqaGpU6fi2rVrOHPmDMsk0T/Iz8/HnTt3ipTH9PR0AEDlypXh6OiI3r17v3+2o62tLTr5n1PpKW/zKmVZJok+U2ZmJvr06YNKlSohKCiIZVJD8CuemlmzZg3WrVuHjRs3olmzZqLjEKmNV69eFXkd4e3bt5GXlwcAqF+/PhwdHeHh4fG+PNauXfsvX0fYvoGpSp9D2d6Gz5Qk+hwKhQJjxoxBQkICLl26hEqVKomORB+JhVKNnD9/HlOnTsXkyZMxatQo0XGIhFAoFEhOTi7yOsLk5GQAQOnSpdG4cWM4Oztj7NixcHR0hL29PcqXL//R9xjqao7Ai0kqyV8oV2CYm7lKrk2k7X777Tfs2rULe/bsQePGjUXHoU/APZRq4tGjR3B2doaNjQ1OnDiBUqVKiY5EpHK5ubm4fft2kZH169evAQDVqlV7v9r47h8bGxsYGBT/Z+GvN1zGhcQ0pa5S6utJ0KJelY9+lzcR/X/Hjh1D9+7dMWfOHPznP/8RHYc+EQulGsjNzUXbtm3x6NEjXL16la9gI62UlpZWZGR9584dFBQUQCKRwNraukh5rFGjxl+OrJUhJT0LnfwikavEx/sYGejhhEdbmFUuq7RrEumChIQEuLi4wM3NDfv37+e+SQ3EQimYQqHAuHHjsG3bNpw7dw7Ozs6iIxEVi1wux4MHD4qMrB8+fAgAKFu2LOzt7T8ojvb29jA2Ni7xrDuvyDAnNFZp1/u1nz0GunDcTfQp3r59ixYtWiAnJwdRUVGoWLGi6Ej0GbiHUrDVq1djw4YNCAwMZJkkjZOdnY1bt24VGVm/ffsWAFCzZk04ODhg2LBh78ujlZWV2qw+DHIxR+rbXHgfu1/sa83q0oBlkugTKRQKjBo1Cg8ePMDly5dZJjUYVygFOnv2LDp06IBJkyZh2bJlouMQ/aPnz58XGVnfvXsXcrkcenp6sLW1hYODw/vi6ODggOrVq4uO/VF2XpFh4b5bKJArPmlPpb6eBAZ6EizqZccySfQZlixZgu+++w6hoaHo27ev6DhUDCyUgjx8+BBOTk5o2LAhjh8/zkM4pDYKCwuRkJBQZGT95MkTAEC5cuXQpEmTD0bWjRs3RpkyZQQnL56U9CzMDYvF2fhU6OtJ/rFYSqCAAhK0tqqKxX3tuWeS6DMcOXIEX375JebNm4effvpJdBwqJhZKAXJyctCmTRs8ffoU0dHRPIRDwmRlZSE2NvaD4njjxg1kZf33od916tT5YNXR0dER9erVg56e9r5WMO5ZBrZfliHi/nPI0rLw5y+QEgBGBRkokMXgyIq5sDI1ERWTSKPFx8fDxcUFrVq1wt69e7X6a4quYKEsYQqFAqNHj8bOnTtx7tw5ODk5iY5EOuLp06dFVh3v378PhUIBfX19NGrUqMjIumrVqqJjC5WZW4CktEzkFchhaKAHyyrGiDh+BD179sTNmzdhZ2cnOiKRxsnIyEDz5s2Rn5+PqKgoVKhQQXQkUgIeyilhK1euRGBgILZu3coySSpRUFCAuLi4IuXx+fPnAIDy5cvD0dERXbp0wbfffgtHR0c0atQIpUuXFpxc/RgbGcCu1off7Dp37ozy5csjODiYhZLoE707hCOTyXD58mWWSS3CFcoSdObMGXTs2BFTpkyBn5+f6DikBTIyMoqMrGNjY5GTkwMAsLCweL/a+G7l0dLSUmXPdtQVw4cPx9WrV3Hr1i3RUYg0yuLFizFv3jyEh4ejd+/eouOQErFQlpCUlBQ4OTmhcePGOHbsmFLe9EG6Q6FQ4PHjx0VWHRMSEqBQKFCqVCnY2dkVGVnzPbiqsX//fvTq1Qu3bt1Co0aNRMch0ggHDx5Ez549sWDBAvzwww+i45CSsVCWgOzsbLRu3Rqpqam4cuUKqlWrJjoSqbH8/Hzcu3evSHlMS0sDAFSqVKnIqmPDhg1haGgoOLnuyM3NhampKWbOnImFCxeKjkOk9uLi4uDi4oK2bdsiLCyMh3C0EAuliikUCowcORLBwcE4f/48mjZtKjoSqZHXr1/jxo0bHxTHmzdvIi8vDwBQr169IuXRzMyMI2s18PXXX+OPP/7AzZs3RUchUmsZGRlwc3NDYWEhoqKiUL58edGRSAU4d1WxFStWYMuWLdi+fTvLpA5TKBSQyWRFHgz+4MEDAICRkREaN26Mpk2bYuTIkXB0dESTJk24YV2NSaVSbNu2DXfu3EHDhg1FxyFSS3K5HMOHD8fDhw9ZJrUcC6UKnT59GjNnzsTMmTMxZMgQ0XGohOTl5eHOnTtFRtavXr0CAFSpUgVNmzaFu7v7+5XHBg0a8OH2GqZLly4wMTFBcHAwFixYIDoOkVr6z3/+g/DwcOzbtw8NGjQQHYdUiCNvFZHJZHBycoKDgwOOHDnCQzha6uXLl0VWHW/fvo38/HwAgLW1dZGRda1atTiy1hLDhg1DTEwMYmNjRUchUjsHDhxAr1698MMPP/CHLh3AQqkC2dnZaNWqFdLT0xEdHY0qVaqIjkTFpFAo8ODBgyLlUSaTAQBKly79/nWE78qjvb09TEz4JhVttnfvXvTp0wd37tyBra2t6DhEauPevXto1qwZOnTogJCQEB7C0QEslEqmUCgwfPhwhISE4OLFi3BwcBAdiT5RTk4Obt269UF5jImJwZs3bwAApqamaNq06Qfl0dramqvQOignJwempqaYNWsW5s+fLzoOkVp48+YNXF1dIZFIcOnSJe6b1BEslErm7+8PDw8P7NixA4MGDRIdh/5FampqkVXHO3fuoLCwEBKJBA0aNPjguY6Ojo6oUaOG6NikRoYOHYrY2FjcuHFDdBQi4eRyOfr27YvIyEhERUXBxsZGdCQqISyUSnTq1Cl06dIFM2fOxNKlS0XHoT+Ry+VISEgoUh4fPXoEAChbtuz7wvju340bN4axsbHg5KTuwsPD0bdvX9y9e5eHDkjn/fjjj/jxxx+xf/9+fPXVV6LjUAlioVSSpKQkODs744svvsDhw4ehr68vOpLOysrKws2bNz8ojzdu3MDbt28BALVq1Sqy6li/fn3+P6PPkp2dDVNTU8yePRvff/+96DhEwuzbtw+9e/fGTz/9xL8LOoiFUgmysrLQsmVLvH79GtHR0ahcubLoSDrj2bNn7/c4viuP9+7dg1wuh76+PmxtbT8ojw4ODjA1NRUdm7TMkCFD3u+7JdJFd+/eRbNmzdC5c2cEBwfzEI4OYqEsJoVCgWHDhiE8PBwXL15EkyZNREfSSoWFhYiLiytSHp8+fQoAMDExKfIeazs7O5QpU0ZwctIFYWFh6NevH+7du8c9Y6RzXr9+jWbNmsHAwACXLl3i0y10FI+lFpOfnx+CgoKwa9culkklefv2LWJjY4uMrLOzswEAZmZmcHR0xLhx496XyLp16/InYhKmW7duKFeuHIKDgzFv3jzRcYhKjFwux7Bhw/Ds2TNcuXKFZVKHcYWyGE6cOIGuXbti1qxZWLJkieg4GkehUODJkydFVh3j4uKgUChgYGCARo0aFRlZ87mepI4GDx78/g1JRLpi4cKF+Omnn3Dw4EF0795ddBwSiIXyMz148ADOzs5wcXHBwYMHeaDjXxQUFODevXtFyuOLFy8AABUqVChyUKZRo0YwMjISnJzo44SGhsLd3R3379+HtbW16DhEKvfuCQeLFy/Gd999JzoOCcZC+RkyMzPRsmVLvH37FlFRUTyE8z/evHmDGzdufFAeY2NjkZubCwCwtLQsUh4tLCz4OkLSaNnZ2ahWrRrmzp2LuXPnio5DpFK3b9+Gq6srunXrht27d/PrN7FQfiqFQoHBgwfjwIEDuHTpEho3biw6kjAKhQIPHz784G0y169fR0JCAgDA0NAQdnZ2RUbWFStWFBucSEUGDRqEe/fu4Y8//hAdhUhlXr16hWbNmsHIyAgXL15EuXLlREciNcBDOZ/I29sbu3btQnBwsE6Vyfz8/Pf7w/48sk5PTwcAVK5cGY6Ojujdu/f78mhrawtDQ0PByYlKjlQqRf/+/REfHw8rKyvRcYiU7t0hnBcvXiA6Opplkt7jCuUnOH78OLp164bZs2dj8eLFouOozKtXr96Xxnf/vnXrFvLy8gAA9evXf7/q+K481qlThyMP0nlZWVmoVq0avv/+e+4pI600f/58LF68GIcOHULXrl1FxyE1ovOFMjO3AElpmcgrkMPQQA+WVYxhbFR04TYxMRHOzs5wc3PD/v37teIQjkKhQHJycpGRdVJSEgDAyMgI9vb2HxTHJk2aoHz58mKDE6mxgQMHIi4uDteuXRMdhUip3h08W7JkCWbPni06DqkZnSyUcc8ysP2yDBH3nkOWnoU//weQADCvXBbtG5hiqKs5rKubIDMzE82bN0d2djauXLmikXsAc3Nzcfv27SLl8fXr1wCAatWqFVl1bNCgAQwMuCuC6FPs2bMHUqkU8fHxqF+/vug4REpx69YtuLq64quvvsLOnTs5kaIidKpQpqRnYW5YLM7Gp0JfT4JC+d9/6u8+3sqqKjJOrcWpfbtx6dIl2NnZlWDiz5OWllZkZH379m0UFBRAIpHA2tq6SHmsWbMmv0AQKcG7sff8+fMxZ84c0XGIiu3ly5do1qwZypQpg4sXL8LY2Fh0JFJDOlMod16RYeG+WyiQK/6xSP4vCRQozM/DIGs9LJ3QR3UBP4NcLseDBw+KrDqmpKQAAMqUKYMmTZp8UB7t7e35xYBIxQYMGICEhARcvXpVdBSiYiksLESPHj1w+fJlREdHo169eqIjkZrSiXlmQEQcvI/d/6zfq4AEeqUMsTtJAvOIOExpL+aBxdnZ2bh169YH5TEmJgYZGRkAgBo1asDR0RFDhw59v+pobW2tFXs9iTSNVCrFgAEDkJiYyG/ApNHmz5+PY8eO4ciRI/yzTP9I61cod16RYU5orNKu92s/ewx0MVfa9f7Kixcv3hfHd+Xx7t27KCwshJ6eHho0aFBkZF29enWVZiKij5eZmYlq1aph4cKFPLxAGuvdfuClS5di1qxZouOQmtPqQpmSnoVOfpHILZAr7ZpGBno44dEWZpXLFvtacrkc8fHxRcrj48ePAQDGxsbv3yTz7h87OzuULVv8exORakmlUjx48ADR0dGioxB9stjYWDRv3hw9e/ZEUFAQ99jTv9LqQvn1hsu4kJj2SXsm/42+ngQt6lXB1jGun/T7srKyEBsb+0FxvHHjBjIzMwEAtWvX/qA4Ojo6ol69etDT01NadiIqObt378bAgQORkJDAUSFplPT0dLi4uMDExAQXLlzgIgZ9FK3dQxn3LANn41OVft1CuQJn41MR/zwDVqYmf/lrnj59WmTV8f79+5DL5dDX10fDhg3h6OgId3f39yPrqlWrKj0rEYnz1VdfoUyZMtizZw++/fZb0XGIPkphYSEGDx6MV69e4cSJEyyT9NG0doXyh323sPVyslJXJ9/R15Pga1cLzP/KFvfv3y9SHp89ewYAKF++fJGRdaNGjVC6dGmlZyIi9dO/f38kJyfjypUroqMQfZQ5c+bgt99+w9GjR9GpUyfRcUiDaG2hbPtbBJLTs1R2ff2sNDxZ9w2ys7MBAObm5kVG1paWltx3QqTDdu3ahUGDBiExMRF169YVHYfoH73bpuHj44OZM2eKjkMaRisL5dvcAtj/cBQq/cQUCkwyTUCzL/77OsLKlSur8m5EpIHevn2LatWqYdGiRTwlS2otJiYGLVq0QN++fbF161YuhtAn08pCeevxa3y14pzK73NwaivY1aqg8vsQkeZyd3dHSkoKoqKiREch+ktpaWlwcXFBhQoVcP78ee6bpM+ilUeI85T4mCB1uA8RaS6pVIorV64gKSlJdBSiIgoKCjB48GC8efMGYWFhLJP02bSyUBoalMynVVL3ISLN1aNHD5QuXRp79uwRHYWoiLlz5+LUqVPYvXs3LC0tRcchDaaVjciyijFUvftD8n/3ISL6J+XKlUP37t0RHBwsOgrRB3bs2IHffvsN3t7e6NChg+g4pOG0slAaGxnAXAlvsvkn5lXKwthIax/jSURKJJVKERUVheTkZNFRiAAA169fx5gxYzBs2DBMnz5ddBzSAlpZKAGgfQNT6OupZp1SXwK0tzFVybWJSPv06NEDRkZGHHuTWkhLS0Pfvn3RsGFDrF27lie6SSm0tlAOdTVXyUPNAaBQASSd2IqEhASVXJ+ItIuJiQnH3qQWCgoKMHDgQGRmZiIsLAxlypQRHYm0hNYWSuvqJmhtVVXpq5T6EqCW5DUO7tgAa2truLu74+LFi0q9BxFpH6lUisuXL0Mmk4mOQjps9uzZOH36NHbv3g1zc3PRcUiLaG2hBIDFfe1hoORCaaCvh11efSCTybBq1SrcvHkTLVq0QIsWLRAaGorCwkKl3o+ItEPPnj059iahgoKC4OvrC19fX7Rr1050HNIyWl0ozSqXxY+97JR6zUW97GBWuSzKlCmDCRMm4M6dO9i7dy9KlSoFd3d3NGjQACtXrkRmZqZS70tEms3ExATdunXj2JuE+OOPPzBmzBgMHz4cU6dOFR2HtJBWvinnfwVExMH72P1iX2dWlwaY3N7qbz9+5coV+Pj4IDg4GBUrVsTEiRMxZcoU1KhRo9j3JiLNt337dgwbNgwymQxmZmai45COePHiBZydnWFqaoozZ85w3ySphFavUL4zpb01lvSzh5GB3ifvqdTXk8DIQA+/9rP/xzIJAC4uLti5cycSEhIwfPhw+Pv7w8LCAmPGjMGtW7eK8ykQkRbg2JtK2rtDODk5OQgNDWWZJJXRiRXKd1LSszA3LBZn41Ohryf5x1Pg7z7e2qoqFve1h9lnPNfy1atXWLt2LZYtW4bHjx+je/fu8PT0RIcOHfiYBiId1bt3b7x48QIXLlwQHYV0gIeHBwICAnDy5Em0adNGdBzSYjpVKN+Je5aB7ZdliLj/HLK0LPz5P4AE/31oeXsbUwxzM4eVqUmx75eXl4ddu3bB29sbN27cgKOjI7y8vDBgwACUKlWq2NcnIs2xbds2fP311xx7k8pt3boVw4cPR0BAACZPniw6Dmk5nSyUf5aZW4CktEzkFchhaKAHyyrGKnsDjkKhwMmTJ+Ht7Y2jR4+iTp06mDZtGsaPH48KFSqo5J5EpF5ev34NU1NT/Prrr5gxY4boOKSlrl69ilatWmHw4MHYsGEDp2KkcjpfKEWJjY2Fr68vtm/fjtKlS2Ps2LGYPn06LCwsREcjIhXr1asX0tLScP78edFRSAs9f/4czs7OqFmzJiIjI1G6dGnRkUgH6MShHHVkb2+PTZs2ISkpCVOmTEFgYCDq16+PIUOG4OrVq6LjEZEKSaVSXLhwAQ8fPhQdhbRMfn4+BgwYgNzcXISEhLBMUolhoRSsVq1aWLx4MWQyGfz8/HDp0iU4Ozujffv2OHDgAORyueiIRKRkvXr1gqGhIUJCQkRHIS3j5eWF8+fPIyQkBHXq1BEdh3QIC6WaKFeuHKZOnYq4uDgEBwcjOzsbPXv2hJ2dHdatW4ecnBzREYlISSpUqIAuXbrwIeekVIGBgVi+fDmWL1+OVq1aiY5DOoaFUs3o6+ujf//+uHjxIs6dO4eGDRtiwoQJsLCwwKJFi5Camio6IhEpgVQqxfnz5/Ho0SPRUUgLXLlyBd988w3GjBmDb775RnQc0kE8lKMB4uLi4O/vj02bNkGhUGDkyJHw8PCAjY2N6GhE9JlevXoFU1NTeHt7Y9q0aaLjkAZ79uwZnJ2dUbt2bURGRsLIyEh0JNJBLJQaJDU1FatXr8aKFSvw4sUL9OrVC15eXmjZsiUfCUGkgXr06IHXr1/j7NmzoqOQhsrPz0fHjh0RFxeH6Oho1K5dW3Qk0lEceWuQqlWr4vvvv0dycjLWrVuH+/fvo3Xr1nBzc8Pu3btRUFAgOiIRfYJ3Y+/Hjx+LjkIaysPDA5cuXUJISAjLJAnFQqmBSpcujTFjxuDmzZs4ePAgjI2NMXDgQFhbW2P58uV4+/at6IhE9BF69+4NAwMDnvamz7Jp0yasXLkSK1asQIsWLUTHIR3HkbeWuHbtGnx8fLBr1y6YmJhgwoQJmDp1Kn9iJVJzX331FTIyMnDmzBnRUUiDREVFoXXr1hgxYgTWrl0rOg4RC6W2kclkWL58OdauXYucnBwMHjwYnp6eaNKkiehoRPQXAgMDMXr0aDx69Ag1a9YUHYc0wNOnT+Hs7Axzc3NERETwEA6pBY68tYy5uTm8vb2RkpKCX375BREREXBwcEDXrl1x7Ngx8OcHIvXCsTd9iry8PPTv3x9yuRwhISEsk6Q2WCi1VIUKFeDp6YmEhARs374dqamp6Nq1KxwcHLB582bk5eWJjkhEACpVqoROnTrxIef0UWbMmIGoqCiEhIRwRZvUCgullitVqhSGDBmC6OhonDp1Cubm5hg5ciTq1q2LJUuW4OXLl6IjEuk8qVSKs2fP4smTJ6KjkBrbsGEDVq1ahZUrV6J58+ai4xB9gHsoddDt27fh5+eHLVu2oFSpUhgzZgxmzJiBunXrio5GpJPS09NRvXp1+Pv7Y/LkyaLjkBq6dOkS2rZti9GjR2PVqlWi4xAVwUKpw549e4aVK1di5cqVePXqFdzd3eHp6QlXV1fR0Yh0Trdu3ZCTk4PTp0+LjkJq5smTJ3ByckK9evVw6tQpGBoaio5EVARH3jqsevXqWLRoEVJSUhAQEIA//vgDbm5uaN26NcLDw1FYWCg6IpHOkEqlOHPmDJ4+fSo6CqmRd4dwJBIJ9uzZwzJJaouFklC2bFlMnDgRd+/eRVhYGBQKBfr27YuGDRti1apVyMrKEh2RSOv16dMHenp6CA0NFR2F1Mi0adMQHR2N0NBQ1KhRQ3Qcor/FkTf9pUuXLsHHxwehoaGoVKkSJk+ejMmTJ8PU1FR0NCKt1bVrV+Tl5SEiIkJ0FFIDa9euxYQJE7BhwwaMHj1adByif8RCSf8oMTER/v7+2LhxIwoKCjB8+HDMnDkTtra2oqMRaZ3169djwoQJePz4MapXry46Dgl04cIFtGvXDuPGjcPKlStFxyH6VyyU9FHS09OxZs0arFixAk+ePMFXX30FLy8vtG3bFhKJRHQ8Iq2QmpqKGjVqYMWKFZg4caLoOCTI48eP4eTkBCsrK5w8eZL7JkkjsFDSJ8nNzcWOHTvg4+ODmzdv4osvvoCXlxf69++PUqVKiY5HpPG6dOmCgoICnDp1SnQUEiA3Nxft2rVDSkoKrl69ypVq0hg8lEOfxMjICCNHjsSNGzdw5MgRVKlSBUOGDEH9+vXh6+uLN2/eiI5IpNGkUikiIyPx/Plz0VGohCkUCkyZMgV//PEHwsLCWCZJo7BQ0meRSCTv3w9+/fp1tG/fHrNnz4aZmRlmzZqFlJQU0RGJNFLfvn0hkUh42lsHrVmzBuvXr8fq1avh4uIiOg7RJ+HIm5Tm0aNHWLFiBVavXo3MzEwMHDgQnp6eaNq0qehoRBqlc+fOkMvlOHnypOgoVELOnz+P9u3bY8KECVixYoXoOESfjIWSlC4jIwMbN26En58fkpOT0aFDB3h5eaFbt248wEP0EdauXYuJEyfiyZMnfFSXDnj06BGcnJzQoEEDnDhxgvvRSSNx5E1KZ2JigunTpyM+Ph67du1CRkYGvvzySzRu3BgbNmxATk6O6IhEau3d2DssLEx0FFKxnJwc9OvXD4aGhggODmaZJI3FQkkqY2BggAEDBuDy5cs4c+YMrKysMG7cOFhaWuLnn39GWlqa6IhEaqlatWpo164dgoODRUchFVIoFJg8eTJiYmIQGhrK1WjSaCyUpHISiQStW7fG3r17cefOHfTp0wf/+c9/YGZmhsmTJyM+Pl50RCK1I5VKERERgRcvXoiOQiqyatUqbNy4EWvXroWzs7PoOETFwkJJJapBgwZYvXo1ZDIZZs+ejeDgYNjY2KBfv364cOGC6HhEaqNv374AwLG3ljp79iymT5+O6dOnY/jw4aLjEBUbD+WQUNnZ2di6dSt8fX1x7949NG/eHJ6enujTpw/09fVFxyMSqmPHjtDT08Px48dFRyElSklJgbOzMxo1aoRjx45x3yRpBa5QklBlypTB+PHjcfv2bezbtw+Ghobo378/bGxsEBAQgMzMTNERiYR5N/ZOTU0VHYWUJCcnB+7u7ihdujR2797NMklag4WS1IKenh569uyJ06dP48qVK2jWrBlmzJgBMzMzzJs3D0+ePBEdkajE9evXDwqFgmNvLaFQKDBx4kTExsYiLCwM1apVEx2JSGk48ia1lZycjGXLlmHdunXIy8vD0KFDMXPmTDRu3Fh0NKIS06FDB5QqVQpHjx4VHYWKKSAgAFOnTsXWrVsxbNgw0XGIlIorlKS2LCws4Ovri5SUFPz88884duwY7O3t0b17d5w4cQL8WYh0gVQqxcmTJ/mYLQ0XGRmJGTNmwMPDg2WStBILJam9ihUrYtasWUhMTMTWrVvx5MkTdO7cGU2bNsW2bduQn58vOiKRyrwbe4eHh4uOQp9JJpNBKpWibdu2WLp0qeg4RCrBkTdpHIVCgVOnTsHb2xtHjhxB7dq1MW3aNIwfPx4VK1YUHY9I6dq3bw8jIyMcOXJEdBT6RNnZ2WjdujVSU1MRHR2NqlWrio5EpBJcoSSNI5FI0LFjRxw+fBixsbHo0qUL5s+fDzMzM8ycORPJycmiIxIp1buxd3p6uugo9AkUCgW++eYb3L59G2FhYSyTpNVYKEmjNW7cGBs3bkRSUhKmTZuGwMBA1K9fH4MGDUJ0dLToeERK0a9fPxQWFnLsrWGWL1+OLVu2YMOGDWjatKnoOEQqxZE3aZXMzExs2rQJfn5+SExMRJs2beDl5YWvvvoKenr8+Yk0V7t27VCmTBkcPnxYdBT6CBEREejcuTNmzJgBb29v0XGIVI7fYUmrGBsbY8qUKbh//z727NmD/Px89OrVC40aNcLatWuRnZ0tOiLRZ5FKpThx4gTH3hogOTkZAwYMQPv27bFkyRLRcYhKBAslaSV9fX24u7vjwoULOH/+POzs7PDNN9/AwsICP/74I168eCE6ItEncXd3R2FhIfbu3Ss6Cv2DrKws9O3bFyYmJti5cycMDAxERyIqERx5k86Ij4+Hv78/Nm3aBLlcjhEjRsDDwwMNGjQQHY3oo7Rt2xbGxsY4dOiQ6Cj0FxQKBb7++muEhobi4sWLcHBwEB2JqMRwhZJ0hpWVFQICAiCTyfD9998jPDwcDRs2RO/evXHmzBk+KJ3U3rux98uXL0VHob/g7++P7du3Y+PGjSyTpHNYKEnnVKlSBfPmzUNSUhLWr1+P+Ph4tG3bFq6urti9ezcKCgpERyT6S+7u7igoKODYWw2dPHkSs2bNwrfffotBgwaJjkNU4jjyJp0nl8tx9OhReHt749SpU7CwsMCMGTMwZswYmJiYiI5H9IE2bdrAxMQEBw8eFB2F/k9SUhKcnZ3h5OSEQ4cOQV9fX3QkohLHFUrSeXp6eujevTtOnjyJa9euoVWrVpg1axbMzMwwe/ZsPHr0SHREovekUimOHz+OV69eiY5C+P+HcCpUqIAdO3awTJLOYqEk+pN37wdPTEzEuHHjsHr1alhaWmLEiBG4ceOG6HhEHHurEYVCgbFjx+L+/fsIDw9H5cqVRUciEoaFkugvmJmZ4bfffkNKSgqWLl2K06dPw8HBAV26dMHRo0d5gIeEqVWrFlq2bIng4GDRUXSer68vduzYgcDAQNjb24uOQyQUCyXRPyhfvjw8PDwQHx+PoKAgpKWloVu3bmjSpAkCAwORm5srOiLpIKlUimPHjnHsLdDx48fx7bffYs6cOZBKpaLjEAnHQzlEn0ChUCAyMhI+Pj44cOAAatasialTp2LChAkcd1GJefToEerUqYPNmzdj+PDhouPonAcPHsDZ2RnNmjXDgQMHuG+SCCyURJ/tzp078PPzw5YtW6Cvr48xY8ZgxowZqFevnuhopANatWqFSpUqYf/+/aKj6JTMzEy0aNECmZmZuHLlCipVqiQ6EpFa4Mib6DM1bNgQa9euRXJyMry8vBAUFARra2tIpVJcvnxZdDzScu/G3q9fvxYdRWcoFAqMGTMGCQkJCA8PZ5kk+hMWSqJiql69On788UfIZDKsXLkSMTExcHNzQ6tWrRAWFobCwkLREUkLubu7Iy8vD/v27RMdRWf89ttv2LVrFzZv3ozGjRuLjkOkVlgoiZSkbNmy+Oabb3D37l2Eh4dDIpGgX79+sLW1xapVq5CVlSU6ImmROnXqoEWLFjztXUKOHTuG7777DnPnzoW7u7voOERqh3soiVQoKioKPj4+2LNnDypVqoSJEydiypQpqF69uuhopAX8/f0xe/ZsPH/+HBUqVBAdR2slJCTAxcUFbm5u2L9/Pw/hEP0FFkqiEvDgwQP4+/tjw4YNKCgowLBhwzBz5kw0atRIdDTSYA8fPoSZmRm2bt2KYcOGiY6jld6+fYsWLVogJycHUVFRqFixouhIRGqJhZKoBL18+RJr1qzB8uXL8eTJE3z11Vfw9PREu3btIJFIRMcjDdSiRQtUq1aNb85RAYVCgQEDBuDIkSO4fPkyfwAk+gfcQ0lUgipVqoQ5c+YgKSkJgYGBkMlk6NChA5ydnREUFIT8/HzREUnDSKVSHD16FG/evBEdRev8+uuv2LNnD7Zs2cIySfQvWCiJBDA0NMSIESMQExODo0ePomrVqhg6dCjq168PHx8flgP6aP3790dubi6fR6lkR44cwdy5czF//nz07dtXdBwitceRN5GauHHjBnx9fREUFIQyZcpg3LhxmD59OszMzERHIzXXvHlzVK9eHeHh4aKjaIX4+Hi4uLigVatW2Lt3L/T0uPZC9G9YKInUzOPHj7FixQqsXr0aGRkZGDhwIDw9PfHFF1+IjkZqytfXF3PnzsXz589Rvnx50XE0WkZGBpo3b478/HxERUXx9DzRR+KPXURqplatWvjll1+QkpICX19fXLhwAU5OTujQoQMOHjwIuVwuOiKpmXdj7wMHDoiOotEUCgVGjRoFmUyG8PBwlkmiT8BCSaSmypUrh2nTpiEuLg67d+9GZmYmevTogcaNG2P9+vXIyckRHZHUhLm5OVxdXfmQ82L65ZdfEBISgq1bt6Jhw4ai4xBpFBZKIjVnYGAAqVSKS5cu4ezZs7CxscH48eNhYWGBn3/+GWlpaaIjkhqQSqU4fPgwMjIyREfRSAcPHsT333+PhQsXonfv3qLjEGkc7qEk0kD379+Hn58fAgMDIZFIMHLkSHh4eMDa2lp0NBIkOTkZlpaWCAoKwuDBg0XH0ShxcXFwcXFB27ZtERYWxkM4RJ+BhZJIg7148QKrVq1CQEAAUlNT0bt3b3h5eaFFixZ8ULoOcnV1Re3atREaGio6isbIyMiAm5sbCgsLERUVxUNNRJ+JP4YRabBq1aphwYIFkMlkWLt2Le7evYtWrVqhefPm2LNnDwoLC0VHpBL0buz99u1b0VE0glwux/Dhw/Hw4UPs3buXZZKoGFgoibRA6dKlMXbsWNy6dQv79+9HmTJlIJVKYW1tjRUrVrBg6Ij+/fsjJyeHp70/0n/+8x+Eh4dj27ZtaNCggeg4RBqNI28iLXX16lX4+Phg9+7dKF++PL755htMnToVNWvWFB2NVKhZs2YwMzNDSEiI6Chq7cCBA+jVqxd++OEHLFiwQHQcIo3HQkmk5WQyGZYtW4Z169YhJycHQ4cOhaenJxo3biw6GqnAb7/9hgULFuDFixcoV66c6Dhq6d69e2jWrBk6dOiAkJAQHsIhUgIWSiId8fr1a6xbtw7Lli3Dw4cP0bVrV3h5eaFjx448wKNFHjx4gHr16mHnzp0YOHCg6Dhq582bN3B1dYVEIsGlS5e4b5JISfhjGZGOqFChAry8vJCYmIht27bh2bNn6Ny5MxwdHbFlyxbk5eWJjkhKULduXTg7O/Mh539BLpfj66+/xpMnTxAeHs4ySaRELJREOqZUqVIYOnQorl27hpMnT6JOnToYMWIE6tati19//RWvXr0SHZGKSSqV4tChQ8jMzBQdRa389NNP2L9/P7Zv3w4bGxvRcYi0CgslkY6SSCTv3w9+69YtdO/eHQsWLICZmRlmzJiBpKQk0RHpM0mlUmRnZ+PgwYOio6iNffv24YcffsCiRYvw1VdfiY5DpHW4h5KI3nv69CkCAgKwatUqvHr1Cv3794eXlxdcXFxER6NP5OzsjLp163L0DeDu3bto1qwZOnfujODgYB7CIVIB/q0iovdq1KiBn3/+GTKZDMuXL8fVq1fRrFkztG3bFvv27YNcLhcdkT6SVCrFwYMHdX7s/fr1a/Tu3RtmZmYIDAxkmSRSEf7NIqIijI2NMXnyZNy7dw+hoaEoKChA79690bBhQ6xZswbZ2dmiI9K/eDf2PnTokOgowsjlcgwbNgzPnj1DeHg4TExMREci0loslET0t/T19dG3b1+cP38eFy5cgL29PSZNmgRzc3P88MMPeP78ueiI9Dfq1auHL774QqdH3j/++CMOHjyIHTt2wNraWnQcIq3GQklEH+Xd+8Hv37+PQYMG4bfffoO5uTkmTJiAe/fuiY5Hf+Hd2DsrK0t0lBIXHh6ORYsW4T//+Q+6d+8uOg6R1uOhHCL6LOnp6Vi9ejWWL1+OZ8+eoWfPnvD09ESbNm34oHQ1kZCQACsrKwQHB6N///6i45SY27dvw9XVFd26dcPu3bv555GoBLBQElGx5ObmIigoCD4+Prh16xacnZ3h5eUFd3d3GBgYiI6n87744gtYW1tj165doqOUiFevXqFZs2YwMjLCxYsX+fpJohLCkTcRFYuRkRFGjRqF2NhYHD58GBUqVMCgQYNgZWUFPz8/ZGRkiI6o06RSKQ4cOKATY+/CwkIMHToUL168QHh4OMskUQlioSQipZBIJOjWrRtOnDiBP/74A61bt8a3334LMzMzfPvtt3j48KHoiDpJKpUiKysLhw8fFh1F5RYuXIgjR45g586dqF+/vug4RDqFI28iUpmHDx9ixYoVWL16NbKysjBo0CB4enrC0dFRdDSd0rRpUzRo0AA7d+4UHUVlQkND4e7ujiVLlmD27Nmi4xDpHBZKIlK5jIwMbNiwAX5+fpDJZOjYsSO8vLzQtWtXHpgoAYsXL8bixYvx4sULlClTRnQcpbt16xZcXV3x1VdfYefOnfwzRSQAR95EpHImJiaYMWMGEhISsHPnTrx+/Rrdu3dHkyZNsGnTJuTm5oqOqNWkUikyMzO1cuz98uVL9OnTB/Xq1cPGjRtZJokEYaEkohJjYGCAgQMHIioqCpGRkahbty5Gjx4NS0tLLF68GOnp6aIjaiVra2s4ODho3UPOCwsLMWTIEKSlpSE8PBzGxsaiIxHpLBZKIipxEokEbdq0wb59+3Dnzh306tULixYtgpmZGaZOnYqEhATREbWOVCrF/v37teq1mfPnz8exY8ewa9cu1KtXT3QcIp3GQklEQtna2mLNmjWQyWSYNWsWdu7cCRsbG/Tv3x8XL14UHU9rvBt7HzlyRHQUpdizZw9++eUXLFmyBJ07dxYdh0jn8VAOEamV7OxsbNmyBb6+vrh//z5atGgBT09P9O7dG/r6+qLjaTQHBwfY2dkhKChIdJRiiY2NRfPmzdGzZ08EBQVx3ySRGuAKJRGplTJlymDChAm4c+cO9u7dCwMDA7i7u6NBgwZYuXIlMjMzRUfUWNow9k5PT0efPn1gZWWFDRs2sEwSqQkWSiJSS3p6eujVqxciIyMRFRUFZ2dnTJs2Debm5vj+++/x9OlT0RE1jlQqxdu3b3H06FHRUT5LYWEhBg8ejFevXiEsLAxly5YVHYmI/g8LJRGpPRcXF+zcuRMJCQkYPnw4/P39YWFhgbFjx+L27dui42mMBg0awN7eXmNPe8+bNw8nTpzArl27ULduXdFxiOhPWCiJSGNYWlrCz88PDx8+xE8//YTDhw/Dzs4OX375JU6dOgVuCf9378beOTk5oqN8kt27d+PXX3/Fb7/9hk6dOomOQ0T/g4WSiDROxYoV8e233+LBgwfYvHkzHj16hI4dO+KLL77A9u3bkZ+fLzqi2pJKpcjIyNCosXdMTAxGjRqFoUOHwsPDQ3QcIvoLPOVNRBpPoVDg5MmT8Pb2xtGjR1GnTh1Mnz4d48aNQ4UKFUTHUzv29vZwcHDAtm3bREf5V2lpaXBxcUGFChVw/vx57pskUlNcoSQijSeRSNCpUyccOXIEN27cQKdOnTB37lyYmZnB09MTMplMdES1IpVKsW/fPrUfexcUFGDw4MF48+YND+EQqTkWSiLSKvb29ti0aROSkpIwZcoUbNq0CfXq1cOQIUNw9epV0fHUwrux97Fjx0RH+Udz587FqVOnsHv3blhaWoqOQ0T/gIWSiLRSrVq1sHjxYshkMvj5+eHSpUtwdnZG+/btceDAAcjlctERhWnYsCHs7OzU+rT3jh078Ntvv8Hb2xsdOnQQHYeI/gULJRFptXLlymHq1KmIi4tDcHAwsrOz0bNnT9jZ2WH9+vVqP/ZVlXdj79zcXNFRirh+/TrGjBmDYcOGYfr06aLjENFH4KEcItIpCoUCFy5cgI+PD8LDw1GtWjVMnjwZkyZNQtWqVUXHKzG3b9+GnZ0d9u3bh549e4qO815qaipcXFxQuXJlnDt3DmXKlBEdiYg+AgslEemsuLg4+Pv7Y9OmTQCAESNGwMPDAzY2NoKTlQw7Ozs4OTlhy5YtoqMA+O8hnK5duyI2NhbR0dEwNzcXHYmIPhJH3kSks6ytrbFy5UrIZDLMnTsXoaGhsLW1Rd++fXHu3Dmtf1C6VCrF3r171WbsPXv2bERGRmL37t0sk0QahoWSiHRe1apV8f333yM5ORnr1q3DvXv30Lp1azRv3hzBwcEoKCgQHVElpFIp3rx5g+PHj4uOgqCgIPj6+sLX1xft2rUTHYeIPhFH3kRE/0Mul+PIkSPw9vZGREQELC0t4eHhgdGjR6NcuXKi4ylVo0aN4OLigs2bNwvLcO3aNbRs2RIDBw7Epk2bIJFIhGUhos/DQklE9A+uXbsGHx8f7Nq1CyYmJvjmm28wdepU1KpVS3Q0pVi4cCGWLVuGZ8+ewcjIqMTv/+LFCzg7O8PU1BRnz55F6dKlSzwDERUfR95ERP/g3fvBExMTMWbMGKxcuRKWlpYYOXIkYmNjRccrNqlUitevX+PEiRMlfu+CggIMHDgQOTk5CA0NZZkk0mAslEREH8Hc3Bze3t5ISUnBL7/8glOnTqFJkybo2rUrjh07prEHeOzs7GBrayvkIeezZs3C2bNnERwcDDMzsxK/PxEpDwslEdEnqFChAjw9PZGQkIDt27cjNTUVXbt2hYODAzZv3oy8vDzRET+JRCJ5f9q7JLNv3boV/v7+8Pf3R5s2bUrsvkSkGiyURESfoVSpUhgyZAiio6Nx6tQpmJubY+TIkahbty5+/fVXvHz5UnTEjyaVSvHq1asSG3tfvXoV48ePx6hRozBp0qQSuScRqRYP5RARKcnt27fh5+eHLVu2oFSpUhgzZgxmzJiBunXrio72jxQKBRo2bIjmzZu/f8i7qjx//hzOzs6oWbMmIiMjuW+SSEtwhZKISEkaNWqEdevWQSaTYebMmdi2bRusrKwwYMAAREVFiY73t96NvcPDw1U69s7Pz8eAAQOQm5uLkJAQlkkiLcJCSUSkZNWrV8eiRYuQkpKCgIAA/PHHH3B1dUXr1q2xd+9eyOVy0RGLeDf2PnnypMru4eXlhfPnzyMkJAR16tRR2X2IqOSxUBIRqUjZsmUxceJE3L17F2FhYVAoFOjTpw9sbW2xevVqZGVliY74nr29PWxsbFR22jswMBDLly/H8uXL0apVK5Xcg4jEYaEkIlIxfX199OnTB+fOncPFixfh4OCAyZMnw9zcHAsXLsTz589FR/xg7J2fn6/Ua1+5cgXffPMNxowZg2+++Uap1yYi9cBDOUREAiQmJsLf3x8bN25EQUEBhg8fjpkzZ8LW1lZYppiYGDg6OuLw4cPo1q2bUq757NkzODs7o3bt2oiMjBTyNh4iUj0WSiIigdLT07FmzRqsWLECT548QY8ePeDl5YU2bdqU+DutFQoFGjRogNatW2PDhg3Fvl5+fj46duyIuLg4REdHo3bt2kpISUTqiCNvIiKBKleujO+++w4PHjzApk2bkJSUhHbt2sHFxQU7duxQ+vj5nyh77O3h4YFLly4hJCSEZZJIy7FQEhGpASMjI4wcORI3btzAkSNHULlyZQwZMgRWVlbw9fXFmzdvSiSHVCpFeno6Tp06VazrbNq0CStXrsSKFSvQokULJaUjInXFkTcRkZqKiYmBr68vgoKCULZsWYwfPx7Tpk1T6XuvFQoFbGxs0LZtW6xfv/6zrhEVFYXWrVtjxIgRWLt2rZITEpE6YqEkIlJzjx49wooVK7B69WpkZmZi4MCB8PT0RNOmTVVyv++++w5r167F06dPUapUqU/6vU+fPoWzszPMzc0RERHBQzhEOoIjbyIiNVe7dm0sWbIEKSkp8Pb2xrlz5/DFF1+gY8eOOHz4MJS9LvBu7B0REQEAyMwtwK3Hr/GH7CVuPX6NzNyCv/x9eXl56N+/P+RyOUJCQlgmiXQIVyiJiDRMQUEBQkND4e3tjStXrqBRo0bw9PTE0KFDlVLiFAoF6n/RCrXbDoRebXvI0rPw528UEgDmlcuifQNTDHU1h3V1EwDApEmTsH79ekRGRqJ58+bFzkFEmoOFkohIQykUCpw7dw7e3t7Yv38/TE1NMXXqVHzzzTeoUqXKZ10zJT0Lc8NicTY+FQp5ISR6+n/7a/X1JCiUK9Daqiqa5N3Gt5NGY926dRg7duznfkpEpKFYKImItMC9e/fg5+eHzZs3Q09PD6NGjYKHhwfq16//0dfYeUWGhftuoUCuQKH847816EmAwrxc2OffxYFlcz8nPhFpOBZKIiIt8uLFC/z+++9YuXIlUlNT0bdvX3h6ev7ro3sCIuLgfez+599YoQAkEnh1scGU9taffx0i0kgslEREWig7Oxtbt26Fr68v7t27h+bNm8PLywu9e/eGvv6HY+ydV2SYExqrtHv/2s8eA13MlXY9IlJ/POVNRKSFypQpg/Hjx+P27dvYt28fDA0N4e7ujgYNGiAgIACZmZkA/rtncuG+W0q994J9t5CSnqXUaxKReuMKJRGRjoiOjoaPjw+Cg4NRvnx5TJw4EXG1u+Lqo7eftGfy3+jrSdCiXhVsHeOqtGsSkXpjoSQi0jHJyclYtmwZNu45hIpDfVR2nxMebWBlaqKy6xOR+mChJCLSUd/tuYadVx9DAYnSr62vJ8HXrhb4oZed0q9NROqHeyiJiHTUhQevVVImAaBQrkDE/ecquTYRqR8WSiIiHfQ2twAyFR+ckaVl/e1rGolIu7BQEhHpoOS0TKh6v5MCQFJaporvQkTqgIWSiEgH5RXIteo+RCQWCyURkQ4yNCiZL/8ldR8iEot/04mIdJBlFWMVHcf5/yT/dx8i0n4slEREOsjYyADmlcuq9B7mVcrC2MhApfcgIvXAQklEpKPaNzCFvp5q1in19SRob2OqkmsTkfphoSQi0lFDXc2V+srFPyuUKzDMzVwl1yYi9cNCSUSko6yrm6C1VVWlr1Lq60nQ2qoqX7tIpENYKImIdNjivvYwUHKhNNCTYHFfe6Vek4jUGwslEZEOM6tcFj8q+X3bi3rZwUzFB36ISL2wUBIR6bhBLubw6mKjlGvN6tIAA124d5JI10gUCoWq375FREQaYOcVGRbuu4UCueKTDuvo60lgoCfBol52LJNEOoqFkoiI3ktJz8LcsFicjU+Fvp7kH4vlu4+3tqqKxX3tOeYm0mEslEREVETcswxsvyxDxP3nkKVl4c/fKCT470PL29uYYpibOU9zExELJRER/bPM3AIkpWUir0AOQwM9WFYx5htwiOgDLJREREREVCw85U1ERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExcJCSURERETFwkJJRERERMXCQklERERExfL/AFk3Iy1YFUx6AAAAAElFTkSuQmCC",
      "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()"
   ]
  },
  {
   "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": [
       "np.float64(2.23606797749979)"
      ]
     },
     "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": [
       "np.float64(2.23606797749979)"
      ]
     },
     "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": [
       "np.float64(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();"
   ]
  },
  {
   "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.15, 0.59, 0.15, 0.48],\n",
       "       [0.46, 0.  , 0.89, 0.24, 0.41],\n",
       "       [0.74, 0.6 , 0.  , 0.53, 0.64],\n",
       "       [0.9 , 0.54, 0.84, 0.  , 0.82],\n",
       "       [0.94, 0.61, 0.94, 0.56, 0.  ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_estimator = FrechetMean(space=total_space, method=\"aac\", max_iter=20)\n",
    "\n",
    "fm = mean_estimator.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();"
   ]
  },
  {
   "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 = GGPCA(space=total_space, 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();"
   ]
  },
  {
   "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 = GeneralizedGeodesicRegression(space=total_space)"
   ]
  },
  {
   "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.18338109],\n",
       "       [-0.16045845],\n",
       "       [-0.02578797],\n",
       "       [ 0.00286533],\n",
       "       [ 0.04011461],\n",
       "       [ 0.        ],\n",
       "       [ 0.02005731],\n",
       "       [-0.18624642],\n",
       "       [ 0.06017192],\n",
       "       [-0.01432665],\n",
       "       [ 0.        ],\n",
       "       [ 0.        ],\n",
       "       [ 0.        ],\n",
       "       [-0.12034384],\n",
       "       [-0.08309456],\n",
       "       [ 0.00286533],\n",
       "       [-0.0487106 ],\n",
       "       [ 0.        ],\n",
       "       [ 0.03438395],\n",
       "       [ 0.16332378],\n",
       "       [ 0.02005731],\n",
       "       [-0.13753582],\n",
       "       [ 0.17765043],\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": [
       "np.float64(16.34454451141818)"
      ]
     },
     "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": "geomstats",
   "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}