{ "cells": [ { "cell_type": "markdown", "id": "ada6ec81", "metadata": {}, "source": [ "# Stratified spaces\n", "\n", "$\\textbf{Lead Author: Anna Calissano}$" ] }, { "cell_type": "markdown", "id": "3af12b42", "metadata": {}, "source": [ "Dear learner, \n", "the aim of the current notebook is to introduce stratified spaces and its implementation within geomstats. " ] }, { "cell_type": "markdown", "id": "4c0d7eda", "metadata": {}, "source": [ "## Spider" ] }, { "cell_type": "markdown", "id": "5c26e954", "metadata": {}, "source": [ "The $k$-Spider consists of $k$ copies of the positive real line $\\mathbb{R}_{\\geq 0}$ glued together at the origin. Within geomstats, we defined the following:\n", "1. Spider Point: a point object defining the ray and the value\n", "2. Spider: the space defined by the number of rays\n", "3. Spider Geometry: by chosing a metric on the rays, we can define a metric on the whole space" ] }, { "cell_type": "markdown", "id": "eb5abd16", "metadata": {}, "source": [ "" ] }, { "cell_type": "code", "execution_count": 1, "id": "72158cdb", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO: Using numpy backend\n" ] } ], "source": [ "import geomstats.backend as gs\n", "\n", "from geomstats.geometry.stratified.spider import Spider\n", "\n", "gs.random.seed(2020)" ] }, { "cell_type": "markdown", "id": "6dc7f1cd", "metadata": {}, "source": [ "We can define a spider with $k=3$ rays (strata) and sample two points from it." ] }, { "cell_type": "code", "execution_count": 2, "id": "85659a76", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spider = Spider(n_rays=3, equip=True)\n", "\n", "spider.n_rays" ] }, { "cell_type": "code", "execution_count": 3, "id": "475990f8", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[r0: 10.028180271833065, r0: 11.079704435029091]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spider_points = spider.random_point(n_samples=2)\n", "\n", "spider_points" ] }, { "cell_type": "markdown", "id": "ac7750b6", "metadata": {}, "source": [ "The points are represented into the SpiderPoint format, where the first input is the stratum and the second input is the value along the stratum." ] }, { "cell_type": "code", "execution_count": 4, "id": "f59ca8e4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n", "[10.02818027]\n" ] } ], "source": [ "print(spider_points[0].stratum)\n", "print(spider_points[0].coord)" ] }, { "cell_type": "markdown", "id": "b62044c5", "metadata": {}, "source": [ "Given a metric $d_{rays}$ on the strata (the rays), we can extend it to the whole space by $$d_{Spider}(s_1,s_2)=d_{rays}(s_1,0) + d_{rays}(0,s_2)$$" ] }, { "cell_type": "markdown", "id": "aeb317b1", "metadata": {}, "source": [ "Given two points on the Spider, we can compute the distance between the two points as well as the geodesic between the two." ] }, { "cell_type": "code", "execution_count": 5, "id": "a6126dd3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(1.05152416)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spider.metric.dist(spider_points[0], spider_points[1])" ] }, { "cell_type": "code", "execution_count": 6, "id": "7881d1da", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[r0: 10.028180271833065] [r0: 10.55394235343108] [r0: 11.079704435029091]\n" ] } ], "source": [ "spider_geodesic_func = spider.metric.geodesic(spider_points[0], spider_points[1])\n", "\n", "print(spider_geodesic_func(0), spider_geodesic_func(0.5), spider_geodesic_func(1))" ] }, { "cell_type": "markdown", "id": "76023d66", "metadata": {}, "source": [ "## Graph Space" ] }, { "cell_type": "markdown", "id": "3358a89e", "metadata": {}, "source": [ "Graph Space is a space defined to describe set of graphs with a finite number of nodes which can be both node labelled or node unlabelled. \n", "\n", "Inspired by: A. Calissano, A. Feragen, S. Vantini, Populations of unlabeled networks: Graph space geometry and geodesic principal components, MOX Report (2020)\n", "\n", "\n", "We consider graphs as triples $G=(V,E,a)$, where the node set $V$ has at most $n$ elements, and the edge set $E \\subset V^2$ has maximal size \n", "$n^2$. The nodes and edges are attributed with elements of an attribute space $A$, which is considered to be Euclidean, via an attribute \n", "map $a \\colon E \\rightarrow A$. Here, the map $a$ allows us to describe attributes on both edges and nodes, as we use self loop edges (diagonal \n", "elements in the graphs adjacency matrix) to assign attributes to nodes. \n", "A graph with scalar attributes is completely specified by a weighted adjacency matrix of dimension $n\\times n$, residing in a space \n", "$X=\\mathbb{R}^{n^2}$ of flattened adjacency matrices. If the attributes are vectors of dimension $d$, the graph is represented by a tensor of \n", "dimension $n\\times n\\times d$, residing in a space $X=\\mathbb{R}^{n\\times n\\times d}$." ] }, { "cell_type": "code", "execution_count": 7, "id": "8452b3d8", "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "import matplotlib.pyplot as plt\n", "\n", "from geomstats.geometry.stratified.graph_space import GraphSpace" ] }, { "cell_type": "markdown", "id": "fa02c396", "metadata": {}, "source": [ "### Graph\n", "Consider a graph with $n=3$ nodes and $A=\\mathbb{R}$ scalar attributes on nodes and edges. It is represented by its adjacency matrix." ] }, { "cell_type": "code", "execution_count": 8, "id": "fd8f3042", "metadata": { "scrolled": true }, "outputs": [], "source": [ "graph_point = gs.array([[10, 3, 1], [3, 2, 4], [1, 4, 5]])" ] }, { "cell_type": "markdown", "id": "0f89275b", "metadata": {}, "source": [ "To simplify the visualization and the access to different methods, the graph can be turned into a networkx graph." ] }, { "cell_type": "code", "execution_count": 9, "id": "edee3873", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8nElEQVR4nO3de3gkd33n+09VdfVNLY1Gc5PGY+wZj4094wsXYwMJJE5ir21ICAs5BE7ObvAlYE6AnE1CloRLsrGTzXI2bA4+BxMmgVxMNoGQJxvA2QyLQ4wBmxADgwfH2B4b25LmPqPWpW9Vdf5o9UiaUZWq1V3dVV3v1/P48WNVq1W2ZzQf/X71+X0Nz/M8AQAAAOtk9vsGAAAAkGwESgAAAHSEQAkAAICOECgBAADQEQIlAAAAOkKgBAAAQEcIlAAAAOgIgRIAAAAdIVACAACgIwRKAAAAdIRACQAAgI4QKAEAANARAiUAAAA6QqAEAABARwiUAAAA6AiBEgAAAB0hUAIAAKAjBEoAAAB0hEAJAACAjhAoAQAA0BECJQAAADpCoAQAAEBHCJQAAADoCIESAAAAHSFQAgAAoCMESgAAAHQk0+8bQLLVHVfH5+s6Pl/Tibm6ZmuOZioNPX5s7sxrLtta0lDW0nDO0qZiVpuGbI0VbGUsfp4BAGAQGJ7nef2+CSTHifmavndkTs+cXNCxuZpOVxrrfq/RQkabi1ntHCvqsm1D2pC3u3inAACgVwiUWNPx+Zq+d3hWBw/P6vBsLbKvc95ITpdtK+myrSWNFgiXAAAkBYESvqZmKvpfTxzXoRMLPf/al2we0o9dvElbhrI9/9oAAKA9BEqc43SlrvufOKED0+W+3odhSC/ePqIf2TWmUo7HfQEAiCsCJc5wXE8PHDqhrz5zSo4bn18WWcvQq3eN6eUvGJVhGP2+HQAAcBYCJSRJlYajv/7OtJ7qw/Z2WJdtHdLr9m6TTTscAIBYIVBCpxbq+otvTenYXPuFG9OQNhbspeOAirZylqmG6+noXE3bSjmZhlRpuDqxeLzQ8bm6TlbqWs+vvO0jOb3pqgm2wAEAiBECZcpNl6v61COTmqs5oT9nrGhrz9aSLt06pK2lnCyz/W1ox/U0Va7qsSPN9ng7xw9tyGf0cy/ZrrEihR0AAOKAQJlic7WGPv7QsypX1w6Tlmno6h0jumpiRFtL2a4+y+h5nqZmqvrW5Iz+ZXIm1MrlWNHWbdecr1yG7W8AAPqNQJlSrufp3n+Z1NMn135m8vLxkq67aFNPzoY8NlfT/3riuB4/Orfmay/bOqQ3XDFOUQcAgD4jUKbUl544rgefPhn4mpFcRm+8clznbcj36K6WHDoxr88emNZ83Q183Q2XbNa1LxjtzU0BAIBVsV+YQs+drqwZJidGcrrlmh19CZOStHOsqFtedr42FYNXRb/4/WPrKhMBAIDuIVCm0D89dSLw+q6xgv79S8/TcJ+b1BuLtm552Q6dN5LzfY3rSV85FByOAQBAtAiUKTNdrurJ4/O+10fzGf3bK8Zjc9Zj3rb0xisnVLQt39d893BZpxbqPbwrAACwXDxSA3rmqwFb3ZYhvfHKcRUCwls/jOQz+rdXbJNf9cbzpK//4FQvbwkAACxDoEyR+Zqjg0dmfa+/7PxRTYz055nJtewcK+ry8WHf69+enFHdCS7wAACAaBAoU+SZUwu+ZzyahmLfln7lhaO+12pO8yxLAADQewTKFHn+dMX32uXjwxrJx3uc4dZSTrs3FX2vPxfw7wcAAKJDoEyR5075B65dY/5BLU52ESgBAIgdAmVKeJ6n6bL/lvCO0Xg+O3m2HQHnYgb9+wEAgOgQKFPC86S6u/oDlFnL0GjMt7tbtg37n0lZaVDKAQCgHwiUkGkYiZmHbSXjNgEASBUCJQAAADpCoEyLgJW9muOq4bMdHjfzdcf3GouXAAD0B4EyJUzD0IjPbG7Xk6ZnktGQfu60f/FmQ0KeAwUAYNAQKFNk+wb/QsuzCTly57lTC77XghrgAAAgOgTKFAkKXAemyvL8xujEhON6evSw/+jI8wiUAAD0BYEyRc4PCFyHZ2t66sR8D++mfY8entXpSsP3elLO0gQAYNDw0FmKbB/JazhnqVxdvdjywFMntWus2PERQrOz0hNPSNWqlMtJu3dLpVJHbynH9fTg0yd9r28ZymqsYHf2RQAAwLqwQpkipmnomvNHfa8/e7qiBw75h7YgBw9K73pXMzyOjEgvfrH08pc3/z4y0vz4u97VfN167P/+MR2bq/lef8UFo4k5SxMAgEFDoEyZl+7YoFzG/3/7l586oSePh9/6PnRIuuEGae9e6aMflZ58sjmVZznPa378ox9tvu6GG5qfF9Z3p8v6xrOnfa+P5DK6fHw4/BsCAICuIlCmTC5j6uodGwJf89ffmdZTIULlvn3Snj3S/fc3/7nh/3jjiuv339/8vH371r7fg4dn9XcHjwS+5uUXjMoyWZ0EAKBfCJQp9MoLR33PpJSkquPqU9+a1Def818VvOsu6fbbpUpl7SB5tkaj+Xm33958n9V4nqevHDqhvz4wHXjo+uYhWy89b6S9GwAAAF1leHE/KwaReO7Ugv7km89rrQE5V04M67qLNmlk2aHh+/Y1w2C37Nsn3Xrr0j+fnK9r//eP6V+PzgV+nm0ZuvVl52tLKdu9mwEAAG0jUKbYwz84pf/5+LE1X5cxDb3iglG94oKNmnzW1J49zRXGbsnnm2Wd8R2OHjh0Qt949vSaQVeSXn/5Np6dBAAgBgiUKeZ5nv720SM6MF0O9fqCbereXz9PBx7KynG698yiZXm64pqa3nzX86o6bqjPufb8DbrhhVu6dg8AAGD9eIYyxQzD0E/u2aoXbQ+3yvfME5a+9dVcV8OkJDmOoW99LafnnrJCvf4VF4zq+ks2d/UeAADA+hEoU84yDb32sq36sd2b1nztNz63QaYZzYK2aXl6+HPB7XPDkF5z6Rb9xMWbOXMSAIAYIVBChmHohy7cqDdeMS474Pid7z9UlOtGE+Rcx9D3Hy76Xs9nTL3lRdv1kjWOPAIAAL3H6EWccdm2ksZHcvrSE8d18PDsimvVeUMnpqIdbXhi0lZ1wVCusHIV9KqJYV23e5OGA446AgAA/cOf0FhhY8HWG64Y17XnL2j/94/rudPNOveJSVtS1NvMhk48b2tid3PE4q6xgn784s0aH85F/HUBAEAnCJRY1Y7Rgn7+6vP0+NE5/fNzp/VcI1z7ulNew9DFm4u65vxR7drkvwUOAADig2ODEMrXvuHoldeEa2F34uvfcHTt1dF/HQAA0D0ESoQyOyuNjEhR/moxDGlmRiqVovsaAACg+2h5I5RSSdq1K9qvcdFFhEkAAJKIQInQbr5ZykT01G0mI910UzTvDQAAosWWN0I7eFDauze69//2AemKvc2tbwAAkBysUCK0PXuk66/v/iqllZF+/CekSy6V6m60z2kCAIDuY4USbTl0qBksK5XuvWc+Lz3yHenCnc1/NiRlLVYqAQBIClYo0ZadO6WPfKS77/nhP1gKk5LkSao5ksuPOgAAJAKBEm277Tbpzjtb/9RZ6vut35Z+/pZzP06oBAAgOdjyxrrt2yf94js91euS64Tfn85kmn99+A9WD5Nns03J4kcfAABiiz+msW633SZ98r5TuvxlVUmSaQX/bNIq81x3XbMx/gu3h/s6dVeqO5R1AACIK1YosW61hqu/e/SwPEnPPpnR/s8U9chXcpp+zpK8pRVLw2geWn7TTdIdd0iXXbb0Hq63GBZDfD3TaK5WUtYBACBeCJRYtx+cXNDDPzh1zscX5g2dnLR1zXljyucN7d4dPAHH85rPS4b5hUgDHACA+Ilo7gnSYGpm9bODCkVPe66x9NLzw6U+w2iGxLq7dgmnVdaxreaKJQAA6D+eocS6uJ6n6Zmq7/WJkXxb72csbmdbIUIiDXAAAOKFQIl1OTZbU90n0VmGtHU41/Z7GkZz5dEO+auy5kiO2/aXAQAAXUagxLpM+mx3S80wmelgP9oym1vgYdAABwCg/wiUaJvneZoK2O7e3uZ292rMxecqw8RSx2MGOAAA/USgRNtmKg3N1Rzf6xMj7W93r6adUOm2muKESgAAeo5AibZNBqxOjhVt5e2Q+9UhtBrgYXbQKesAANAfBEq0ze+4IKn9dncYrQZ4JsSv1laopKwDAEDvECjRlkrd0Yn5uu/17V3a7j6bYTQDZdgGeN0lVAIA0CsESrQlqIwzlLU0ko/2rHwa4AAAxA+BEm0JOi5oYiQnowczEWmAAwAQLwRKhNZwPR0pR3tcUFg0wAEAiA8CJUI7Uq7K8Qlltmlocynb0/uhAQ4AQDwQKBFaULt7fCQnswfb3WejAQ4AQP8RKBGK53mB509GcVxQWDTAAQDoLwIlQjkxX1e1sXoKMySND0dzXFA7aIADANAfBEqEErTdvbmUVTbMnnMP0AAHAKD34pECEHtB50/2st0dBg1wAAB6i0CJNc1VGzpdafhen4hoOk4naIADANA7BEqsKaiMM5LPqJSLdjrOetEABwCgNwiUWFPQ85NRze7ulvU0wH26RwAAwAeBEoFqjqujszXf6/08Lqgd7TTAGzTAAQBoC4ESgQ7PVOWXq3IZU2NFu6f30wka4AAARINAiUCTAdvdEyM5GX2YjtMJGuAAAHQfgRK+XM/TdIKOCwqLBjgAAN1FoISvY3M11X2SlGlIW2MwHWe9aIADANA9BEr4mjztvzq5bTinTJglvhijAQ4AQHcQKLEqz/PWOC4omdvdq6EBDgBAZwiUWNVMtaG5muN7fTzm50+2iwY4AADrR6DEqqYCtrs3Fm0V7JBLeglitlHWoQEOAMASAiVWFXRcUNyn43SiVdahAQ4AQHgESpyjUnd0Yr7ue32Qnp9cDQ1wAADaQ6DEOaYCzp4sZi2N5DM9vJv+oAEOAEB4BEqcI7jdnbzpOJ2gAQ4AwNoIlFjBcT0dLvuvUE4M+Hb3amiAAwAQjECJFQ7PVuX4hCHbNLSllO3tDcUEDXAAAPwRKLHC1Gn/7e5tIzmZKdruPlurrGPRAAcAYAUCJc5oTsfx3+4e9HZ3GK2yDg1wAACWEChxxsmFuio+VWVD0vjw4J4/2Q4a4AAArESgxBmTAdNxNpeyyoZZlksRGuAAADSREHBG8HFBbHevhgY4AAAESiyaqzV0utLwvT4xwOMWO0UDHACQdgRKSJKmAra7R/IZlXKDPx2nE+02wKs0wAEAA4RACUnSZMB2N6uT4bTTAJdogAMABgeBEqo7ro7O1nyv8/xkeDTAAQBpRKCEpmeq8tt9zWVMjRXtnt7PIKABDgBIEwIl1tzuNlI8HacTNMABAGlBoEw51/M0XWY6TlRogAMA0oBAmXLH5mqqO6snGNOQtpayPb6jwUMDHAAw6AiUKRd0XNC24ZwyFr9EuoEGOABgkJEWUszzvDWen2S7u5togAMABhWBMsXK1Ybmao7vdc6fjAYNcADAoCFQpthkwHb3xoKtgh0y9aBtNMABAIOEQJliQdvd2zewOhk1GuAAgEFBoEypSt3Rifm673Wen+wNGuAAgEFAoEypoLMni7alDflMD+8m3WiAAwCSjkCZUpOng7e7mY7TWzTAAQBJRqBMIcf1dHi25nud7e7+oQEOAEgiAmUKHZmtyvF5EC9jGtoyxHScfqIBDgBIGgJlCgUdFzQ+nJMZpnaMSNEABwAkCYEyZTzP01TgcUFsd8cFDXAAQFIQKFPm5EJdFZ82h6HmCiXigwY4ACAJCJQpE7TdvXkoq2zY5IKeoQEOAIg70kPKBG13TzAdJ9ZogAMA4opAmSJztYZOVxq+17dzXFDs0QAHAMQRgTJFpmb8t7uHcxmVckzHSQIa4ACAuCFQpkjgdJwRtruThAY4ACBOCJQpUXdcHZ0LmI7DcUGJQwMcABAXBMqUmC5Xfbc9cxlTm4p2b28IXUEDHAAQBwTKlJgK2O6eGM7JMJiOk2Q0wAEA/USgTAHX8zRV9i/ksN09GExDytEABwD0AYEyBY7P1VR3Vk8OpiFtK2V7fEeIikEDHADQBwTKFJgMOC5oaymnjMUvg0FCAxwA0GskiQHneV7wcUFMxxlINMABAL1EoBxw5WpDczXH9/oE03EGFg1wAECvECgHXNB298aCrYIdshqMxKIBDgCIGoFywAUeF8R2d2rQAAcARIlAOcAqdUfH5+u+17ez3Z0qNMABAFEhUA6w6YCzJ4u2pQ35TA/vBnFAAxwAEAUC5QCbnAnY7h5hOk5a0QAHAHQbgXJAOa6nw+Wa73Wm46QbDXAAQDcRKAfUkdmqHJ+9yoxpaMsQ03FAAxwA0B0EygE1FXBc0PhwTlaYZgZSgQY4AKBTBMoBtNZ0HI4LwtlogAMAOkGgHEAnF+qq+DzwZkiaGOb5SZyLBjgAYL0IlAMoaLt701BW2bD1XqQODXAAwHqQLAZQ0HSc7Wx3Yw00wAEA7SJQDpj5mqNTlYbv9Qmm4yAkGuAAgLAIlAMm6DDz4Zyl4RzTcRAeDXAAQBgEygETtN3N6iTWgwY4AGAtBMoBUndcHZnzn46znek4WCca4ACAIATKAXK4XPVdGcpapjYV7d7eEAYKDXAAgB8C5QCZDDguaGIkJ8NgOg46QwMcALAaAuWAcD1P0wGFHLa70U00wAEAyxEoB8TxuZpqzup/YpuGtK2U7fEdYdDRAAcAtBAoB0TQdvfWUk4Zi//V6D4a4AAAiUA5EDzPW+O4IKbjIDo0wAEABMoBUK46mq05vtc5fxJRowEOAOlGoBwAUwFlnI0FW8Ww7QmgAzTAASC9CJQDYJLtbsQIDXAASB8CZcJVG46Oz9d9r3NcEPqBBjgApAuBMuGmAtrdBdvUhnymh3cDLKEBDgDpQaBMuKBAuX0kz3Qc9BUNcABIBwJlgjmup+ly0LhFtrvRf4Yh2RYNcAAYZATKBDsyW5Xjs5yTMQ1tYToOYmQ9DXC2wAEgGQiUCRa03b1tOCcrzMNrQA+12wAnVAJAMhAoE8rzvMDjgrZzXBBiigY4AAweAmVCnVpoqBJwKvQ4z08ixmiAA8BgIVAm1GTAdJzNQ1nlwjYggD6hAQ4Ag4PUkVBB4xbZ7kZS0AAHgMFAoEyg+ZqjUwsN3+sTTMdBwtAAB4BkI1AmUNDq5HDO0nCO6ThIHhrgAJBcBMoEmgw4LojDzJFkNMABIJkIlAlTd1wdmQ0etwgkGQ1wAEgeAmXCHC5Xff/wzFqGNg3Zvb0hIAI0wAEgWQiUCbPWdrdhMB0Hg4EGOAAkB4EyQVzP0zTHBSFlaIADQPwRKBPk+FxdNWf1PylNozm/GxhENMABIN4IlAkSdFzQllJOGYv/nRhcNMABIL5IIAkSNG6R7W6kAQ1wAIgnAmVClCsNzVYd3+ucP4m0oAEOAPFDoEyIoNXJ0UJGxbAPmAEDgAY4AMQLgTIhgre7WZ1EOtEAB4B4IFAmQLXh6vhc3fc6291IMxrgANB/BMoECDp7smCbGi1keng3QPzQAAeA/iJQJgDTcYC10QAHgP4hUMac43qaLvsHSo4LApbQAAeA/iBQxtzR2aocnz/xMqahLSUCJbAcDXAA6D0CZcwFbXdvG87JCrO/B6QQDXAA6B0CZYx5nhc4bpHtbiAYDXAA6A0CZYydWmhooe6/FzfOcUHAmmiAA0D0CJQxFrQ6uXkoq1zYh8SAlKMBDgDRIpHEWNB0nAm2u4G20AAHgOgQKGNqvubo1ELD9zrjFoH20QAHgGgQKGMqaLu7lLM0nGc6DrBeNMABoLsIlDEVdFwQq5NA52iAA0D3EChjqO64OjobPG4RQOdogANAdxAoY+hwuepbBshahjYN2b29IWCAracBTlkHAFYiUMbQVMB298RIXqbBdBygm9ptgBMqAWAlAmXMrDUdh+OCgGjQAAeA9SNQxszxubpqzupLH6bRnN8NIDo0wAGgfQTKmAk6zHxLKSfb4n8ZEDUa4ADQHtJJzARtd29nuxvoGRrgABAegTJGypWGylXH9zrHBQG9RQMcAMIhUMZI0Hb3aCGjYtg9OABdQwMcANZGoIyRtY4LAtAfNMABIBiBMiaqDVfH5mq+1xm3CPQfDXAAWB2BMiamA7a7C7ap0UKmh3cDwA8NcAA4F4EyJibX2O42mI4DxAYNcABYiUAZA47r6XA5KFByXBAQNzTAAWAJgTIGjs7W1PD5k8YyDW0tESiBOKIBDgBNBMoYCDrMfHw4JyvMEgiAvqABDgAEyr7zPC/w/Em2u4FkoAEOIM0IlH12aqGhhbr/cgWBEkgOGuAA0opA2WdB292birZyGabjAElCAxxAGhEo+yxoOs72DRxmDiQRDXAAaUOg7KP5mqOTC3Xf64xbBJKLBjiANCFQ9lHQdncpa2k4x3Y3kGQ0wAGkBYGyj4K2uyc2MB0HGBQ0wAEMOgJlnzQcV0dmA56fpN0NDBQa4AAGGYGyTw6Xq77PS2UtQ5uGsr29IQCRowEOYFARKPtkMmC7e3wkL5PtbmAg0QAHMIgIlH3geV7wcUFsdwMDjQY4gEFDoOyD4/N11XyqnIYhbRsmUAKDjgY4gEFCoOyDydP+xwVtHcrKtvjfAqQFDXAAg4Dk0gdB509OMB0HSB0a4ACSjkDZY+VqQ+Wq43t9O9NxgFSiAQ4gyQiUPTYVsN09ms+oGHaZAsDAoQEOIKkIlD0WdFwQ290AaIADSCICZQ9VG66OzdV8r09wXBAA0QAHkDwEyh6aLvtvd+czpjYW7B7eDYC4owEOICkIlD00dTrgMPMNeRlMxwFwllYDPMx3h4ZLWQdAfxAoe8RxPU2XA56fZLsbgA/TCB8qXRrgAPqAQNkjR+dqavg8OW+ZhraWCJQA/NEABxBnBMoeCTouaNtwVlaYPyUApBoNcABxRaDsAc/zNBkwHYfDzAGERQMcQBwRKHvgdKWhhbr/d3SenwTQLhrgAOKEQNkDQauTm4q2chmm4wBoHw1wAHFBoOyBoOOCmI4DoBPtNsBrDqESQPcRKCO2UHd0cqHue307290AOtROA5yyDoAoECgjNhWw3V3KWhrOZXp4NwAGFQ1wAP1EoIzY5Brb3UzHAdAtNMAB9AuBMkINx9WR2YBxi2x3A4gADXAAvUagjNDh2ZrvlpJtGdo0lO3tDQFIDRrgAHqJQBmhyYDpOBPDOZlsdwOIEA1wAL1CoIyI53maLnNcEID+ogEOoBcIlBE5Pl9XtbH60+6GIY0P8/wkgN6gAQ4gagTKiEwFbHdvGcrKtvhPD6B3aIADiBKpJiKTMwHtbra7AfQJDXAAUSBQRqBcbahcbfhen+C4IAB9RAMcQLcRKCMQNB1nQz6joSzTcQD0Fw1wAN1EoIxA0HQctrsBxAUNcADdQqDsslrD1fG5mu91trsBxAkNcADdQKDssulyVX7fa/MZUxsLdk/vBwDW0mqAhy3r0AAHcDYCZZcFTscZyctgOg6AmLLabIDXea4SwCICZRe5bvB0nO0b2O4GEG/tNMAdjwY4gCYCZRcdnaup4fNwkWUa2loiUAKIPxrgANpFoOyioO3ubaWsrDBVSgCIARrgANpBoOwSz/M0xXQcAAOEBjiAsAiUXXK60tB83fG9Pj7MdjeA5KEBDiAMAmWXBE3HGSvayttWD+8GALqrVdYJgwY4kD4Eyi4JnI4zwnY3gORrp6xDAxxIFwJlFyzUHZ1cqPte57ggAIOCBjiA1RAouyCojDOUtTScy/TwbgAgWjTAAZyNQNkFQccFbWc6DoABRAMcwHIEyg41HFdHZv1XKCfY7gYwoGiAA2ghUHbo8GzN96du2zK0eSjb2xsCgB6jAQ6AQNmhqYDt7vHhnEy2uwGkAA1wIN0IlB3wPE9TZY4LAgCJBjiQZgTKDpyYr6vaWP2BIEPS+AjPTwJIFxrgQDoRKDswGTAdZ0spK9viPy+A9Gk1wDMhvgW2QiVlHSDZSDwdmGI6DgCsyjCagTJsA7zuEiqBJCNQrtNstaGZasP3OscFAQANcCAtCJTrFLTdvSGf0VCW6TgAINEAB9KAQLlOQdvdE2x3A8AKNMCBwUagXIdaw9WxuZrv9e1sdwPAOWiAA4OLQLkO0+Wq/L7H5TOmNhbsnt4PACQFDXBgMBEo1yHo+cmJkZwMpuMAgK/1NMB9jvwFEBMEyja5rqfpmYDnJzfw/CQAhNFOA7xBAxyINQJlm47O1dTweajHMqRtJZ6fBICwaIADg4FA2aapgO3ubcM5WWGeNgcAnEEDHEg+AmUbPM/TJMcFAUDX0QAHko1A2YbTlYbm647v9YkRtrsBYL1ogAPJRaBsQ9B291jRVt4O+XQ5AGBVNMCBZCJQtmEyoN29ne1uAOgaGuBAshAoQ1qoOzo5X/e9PsF0HADoKhrgQHJk+n0D/eB5nlyv+Y3HkyfHbf40bMiQYTS/iZ19OPlUwOrkUNbSSC6V/ykBIFKtUFlz5DuhrKXVAM9aza1zAL0z0CnI8zzVHE91x1Nj8e91x/M9R3K5jGnItgxlLEO2aej4XE2GVv+GNjGSZzoOAESk1QCvu2s3u1tlHTtkYxxAdxieN1gbBJ7nqVJ3NV9zNV93u7r94bieTi3UdGKhplOV+pn3fvVFY9rKgeYAECnPa25thy3h2GZz9wlA9AZmhbLacDVbcboeIpezTEObhnLaNJQ7Ey6Pz9e0eSgbzRcEAJxhGFLGaD5TWQ8RKutuc8UyzDFEADqT+BXKhuPp1EJD87X+nRtRsE2NFi3Z/CgMAD3Rel4yDGvxKCKeTAKik9hA6bqeTlcclSshv6P0QClnakMhw/hFAOgB11s8LijEa83FQ9MJlUA0Ehkoqw1Xx8p1OTG8c9OQNpds5cOeygsAWDfPC1fWkZpb5TTAgWgkLlDOVR2dmGuE+ok0iGU0vwGd/T6W2Z1RXmPFjEp5JucAQNTaDZU0wIHuS0yg9DxPMxVHpxfa3+LO26ayVvMYoNZRQGbAj6iut3TM0GNHZpUxTG3I2zLb/A40nLc0WrA4UggAIkYDHOivxLS8Ty84mmnjecmCbaqYNVWwzbaDoGkYymYM1ZyGDp2YX/yYNJrPaqyY1Ya8Heo5yXLFked5Ghuy2/r6AID20AAH+isRgXK+Fj5MDmVNbShmlOnCfsby6TiuJ51YPIMyYxraOTakjYW1jwuarbrKZhyVcmx/A0DUrMXiTZgGeGNxVCMNcKBzsf/ZrO54Oj7XWPN1uYyh8RFbm0p2V8KkJE3OVFb9eMP1ZJmeJjbYKoQo35yca6gWdh8GANARZoADvRfrQOl5no7N1tf8jT6St7R12Fa2i3sXtYarY7M13+sTI3nZlqnNpYw2FIJXHz1Jx2YbcvmOBQA90QqVYdYXWmda8i0aWL9YB8rZqqv6GmcDjQ1lNFrMdL34Ml2u+jbJcxlTY8Xmc5GGYWhDIaPNpUzgT8MN14vVmZkAMOiMxbMnrRB/PLRmgIdpigM4V2wDped5KleCt7rHipnInk2c8tnulqTtI7lzAmwxa2lTKfiR1HLFYZUSAHrIWJySE2YDqxUqu3F0HJA2sQ2U8zU38PiHoaypoVw0t++6nqaXFXLONjGSX/Xjxayl4YCzJ11PmqvynQoAeqkVKsPOm6i74Y8fAtAU20AZ1OrOmIY2DnV/m7vl6FxNdZ99D8uQtg7nfD93tGDJDthfmak0lJCjPwFgoFhm87nKMBru4lhHvl0DocQyULYOFfczUrACDybvVNB299bhXGCLvPlMpf93LMfVms+FAgCiQQMciEYsA2U1YK/BMprb3VHxPE+TAdvd2322u5cr2GZg6Kw2+O4EAP1CAxzovsQFymIu2lGGM5WG5gNOxJ0Y8d/ubjEMI/D5zqB/PwBA9NptgFdpgAOBYhoo/X/X5jLRjjMIWp0cK9rK2+EewMkFVAoJlADQf+00wCUa4ECQWAZKJ+DHwKCg1g1Bz0/6tbtXkw0Ivo4rijkAEAM0wIHuiGWgDIpaXZqquKpK3dGJ+brv9e0htrtbTMOI9F4BAN1DAxzoTCwDZb9MBWx3D2UtjeSDDy4HACQXDXBg/QiUy0wGbnefOx0niOd5fKMBgIShAQ6sTywDZdBNRXWGY8P1dLjc/nQcP47rv3XPTjgAxBcNcKB9sQyUdmBDOprftUfKVd9vCLZpaEsp29b7BTW5bcuI9OgjAEBnaIAD7YlloAw6GqhSj+Z3bFC7e3wk1/ZknkpAoMyFrRMCAPqGBjgQXiyTTdDRQAt1t+vb3mtNx2l/u9vTfDUgUEZ8liYAoHtogANri2WgzGaMwOcMy5VG177W7Kz05a839N1HLB16LKOF+ZVf2ZA0Phz+uKDm/TmBRx9FfZYmAKC7aIADwWJ5Do5pGCrmTM35rPLNVl0N5dx1B7ODB6V77pG+8AXpqackz7MlbW5eNDxtO8/RS15V1fVvnNdLrjKVbePr1BquyhX/0Y0F25TFAZUAkDitUFl31y7htBrgWau5dQ4MOsOL6ciWuuNq6rT/IeOWKY2PZNsKZ4cOSW97m7R/v5TJSI2AhU7T8uQ6hn7oRxr6s09ktHPn2u/vep6mT9fVCPhOs23EZoUSABLM85pb22Gfvgp7DBGQZLFNNrZlqhDwJLTjSsdn63JD5uF9+6Q9e6T772/+c1CYlCTXaf7u//pXLO3Z0/z8IJ7n6fhsIzBM5jIGYRIAEo4GOHCuWKebDYXgp6ArDU9HZuqBs78l6a67pNtvlyqVtYPk2RzHUKXS/Py77vJ5jevpSLmuhTUa6BsKsXzCAADQJhrgwEqx3fJuOTXf0EzAM4lSc/t7S8le9VnHffuaYbBb9u2Tbr116Z/rjquj5eCVSUkq5SyNDREoAWDQtJ6XDMNaDKI8V4lBE/tA6XmejpbrqoQ40Hwoa2pDMaPM4sMqhw41t7kr/kdMti2fb5Z6XnCBp9MLDc0GHA/UkrUMbRuxOcwcAAaU6y0eFxTitebiJB7+SMAgifWWtyQZhqFNJVtWiDudq7maOlXTqfmG6o6nt72t/S3utTQanm693dXkqVqoMGka0uYSYRIABhkzwJF2sV+hbKk1XB0p10PPS338MUPX/1B74xLb8cWv1nTxC4NvxjCkrSWbyTgAkBI0wJFWiUk62Yyp8ZGsbCvc77x7P2nKsqLJypbl6c8/EfyfrnmsEWESANKEBjjSKlFpJ2MZ2jZsKx9idOH9+005TjQ/9jmOoX/8ov9/umzGWAy/ifrPCwDoAhrgSKPEJR7TNLRl2NZI3v9Iodmy9INnot1DeOZpQ3Oz5368lDO1bdhmGg4ApBwzwJEmiQuUUrOoM1rMaGJDVsXsuf8KzzxtyPOiDXSeZ+jpQ0tfI2sZGh+xNTZEAQcA0MQMcKRFog9GtC1Dm0u2qg1Xp+Ybqi4eLVSr9ubr16rNcyjnKs1W+eRJKZ81VcxaKtimCllLxaypDFvfAJBazABHGiQ6ULbkMqa2jWRVa7iar7kqFHrz413FqevU3NJptp6khZqrhdrKh2FsyzgTLgtZS0XbVM42WckEgJQwFs+eDNMA9yRVHRrgSJbEHBvUjnLZ04YNinTb2zA8fflAWcWh9X6+VLAXVzOzlgqLK5s8ewkAg8vzmoEybAnHNhXqHGag3wZihfJsw8OGdu2Snnwyuq+x4wXuusOk1PymMl9rrqhK9TMfz561mlnIWsplDFYzAWAAGIaUMZrPVNZDhMq621yxDHsMEdAvAxkoJenmm6WPfrT7k3Kk5jmUr/zRCN5YUs3xVFto6PTC0sdMQ0urmHbz7wVWMwEgsazF0YthZoA3Fos6zABHnA3klrfUnLe9d2907//p/bPaubu/B4flMuaZcFnMmirYlrKsZgJAYnitMYwhXssMcMTZwAZKSbrhBun++7u7SpnJSNdd5+lvP9cs3yzUHc3XXC3UHNXDztqKkLV8NbP1fKZtymQ1EwBiyfPCNcCl5lY5DXDE0UAHykOHpD17pEqle++ZzzdXP3fuPPdaw3HPhMuFmqv5uqNKzQ31k2fU8ra54iijQtaSbbGaCQBxwAxwJN1AB0pJ2rdPuv327r7frbeGf73nearUXc23Qubi3xthfhSNmGUaS+WfxcZ5PmvKJGQCQM/RAEeSDXyglKS77pLe977uvM+v/3rn7yNJ9Yar+fqy1cyao0qYyl8P5M8cZ7T0d+aSA0BvOG64BrjULOrQAEccpCJQSs2VxXe+s/k8ZTvPVGYyzb/uvru9lcn1cD1PldYq5mLYnK+5cmKwmplZvpq5GDTzHM4OAJFoTcwJwzJogKP/UhMopeYzlW97m7R/fzMkBgXL1vXrr5c+9rHVn5nsBc/zVHe8M+FyYTFsxmE109DiqMllRxkxahIAuoMGOJIkVYGy5eBB6Z57pPvuax5+vvy/gGFIF10k3XSTdMcd0mWX9e8+g7iut2wVc2nbPAaLmYyaBIAuoQGOpEhloFxudlZ64gmpWpVyOWn3bqlU6vddrY/neao1vBVHGS3UHFUb/f9f3Bo12QyaiyuatqWMxXc+AAhCAxxJkPpAmQaO6604ymhhMWzGYTWTUZMAsDYa4Ig7AmVKeZ6nasM7s4rZWtGsxeBwdkZNAsDqaIAjrgiUWKHhNLfMl87MbJaA4vCrJJcxzqxiMmoSQFrRAEccESixJs/zVK0vnZvJqEkA6C8a4IgbAiXWLc6jJnMZc6llzqhJAAOIBjjihECJrjp71GRrRZNRkwDQfTTAERcESvRE3XFXzDKP66jJ1rmZNk+yA0gIGuCIAwIl+ibuoyaXzzJn1CSAuKMBjn4iUCJWGDUJAOtHAxz9QqBEIqw2anKh5oR+bihKK0ZNLoZNVjMB9AsNcPQDgRKJtXzU5PLnM6thHySKEKMmAfQTDXD0GoESAydZoyZN5TKsZgLoPhrg6CUCJVIhEaMmz2qbM2oSQKdogKNXCJRINcf1zjozk1GTAAYPDXBEjUAJnCXOoyZNQyvPzGTUJICQaIAjSgRKICRGTQJIOhrgiAqBEujA0qhJd8XzmYyaBBBXNMARBQIlEIH4j5o0l57PZNQkkDo0wNFtBEqgR1qjJhfqS89lMmoSQL/QAEc3ESiBPkrCqMmCba14PpNRk8BgoQGObiBQAjG0ctTk0vOZMSiaM2oSGEA0wNEpAiWQEJ7nqeYsmwLEqEkAXUQDHJ0gUAIJl4RRk8ufz2TUJBBfNMCxXgRKYACtOmqy7qjW6P9vd0ZNAvFGAxzrQaAEUiQxoyYXt88ZNQn0T8OlAY7wCJRAyp09arL1fCajJgG02wC3DLbA04pACWBVjeWHs9cXD2eP4ajJ1vOZjJoEokEDHGEQKAGEFvdRk4WsqaK98nB2VjOBztEAx1oIlAA6dvaoydYB7XHAqEmgO2iAIwiBEkAkkjRqstU6Z8scCEYDHH4IlAB65pxRk4vnZsZm1KS9NGKSUZOAPxrgOBuBEkDfMWoSSB4a4FiOQAkglpIzarK1Zc6oSaQPDXC0ECgBJMqZUZNnrWjG4NFMRk0ilWiAQyJQAhgAnuep1mhNAYrvqMnlK5qMmsQgoQEOAiWAgRXnUZPZjLG4ismoSQwGGuDpRqAEkCqe56naWH44e5xHTTaLQBzOjiShAZ5OBEoAEKMmgW6iAZ4+BEoA8NEaNbnUMnc0X3fViMFqpmVqcbucUZOIJxrg6UKgBIA2xX3U5JmWud06nJ3VTPQHDfD0IFACQBesNmpyoeaqEYPzjFYbNZm3TZn8yY0eoAGeDgRKAIjI8lGTCzVX8zEfNVnImrJpSCACg9QAn1uo6vCxGVVrDc1Xanp68rh27diiQs5WLpvR+OYRFfLZft9mzxEoAaDH4j9q0lzcLmfUJLorSQ3w8lxFD33nkB47NK3vP31Ejz9zWI8/fUSTR06t+bnnj2/UJRdu08UXbNUlF27TpTvHde2VO1UsDG7QJFACQAwwahJpEecG+Onygr7wTwf02f2PaP/XvqdqrdG19y7ms7rxVXv1+p94sW561V4NFXJde+84IFACQIwxahKDKE4NcM/z9A8PHtTHP/2A9n/tMdXq3QuRfgp5Wzf+8OW6402v1quuvjjyr9cLBEoASJjYj5q0zzqcnVGTWEUcGuCPfO9ZvffDf6Mvf+Px7r5xG25+9eW6690/rUt3jfftHrqBQAkAA6K1mjnPqEkkRL8a4M9On9Rv3v13+tTnH+78zbrAsky99adfofe/47XaOjbc79tZFwIlAAywuI+aXL6KWcyaytusZqZNrxvgn//yAb31N/5E5bnK+t8kIptHS/rUh25N5DY4gRIAUqixWABqHWXUCptx+AOBUZPpFHUD3PM83f2pf9Sv/dfPqpPoY2cs7dyxWafL8zp8vHzm4+dtG1UhZ+vpyeNqdFCmszOW/r/3v1k/91MvX/d79AOBEgAgKQGjJpcdZdRsmjNqctBE1QBvNBz9yoc+o4/91QNt3c9IKa+bXnW5XnzZ+brkwm265IJtumD7mDIZy/dz6nVHh54/psefPqzHnzmibz76jP7+gUc1X6m19bXfc+sN+uA7XivTTMbZsARKAEAgRk2il7rdAPc8T2//rXv1p3/79VDvOVLK66euu0qv/4kX6cdffqlyWTvczQSYX6jpH756UJ/d/4i+8E8HNLcQLlz+0r/7cf3u//X6jr9+LxAoAQBtc8+sZjJqEt3XzQb4xz/9gN71O3+55vvkshn94luu06/ccr1Gh4vt3XAbjp4o664/vE/7PvMVOc7aP5j9+X+5RW+4/iWR3U+3ECgBAF3BqEl0Uzca4N989Bn92Fs/vObZkm+84SW6890/rQu2j63/htv0+NOH9R9//7O674FHA19XKub04L3v0SUXbuvRna0PgRIAEKnloybPbJ3XHYVYnIkcoybjrZMGeHmuoqt/5nf0g6kTvq83TUMf+pU36B1v/tHOb3YdPM/T7/7h3+u37/l84Ov27p7Qg/e+pyvb71EhUAIAeq61mjkf81GThaypot08P5NRk/2zngb4nfd8Xnd97D7f15WKOf3pf36rbnrV5V26y/X77/d9Q2/74L2BK6m//2s/ozt+9kd6eFftIVACAGIjzqMmbcs4s4rJqMnea6cBXqlUtec179fJmXnf1/zN/3OHbnzV3i7dXef+4vMP65b3/anv9fPHN+rR//Gbsm3/hnk/Zfp9AwAAtFimoVI+o1J+6WNnj5pcqDdXNHs9arLueDq90NDphaWPMWqyd6zF4k2YBvgff/bBwDD5H2+7MVZhUpLe/Jpr9LVvPaWPf+Yrq15/dvqkPv0/v6m3vPaaHt9ZOKxQAgASKRGjJs8UgRg12S1rNcAdx9WVP/UBTR45ter1l11xoe7/xH+QFcNCVqVa17U/+5/1+NOHV71+xcXn6eG/em+P7yocVigBAInkt5oZh1GTtYanWqOhU8s+xqjJ7jCMZvnGrwH+rcd+4BsmJem9t98YyzApSfmcrV996/W6/YN/vur1A99/Xk89e1S7zt/S4ztbG4ESADAwDMNQ3raUty1paKkRG4dRk64nzVUdzVUdSfUzH2fUZPuMxbMnV2uAP/ztQ76fd/nF23XjD8drq/tsb7rpZfqtj35ez02fXPX6Q985RKAEAKAfMpah4UJGw4WlP/ZWGzW5UHd7vppZbTTb7Sfnlxq+jJpcm2FItiUZZzXAHz7wlO/n/OzN18Q+qNu2pTfdeLX+6yf3r3r9698+pDe/Jn7PURIoAQCpZBjG4qqgpTEtrWa2Rk0ufz6zUu/taqbjSrNVR7PVlQ2Us0dNNg9nT/dqZsZsHmzeaoB/48DTvq996RU71XDDzwDvl2uv3Ol77eED/iuw/USgBABgGdsyZRdMjSxbzYzLqMlKvTl56OTc0mpma9Tk8ucz0zZqstUAr9RdPX949a3ijGXqRZe+QA1X8kLMAO+na6/yD5TPTPof1N5PBEoAANZgGs3WdjFradPix1aMmqwvHc7e61GTDddTueKoXFlazUzjqEnTkLIB/3qlobzyueZKtONJnhs8A7yftmws+V5z43Ao6yoIlAAArINhGMpmDGUzpjYs+3gcRk16UvMe6q5OzC19fNBHTbbzr+EuHj9kLxvXiPUjUAIA0EWmaWgoZ2kotzTRJC6jJuuOp/qCo5mFlauZhezgj5p0HFee560Iz56aoTIbs1BZb4Q4vT1mCJQAAERs+WrmaHHp447rqVJffjh770dNepLma82zO5dL4qhJ0zRUzGc1X6mdc608V9GzUyf0gu2bzrlWc5ZmgMfBgcef9702vPzg1RghUAIA0CeWaWgol9FQbuljcR81mbdbATN+oyYNw9CVL9yhr3979aODHj5waNVAKTVb4p7i0QB/6Dv+Te4XvXBHD+8kPAIlAAAxYhiGcrahnG1q47LD2c8eNblQc7VQ7+1qpustX81cOpw9TqMmr73yQt9A+Y8PPaY3/purfT83Lg3wL339Md9r1wQcKdRPBEoAABIgaNTk8ucy4zpqshU2o17NfPlVu/QHf/alVa99+u//Wb9xx2s1sWXU9/P73QD/10PT+sIDj/pef8VVu3p4N+ERKAEASKjloyY3JmjU5PLnMrs9avKHXnyRLMuUs0qtvlZv6J6/+Ef91rt+es17X6sBPjsrPfGEVK1KuZy0e7dU8j/tJ7Tf/5MvyvNW/z9VzGf10r0v6PyLRIBACQDAgAk1anIxbPZr1OSpiEZNbhkb1ut+7Cp9dv8jq17/w7/6sv63G6/W3kuCn0VcrQF+8KB0zz3SF74gPfWUtDz3GYa0a5d0883S298u7dnT9q3rK998Qvd+7mHf6z/3k9eqkM+2/8Y9YHh+MRgAAAy8OIya9HNm1KS9FDTDrGY+8r1n9cq3/J7v9V3nb9aX/+w9Kg0VfV+z3HPPSO+4Q9q/X8pkpEbD/7Wt69dfL33sY9LOkI88Th09rVe+5fc0fWxm1euWZeq7f/sBXXje5nBv2GMESgAAsMLyUZPLn8/s9ajJ1YQdNfmT77hbX/yaf7nlxh/eqz//L7fJsm3f10jSJ/5I+g+/1AyJQUHynPvMNP/6yEek224Lfu3M7IJe/6579NVHnvR9zZtuulqf/J2fD38DPUagBAAAa/I8Tw3Ha65i1p2+jZpczWqjJr/zvad1/a3/LfDzrr1yp/777/+CRjcMr3r9935X+s0PdH5/d94p/cZvrH7tB1Mn9IZ336Pvfn/S9/Mty9TDf/le7bloovObiQiBEgAArFscRk36+eif/r3u/Zt/CnzNBds36VMfuk17Lzl/xTb/J/5Iesfbu3cv+/ZJt9668mMPPvKk/vdf3afDx8uBn3vnu1+nX/7567t3MxEgUAIAgK46e9Rk6/nMXo+abDiO3v2BP9K3Dz4d+DrDMPRzP3mtfv3tr9X4llE9fUh68ZVSpdK9e8nnm6WenTulQ88d0/s/8j/01//wL2t+3mt/5Ar91Yd/IdYTiiQCJQAA6JF+jJo8dmJGt/zy3TpxanbN1+Zztm79mVfrwb+5Wf/8UK6tZybXksl4uvqaql75us/p459+INS87p07Nuurn3qPRofDlYf6iUAJAAD6ZsWoyXoraHZ31ORjTzynX/5Pn9Tp8vyar3Ur46o99r6ufe2zZS/9bZn5w2u+bvvWUX3hnl/UC3eOR3Yv3USgBAAAsdPtUZPPTR3Xr975J3p28ljg6+rPvVHOsVdJstb3hQI5sjY/IHvHZwJf9aJLd+gzf/B2nbd1NIJ7iAaBEgAAJEKnoyZnZhf0vt+7V//y3dVnfUtS5eAHpdqWbt3yubJHlN/zn3wvv/ZHr9Qnf+ffa6iQi+4eIkCgBAAAidZwvMXJPytXNFcLOPV6Q/v+4ov6y797UI2znmP0nJyqB/5vNQ8iioqn3BW/LMOqrfhoIW/r1279N/qVt94gyzIj/PrRIFACAICBs9aoyeenjutj9/6DvvTggTOf486fp9rj74383rKX/K7M4vOSmg3z/+OnrtX773iNdmzbGPnXjgqzvAEAwMAxDKM5EzxraUxL03BaoyZ3bMzpig/+Oz307af03z5xX/NoIa9HsWjx61z/yst057tepytfGDxXPAlYoQQAAKnmep6+9a/P6//94yf0xx/60ci/3i994Kv6P2/Zrd0XbI38a/UKgRIAAEDS7Kw0MiJFmYwMQ5qZkUql6L5GPyTvqU8AAIAIlErSrl3Rfo2LLhq8MCkRKAEAAM64+WYpE9GjlJmMdNNN0bx3v7HlDQAAsOjgQWnv3mjf/7LLonv/fmGFEgAAYNGePdL113d/lTKTab7vIIZJiRVKAACAFQ4dagbLSqV775nPN1cnd+7s3nvGCSuUAAAAy+zcKX3kI919z7vvHtwwKREoAQAAznHbbdKdd3bnve66S7r11u68V1yx5Q0AAOBj3z7pne+UGo3mX2FlMs2/7r578MOkxAolAACAr9tuaz77eN11zX9eq6zTun7ddc3PS0OYlFihBAAACOXgQemee6T77pOefHLlRB3DaB5aftNN0h13DG6b2w+BEgAAoE2zs9ITT0jVqpTLSbt3D+YEnLAIlAAAAOgIz1ACAACgIwRKAAAAdIRACQAAgI4QKAEAANARAiUAAAA6QqAEAABARwiUAAAA6AiBEgAAAB0hUAIAAKAjBEoAAAB0hEAJAACAjhAoAQAA0BECJQAAADpCoAQAAEBHCJQAAADoCIESAAAAHSFQAgAAoCMESgAAAHSEQAkAAICOECgBAADQEQIlAAAAOkKgBAAAQEcIlAAAAOgIgRIAAAAdIVACAACgIwRKAAAAdIRACQAAgI4QKAEAANARAiUAAAA6QqAEAABARwiUAAAA6AiBEgAAAB0hUAIAAKAjBEoAAAB0hEAJAACAjhAoAQAA0JH/HxvQZL9rHBKrAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "graph_point_nx = nx.from_numpy_array(graph_point)\n", "edges, weights = zip(*nx.get_edge_attributes(graph_point_nx, \"weight\").items())\n", "pos = nx.spring_layout(graph_point_nx)\n", "nx.draw(\n", " graph_point_nx,\n", " pos,\n", " node_color=\"b\",\n", " edgelist=edges,\n", " edge_color=weights,\n", " width=5.0,\n", " edge_cmap=plt.cm.Blues,\n", ")" ] }, { "cell_type": "markdown", "id": "d565bc0b", "metadata": {}, "source": [ "### Set of Graphs: GraphSpace\n", "Graphs can have different numbers of nodes and different node labels or order. We assume the existence across the populations of at most $n$ distinct nodes and we add fictionally null nodes to smaller networks, so that all graphs can be described by a fixed-size adjacency matrix. Graph Space is initalized by the maximal number of nodes in the set. " ] }, { "cell_type": "code", "execution_count": 10, "id": "050724e5", "metadata": {}, "outputs": [], "source": [ "total_space = GraphSpace(n_nodes=4)" ] }, { "cell_type": "markdown", "id": "f442a08f", "metadata": {}, "source": [ "Within GraphSpace, we can sample points from random adjacency matrices, we can check if the points belongs and we can return a set of adjacency matrices." ] }, { "cell_type": "code", "execution_count": 11, "id": "5748f785", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ True, True])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "points = total_space.random_point(2)\n", "\n", "total_space.belongs(points)" ] }, { "cell_type": "markdown", "id": "a7cba4c8", "metadata": {}, "source": [ "### Labelled or Unlabelled?\n", "Graphs can be considered labelled or unlabelled, meaning that the observation in the set can share the same nodes (labelled) or not (unlabelled). We can describe both cases into Graph Space by using a little trick.\n", "\n", "To deal with unlabelled nodes, alignment two graphs corresponds to finding optimal permutations of their nodes. By allowing node permutation, a concept of equivalence class is introduced (i.e., every network is associated to the set of \n", "networks obtained by permuting the original one). In geometrical terms, this can be described as a quotient space obtained by applying a permutation \n", "action to the set of adjacency matrix. \n", "\n", "In details, the group of node permutations $T$ can be represented via permutation matrices, acting on $X$ \n", "through matrix multiplication. The binary operation: \n", "\n", "$$\\cdot \\colon T \\times X \\rightarrow X, (T,x)\\mapsto Tx $$\n", "\n", "thus defines an action of the group $T$ on $X$. The obtained quotient space $X/T$ is called graph space, and \n", "each element of $X/T$ is an unlabelled graph $G$, represented as an equivalence class $[x]=Tx$ which contains all the flattened adjacency matrices \n", "in $X$ which can be obtained from $x$ by permuting nodes. The map $\\pi \\colon X \\to X/T$ given by $\\pi(x) = [x]$ can be thought of as a \n", "projection of the Euclidean total space $X$ onto the graph space $X/T$, and the total space $X$ plays a similar role relative to graph space, \n", "as the tangent space does for manifolds, by providing a Euclidean space in which approximate computations can be carried out and projected back onto \n", "the space of interest -- in our case the graph space $X/T$.\n", "\n", "\n", "To deal with labelled nodes, we restrict the set of permutation matrices to the identity: $T=\\{I\\}$" ] }, { "cell_type": "markdown", "id": "ac7b3838", "metadata": {}, "source": [ "### Graph Space Metric\n", "\n", "To define a metric on graph space, we need to choose a metric on the total space.\n", "Any metric $d_X$ on $X$ defines a quotient pseudo-metric\n", "\n", "$$d_{X/T}([x_1],[x_2])=\\min_{t\\in T}d_X(x_1,t^Tx_2t)$$\n", "\n", "on $X/T$. Since the permutation group $T$ is finite, $d_{X/T}$ is a metric, and the graph space $X/T$ is a geodesic space. In the implementation, we suppose that the default metric in the total space is the Frobenius metric between adjacency matrix." ] }, { "cell_type": "markdown", "id": "b420d7ca", "metadata": {}, "source": [ "In `geomstats`, we can equip the total space with a group action, and then get the corresponding quotient structure:" ] }, { "cell_type": "code", "execution_count": 12, "id": "22aa982f", "metadata": {}, "outputs": [], "source": [ "total_space.equip_with_group_action() # permutations by default\n", "total_space.equip_with_quotient_structure();" ] }, { "cell_type": "markdown", "id": "a7c79820", "metadata": {}, "source": [ "The graph space becomes then available under `total_space.quotient` and is equipped with a quotient metric." ] }, { "cell_type": "code", "execution_count": 13, "id": "b9d9554c", "metadata": {}, "outputs": [], "source": [ "graph_space = total_space.quotient" ] }, { "cell_type": "markdown", "id": "f1af545d", "metadata": {}, "source": [ "Let's compare the difference between distances in the total and graph spaces:" ] }, { "cell_type": "code", "execution_count": 14, "id": "d0d87f62", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dist quotient <= dist total:\n", "1.3431542528099538 <= 1.4460818164862586\n" ] } ], "source": [ "point_a, point_b = points\n", "\n", "dist_total = total_space.metric.dist(point_a, point_b)\n", "dist_quotient = graph_space.metric.dist(point_a, point_b)\n", "\n", "print(f\"dist quotient <= dist total:\\n{dist_quotient} <= {dist_total}\")" ] }, { "cell_type": "markdown", "id": "bea5cc08", "metadata": {}, "source": [ "### Graph to Graph Alignment\n", "\n", "The metric on Graph Space relies on the optimization along the quotient fibers. In this context the optimization problem is known as graph matching (or network alignment) and corresponds in finding a match between the two sets of nodes which minimizes the distance between the corresponding adjacency matrices. The distance function needs an aligner input, which solves the minimization problem by returning the second input graph optimally aligned. One of the available aligners are:\n", "\n", "1. 'FAQ': the Fast Quadratic Assignment Matching implemented in `scipy.optimize.quadratic` assignment which is the state of the art in the matching literature based on the Frobenius norm." ] }, { "cell_type": "markdown", "id": "3f87b520", "metadata": {}, "source": [ "The aligner algorithm can be set in the object `total_space.aligner`, which connects the total and the quotient spaces (it has a similar role as `FiberBundle`, but with less structure)." ] }, { "cell_type": "markdown", "id": "94abfdb2", "metadata": {}, "source": [ "We can align a set of points using the following function, which returns the permuted graphs:" ] }, { "cell_type": "code", "execution_count": 15, "id": "26f349ca", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "array([[[-0.34330033, -0.35911276, 0.25708028, 0.23632492],\n", " [-0.14433691, -0.15890698, 0.16680305, -0.28289936],\n", " [ 0.06142698, -0.37582122, -0.18026352, 0.45321387],\n", " [-0.36264321, 0.0694131 , 0.47566548, 0.00336706]],\n", "\n", " [[ 0.16766421, -0.4658085 , -0.04388063, -0.34414864],\n", " [-0.02395103, -0.33029756, 0.39625834, -0.12660624],\n", " [-0.12030707, 0.35831659, 0.14606105, 0.0834617 ],\n", " [ 0.16835003, -0.32220738, 0.34924802, -0.05762742]]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aligned_points = total_space.aligner.align(point=points, base_point=points)\n", "\n", "aligned_points" ] }, { "cell_type": "markdown", "id": "880735ce", "metadata": {}, "source": [ "The permutations of the nodes computed by the `align` function are saved in the `perm_` attribute. Notice that only the last run output is saved." ] }, { "cell_type": "code", "execution_count": 16, "id": "adc454c6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 2, 3],\n", " [0, 1, 2, 3]])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "total_space.aligner.perm_" ] }, { "cell_type": "markdown", "id": "e8b43262", "metadata": {}, "source": [ "### Graph to Geodesic Alignment\n", "In different algorithms for graphs, we need to align a graph to a geodesic. Given a point $[X] \\in X/T$, a $\\gamma : \\mathbb{R} \\rightarrow X$, a domain $[s_{min}, s_{max}] \\in \\mathbb{R}$, the alignment with respect to a geodesic is performed as follow:\n", "\n", "1. Sample $s_i\\in [s_{min}, s_{max}]$\n", "2. Compute $\\gamma(s_i)$\n", "3. Compute $d_i = min_{t\\in T}d_X(\\gamma(s_i), t^TXt)$ is minimum\n", "4. Select the $t^TXt$ corresponding to the $d_i$ minimum \n", "\n", "The algorithm is described in: Huckemann, S., Hotz, T., & Munk, A. (2010). Intrinsic shape analysis: Geodesic PCA for Riemannian manifolds modulo isometric Lie group actions. Statistica Sinica, 1-58. " ] }, { "cell_type": "markdown", "id": "5f309ece", "metadata": {}, "source": [ "To perform the alignment between the geodesic and the point, we need to to define which methodology to adopt. This is specified in the `set_point_to_geodesic` function." ] }, { "cell_type": "code", "execution_count": 17, "id": "fef23c78", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "total_space.aligner.set_point_to_geodesic_aligner(\"default\", s_min=-1.0, s_max=1.0)" ] }, { "cell_type": "code", "execution_count": 18, "id": "9be2c5c8", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.09464325, -0.22793224, 0.18702875, -0.35236845],\n", " [-0.13903918, -0.25412397, 0.29427821, -0.13448654],\n", " [ 0.13459113, 0.03203312, 0.00102791, 0.07366849],\n", " [ 0.19856109, -0.33860977, 0.30828458, -0.18459316]])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geodesic_func = graph_space.metric.geodesic(points[0], points[1])\n", "\n", "total_space.aligner.align_point_to_geodesic(geodesic=geodesic_func, point=points[1])" ] } ], "metadata": { "backends": [ "numpy" ], "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 }