{ "cells": [ { "cell_type": "markdown", "id": "f2d7b3c9", "metadata": {}, "source": [ "# Steiner Trees (STs)\n", "\n", "A Steiner tree is a concept stemming from combinatorial optimization that seeks to connect a given set of vertices (or nodes) in the shortest possible total distance, but with the added provision that additional nodes not in the original set can be included if they reduce the overall connecting distance. \n", "\n", "Given a weighted graph $ G = (V, E) $ and a subset $ T \\subseteq V $ of terminal vertices, the goal of the Steiner Tree problem is to find a tree of minimum weight in $ G $ that spans all vertices in $ T $. This tree can include vertices not in $ T $, and these additional vertices are called Steiner vertices. The objective is to connect all the terminals using the shortest possible total edge weight.\n", "\n", "One notable application of Steiner trees is in the analysis and modeling of biological networks. By using Steiner trees, researchers can identify minimal subnetworks that connect a given set of proteins or genes, potentially pinpointing previously unobserved or unknown elements (analogous to the Steiner points) that play crucial roles in these networks. Such insights can help in understanding disease mechanisms, drug targets, and the fundamental underpinnings of biological processes.\n", "\n", "The Steiner tree problem, like many other combinatorial challenges, is an NP-hard optimization problem. This implies that finding optimal solutions for large instances of the problem can be computationally hard. To tackle this, numerous heuristic methods have been developed which offer reasonably good solutions within an acceptable timeframe. However, many of these methods are typically tailored for specific cases, for example undirected graphs. Certain assumptions inherent in these methods can make them difficult to modify or adapt for different problem scenarios.\n", "\n", "CORNETO stands as a robust tool in this realm, offering flexible modeling capabilities for Steiner trees. By modelling steiner tree problems as Integer Linear Programming problems, CORNETO offers a more advanced and flexible approach to modeling Steiner trees. This allows users to model and solve Steiner tree problems in a variety of ways, including:\n", "\n", "- **Directed/undirected edge support**: Unlike many tools that restrict you to either directed or undirected edges, CORNETO supports both simultaneously. This is particularly useful when modeling systems where interactions can be both reciprocal (undirected) and one-way (directed), which is often the case in biological networks.\n", "\n", "- **Solution flexibility**: With CORNETO, users can decide how strict they want to be with their Steiner tree solutions.\n", " - For cases where only Steiner tree solutions are acceptable, CORNETO can enforce infeasibility for non-Steiner trees. This means every solution provided, even if not the best possible, will always be a valid Steiner tree.\n", " - Alternatively, if users prefer to explore a broader range of solutions before arriving at an optimal Steiner tree, CORNETO permits non-Steiner tree solutions. Only the optimal solution, once identified, is guaranteed to be a true Steiner tree.\n", "\n", "- **Problem extensions**: By transparently modeling the Steiner tree problem as an Integer Linear Programming task in CORNETO, users gain the flexibility to introduce additional constraints or modify the objective function as they wish. For instance, they can set a limit on the number of Steiner points included in the solution or define a maximum distance between any two nodes within that solution. This adaptability enables users to customize the problem according to their unique requirements.\n", "\n", "In this tutorial, we will demonstrate how to model and solve undirected Steiner tree problems using CORNETO. We will use a small example network to illustrate the process, and we will compare the solutions with those obtained with the very well known `networkx` graph library." ] }, { "cell_type": "markdown", "id": "0488a868", "metadata": {}, "source": [ "## Heuristic steiner tree from NetworkX\n", "\n", "First, we will create a random undirected graph with random weights between 1 an 10. Then we will run the heuristic steiner tree algorithm from NetworkX to see which solution it finds." ] }, { "cell_type": "code", "execution_count": 1, "id": "53120fb1", "metadata": {}, "outputs": [], "source": [ "import networkx as nx\n", "import numpy as np\n", "from networkx.algorithms.approximation import steiner_tree\n", "\n", "# Generate a random graph with networkx\n", "\n", "\n", "def create_random_graph(n_nodes, prob, seed=0):\n", " np.random.seed(seed)\n", " random_graph = nx.erdos_renyi_graph(n_nodes, prob, seed=seed)\n", " weights = np.random.randint(1, 10, len(random_graph.edges))\n", " edges_r = list(random_graph.edges())\n", " for i, weight in enumerate(weights):\n", " edge = edges_r[i]\n", " random_graph[edge[0]][edge[1]][\"weight\"] = weight\n", " return random_graph\n", "\n", "\n", "random_graph = create_random_graph(50, 0.1)" ] }, { "cell_type": "code", "execution_count": 2, "id": "7e3f6074", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cost: 45\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABiUElEQVR4nO3dd3hUZd4+8PucacmkF0IgZAgQAhqKigiilKyIq2AEFCYquK6yvrq67vLa+LGva1tRV1fcFcsqrAXRDCII4krT0AQCNtoKJJRMaGmTMiWZZOac3x8hY4b0zCRTcn+uK5fknDPPeSaJ5OY55/s9gizLMoiIiIiIOkn09QSIiIiIKLAxUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyiNLXEyAiIiIKNFa7A6fKrKh1SFArRaTEhSFM03NjVc9950REREQdkFdkxopcI3KOFsNoskFutE8AoIvVImNIAu4co8Pg3hG+mqZPCLIsy20fRkRERNQzFZpsWLjmIHbkl0IhCnBKLUenhv3jU+OxaMZwJMdqu3GmvsNASURERNSC7H1GPLXuMByS3GqQvJhCFKAUBTyTmY6s0bounKF/YKAkIiIiasaSnDy8sumYx+M8OiUND2UM9sKM/BervImIiIgukr3P6JUwCQCvbDoGwz6jV8byV1yhJCIiImqk0GTD5MXbYHdIkOw2VO7KRm3RSdQWHYdUXYWoa25H9Pg7W3y9LMsoWrEA9tOHEXHFVMROeQAapYgt8ycG7T2VXKEkIiIiamThmoNwXLhfUqo2w/zTRsjOOmjTxrbr9eYf1sNRcc5tm0OSsXDNQa/P1V8wUBIRERFdkFdkxo78UlcBjiIqAcl/ykbinS8ieuJv2ny9o6IIFds+ROz197ttd0oyduSXIr/Y3CXz9jUGSiIiIqILVuQaoRAF1+eCIEAQhFZe4a5swxKEpFwG7ZBxTfYpRAEf7QnOeykZKImIiIguyDla3KH2QI2Z92+E/dyxJquTDZySjJxjxZ5Mz28xUBIREREBsNgdMJpsnXqtw1yK8m/+jZiM30IZEdficcYyG6x2R2en6LcYKImIiIgAFJRZ0dnWN6YNb0CdMADhI29o9TgZwKkyayfP4r8YKImIiIgA1DqkTr3OemQnqk/+gJiM30K2WyHVWCDVWAAAstMBqcYC2fnLqmRnz+PPlL6eABEREZE/UCs7t85WV1IASE6c//CRJvss+zfCsn8jes38M7RpV3t0Hn/GQElEREQEICUuDALQ4cve4cMnI0Q3vMn2ok8WInTwWERemQlVr/4AAOHCeYINAyURERERAEeNFVHKOlQ4VG7bq49/B6muBnJtNQCgrqwQ1iM7AQChg66EMro3lNG9mx1TGRGHkP4jXJ/r4rQI0wRf/Aq+d0RERETUTrIsY9euXXj33XexcuVKhF77G0RccRMg/HJZumzjm3BW/dLux3ZkJ2wXAmXS/csgRoe061wKUUBGWoJ334Cf4LO8iYiIqMcpKSnBhx9+iKVLl+LIkSMYMGAA5s2bh4mZety54kiXnXfL/AlITYjosvF9hSuURERE1CNIkoQtW7Zg6dKl+PzzzyEIAmbOnIklS5YgIyMDoli/Kjk+tRS7TpR1usF5cxSigHED44IyTAJcoSQiIqIgd/r0abz33ntYtmwZCgoKkJ6ejnnz5mHu3LmIi2vahLzQZMPkxdtg92J7H41SxJb5E5Ecq/XamP6EgZKIiIiCTl1dHb788ku8++672LBhA0JDQ5GVlYV58+ZhzJgxbT6fO3ufEQtWH/TafF6aORz60TqvjedveMmbiIiIgkZ+fj6WLl2K999/H0VFRbjqqqvw9ttvQ6/XIzIyst3jZI3WodRixyubjnk8p8emDAnqMAlwhZKIiIgCXE1NDT777DMsXboUW7duRXR0NObOnYt58+ZhxIgRbQ/Qiux9Rjy17jAcktyheyoVogClKODZzPSgD5MAAyURERF5wGp34FSZFbUOCWqliJS4sG7rs3jgwAEsXboUH330EcrLyzFp0iTMmzcPM2fORGhoqNfOU2iyYeGag9iRXwqFKLQaLBv2j0+Nx6IZw4P2nsmLMVASERFRh+QVmbEi14ico8UwmmxuT5YRAOhitcgYkoA7x+gwuLd3q5rNZjOys7OxdOlS7N27F71798Zvf/tb3HPPPRg8eLBXz3Ux1/s+VgxjWTPvO06LjLQEzBmrC9pq7pYwUBIREVG7+GqlTpZl5ObmYunSpcjOzkZ1dTVuvPFGzJs3D1OnToVKpWp7EC/z5cqsP2KgJCIiojZ5ei/hM5npyOrgvYRlZWVYvnw5li5disOHD6N///649957cffddyM5Obmjb4G6EAMlERERtWpJTp5Xqp0fnZKGhzJavywtSRJycnKwdOlSrF69GrIsY/r06Zg3bx4mT57saj5O/oWBkoiIiFrUXf0Yz549i/fffx/Lli3DiRMnMHToUMybNw933XUXevXq5bXzU9dgoCQiIqJmdfUTYxwOB7766issXboUX375JdRqNWbPno3f/e53GDduXJvNx8l/MFASERFRs+Yuy231mda154+j4ttPUHvuGKQaKxSRvRB26UREjpkBURXS7GsUooDL+mgx9NxmvPfeezh37hxGjRqFefPm4fbbb0dUVFRXviXqIgyURERE1ERekRnXv7a9xf21pUacf/9PUMYmIerq2RBDI2EvPITKXQaEDhqNhNuebHX8qk8ew+1TMzBv3jxcfvnl3p4+dbOeW99ORERELVqRa2y1NZD1v9sgO2rRa8ZCqGL6AABCU0bCaS2H5acNcNZYoAgJb/a1ImQ88OrHeH7mZV01fepmLJUiIiKiJnKOFrfaHkgQFQAAUePeX1LUhAGCCEFsec1KgoCdJ8q9M1HyCwyURERE5MZid8BosrV6TPiw6yBqwmDa+CbqKs5Dsttgy98L808bEHHFVIjq5u+hbGAss8Fqd3hz2uRDvORNREREbgrKrGirwEIZ3RuJd72C4tXP4+zb81zbI0bdjJjJ97V5DhnAqTIr0vuyCCcYMFASERGRm9p2tAlyVBSheNWzUITFIHr6/4NCGwX7uaOo3GWAVFeD+Jv+6JXzUGBgoCQiIiI3amXbd8SVb3sfkr0afX77uuvydohuGBShkSj7zz8QPuxXCNEN9/g8FBj4nSQiIiI3KXFhaKuleG3RSajik5vcK6nuU/9oxbqSglZfL1w4DwUHBkoiIiJyE6ZRIjlW2+oxivBY1JUaIdVWu223nzlSvz8ivtXX6+K0CNPwQmmwYKAkIiIiF5vNhmXLlqHox68hS84Wj4scfQskWxWKsv8P1p93oPrUflTuWonyb5ZBFa9D6KBRLb5WIQrISEvoiumTj/BJOURERITjx4/jrbfewr///W9UVFTguplzkDdY3+pragoOoHLPp6grPgXJboMiMh7a1KsQefUsKEIjW33tlvkTkJoQ4c23QD7EtWYiIgBWuwOnyqyodUhQK0WkxIXxchwFPUmS8NVXX+GNN97Ahg0bEBMTg3nz5uGBBx7AgAED2nyWd0j/EQjpP6JD51SIAsYNjGOYDDJcoSSiHiuvyIwVuUbkHC2G0WRz67snANDFapExJAF3jtFhcG/+8qPgYTKZ8O9//xtvvfUWTpw4gVGjRuHBBx9EVlYWQkNDXccVmmyYvHgb7F5s76NRitgyf2Kb92hSYGGgJKIep9Bkw8I1B7Ejv7TVZxUDcO0fnxqPRTOG85cgBbQffvgBb7zxBj7++GNIkgS9Xo8HH3wQV111FQSh+bru7H1GLFh90GtzeGnmcOhH67w2HvkHBkoi6lGy9xnx1LrDcEhyq0HyYgpRgFIU8ExmOrL4y5ACiN1ux6pVq7BkyRLs2bMHOp0O999/P+69914kJLSvMGZJTh5e2XTM47k8NmUIHsxI9Xgc8j8MlETUY3jrl+KjU9LwUMZgL8yIqOsUFhbi7bffxrvvvouSkhJMnjwZDz74IKZNmwalsuP3B3v6j7FnM9O5MhnEGCiJqEfgZTvyle4s+JJlGd988w3eeOMNrF27FuHh4fjNb36D3//+9xg6dKjH4/N2EWoJAyURBb2WCgvM+zfC9NXrEFQh0D2yyrW9pvAwrAe/Rm3RcdSWFgBOB5LuXwZldG/XMSwsoNZ0d8FXVVUVPvjgA7z55ps4cuQIhg0bhgcffBBz5sxBeHi4x+NfzPX+jhXDWNbM+4vTIiMtAXPG6ljN3UMwUBJR0Guu9YnDXIqzSx+EqNJAstvcAmXFzo9hObAF6t4DIdmtsBsPNgmUDa1Plt87plvfC/m37l7BO3z4MN544w0sX74c1dXVmDlzJh588EFMmDChxSIbb2PLLQIYKIkoyOUVmXH9a9ubbC/+9BlAECCGRMB29Fu3QCnLEgSh/kFilbmrUZHz7yaBsgGbM1OD7ir4qqurw9q1a7FkyRJs27YNiYmJuO+++3DfffchKSnJk7dA1Gl89CIRBbUVuUYoRPeVGsuhHNQUHkLslN83+5qGMNkWhSjgoz1Gj+dIgW9JTh4WrD4Iu0PqUJgEAKckw+6QsGD1QSzJyWvxuPPnz+PZZ59FSkoKZs2aBafTiezsbBQUFOCZZ55hmCSf4po0EQW1nKPFbr/gndYKlH/9LmIm3Q1lZLxHYzslGTnHivE00j2dJgWw7H1Gr3QPAIBXNh1Dr3CNq+BLlmV8++23eOONN/DZZ59BpVJhzpw5+P3vf4+RI0d65ZzewMvexO82EQUti90Bo8nmts206U2oYpMQfvlNXjmHscwGq93BX549VKHJhqfWHW6yvaWCr4IXp7U4ljK2H5Luext/WXcYl/XRYvtXa7BkyRIcOHAAgwcPxt/+9jfcfffdiI6O7oq30mF80hQ1xr8BiShoFZRZ3X7JWY98C1v+XvT57T+9VrAgAzhVZkV63yivjEeBZeGag3BcdInbYS5F+Tf/hiI8FpLd/R80iXNfaTKG/exRlH/9LrRpYwEAtXUO/GrBuzj38Z8xbdo0vPzyy5g8eTJE0T/uUmtP4ZEMoMBkw/LcAry/+xRbB/UADJREFLRqG7UJkmqrYdr8FiJH3QxleCykGgsAQJYc9ftrLICohKgO8eg81HPkFZmxI7+0yXbThjcQkpzuKvhqTJPUtBek+acNAASEj5wCAJAFEarkEcj5/mdMvHxIl8y9sxoXHgFo837Rhv27TpRh8uJtfNJUEGOgJKKgpVb+sqIj2aogWStQtXcNqvauaXJs4WtZCB08Fgm3/p9H56Geo6Hgq3Goaij46jvvLVRsX97mGJLdBtuRndDohkEV09e1XSEKyCl0YOLlXTL1TvHkSVPOC5XvC1YfRKnFzidNBSEGSiIKWilxYRBQf/lNER6D3rcvanJM5Z5VsBceQsKspyFqIzt8DuHCeajn8UbBl/Xn7ZDralyrk66x/KzgqysLjyg4MFASUdAK0yihi9WiwGSDoFQjpP+IJsdYDn4NCKLbPqetEjXG+sc01pWcAgBUn/geojYSCm0UQnTDXcdGirUwFZ9DWHJy174Z8iveKviyHNgMUROGsCHXNNnnLwVfLRUeeeIv6w5j3KB43lMZRBgoiSioZQxJwPLcgg71BqwrKUDp5y+6bTNtehMAoEkehsQ7L+yTJZz9fjN0upkYN24c9Ho9Zs2ahT59+nht/uSfvFHwVVtSgNqzRxFxxVQISnWT/f5S8NVc4RHQciU7AMhOB8zffwHLwS1wlJ8DFCqo45MRnXEvQvpdAockY+Gag3zSVBBhoCSioHbnGB3e332qxf3x0+YD0+a7bQvpPwL9F6xve3BBxOY3/4JDu8bBYDDg0UcfxZ/+9CdMmDABer0et956KxISEjx8B53H3oDeZzabcfz4cXyz/wQADYDOF3xZDmwCAISPvKHF8z2+YCFSY5To06dPk4/Y2Ngur/xuqfCotUp2WXKiZPXzqDn9X0SNuRWapKGQ6uyoPZ8Pua4GQP0l/R35pcgvNvNJU0GCj14koqDX3LO8PdXcs7zLy8vx+eefw2AwYMuWLQCAjIwM6PV6zJw5E7GxsV47f0vYG9AzsiyjuLgYx48fb/KRn5+PkpISAIAqYQD63vM6AMBRUYQzb9/b6rgXF3zJzjqcXvIbKKMS0Ofu11p8Xb9Dy2E6fgDnz59HeXm52z6VSoXevXu7hczExMQmwbN3795QqVSd+no8ve5wsyv8rT26tGrfWpR/swyJc/7WbFV7A4UoYO6Y/ng60z/uEyXPMFASUdArNNkwefE22L3Y3kejFLFl/sQW7wErLS3F6tWrkZ2djW3btkEURVx//fXIysrCLbfcgqgo717GbE9vwAYN+3tqb0CHw4HCwkLk5+c3CY0nTpyAxWJxHdu7d28MGjSoyUef5BRkvPEDZACyoxb2M0eanOfigi91rxTXPuuRnSj9/EXETvk9Iq5o/p5LAcChp29wrSrX1NTg/PnzOHfuHM6dO+f258YfxcXFkCT3n/X4+PgmQbO58BkW5l5gNvHlHBRcdK+o5VAOTJvedFWyXxwoT791L5RRCUi844U2vxf947TY9mhGm8eR/+O1DyIKesmxWjyTmY4Fqw96bcxnM9NbDWLx8fG47777cN999+H8+fNYtWoVDAYDfvOb30CtVuPGG2+EXq/HzTffjPDwcI/mwt6ATdlsNpw4caLZlcZTp07B4ai/HK1QKKDT6TBo0CCMGzcOc+fOdYXGgQMHtvq96UzBl2vfgc0QlBqEXTqx5fHjtG63KISEhCAlJQUpKSmtvnen04mSkhK3kNk4fB47dgzbtm3DuXPnYLfb3V4bHh7uCpcJfZNR0D8LaHRPaFuV7I6qEjgri6BNvQrl2z6AZf9mSNVVUMX1Q+SYWxE+/Dq34/2l8Ig8x+8gEfUIWaN1KLXYvdL65LEpQzrU8iQxMREPPfQQHnroIZw+fRqffvopDAYD7rjjDoSGhmLq1KnQ6/W46aaboNV2bLWwp/YGlGUZJpPJ7XJ049B47tw517GhoaGukJiZmem20ti/f/9OXw7uTMEXUB+6ak7+iLD0SRBDmm85pRAFZKR17v5bhUKBxMREJCYm4vLLW25kKcsyKisrWwyepyrr3MIk0HYlu9NcBgCwHPoayoh4xE65H6JGC8tPG1H25WLIzjpEXPbrX+YA/yg8Is/xkjcR9SiNV/M6EgQUogClKODZzHSv9c87deoUVq5cCYPBgB9++AFhYWHIzMyEXq/Hr3/9a2g0mlZfn73P6NVV15dmDver3oCSJOH06dPNrjIeP34clZWVrmPj4uJcITE1NdUtNCYmJnrtUZuN5RWZcf1r270+boMt8yf4tGDlR2M5Zry1y/W59ci3KP3iZfT57T+hjq//OSldv9jtknfN6Z9R9NFjgEKJpPvegTKqPhTLsozz7/8JTlsl+j34vtt51jwwDpfrYrrnTVGX4QolEfUoWaN1uGZQfIfvNxw3MM7r9xumpKTg8ccfx+OPP468vDxXuPzkk08QGRmJ6dOnQ6/XY/LkyVCr3dvKBEtvQLvdjpMnTzZbAHPy5EnU1tYCAARBQHJyMgYNGoQrrrgCs2bNcguN3r4ntT0G947A+NT4Liv48nX1s9uTptpZya4IrZ+zKrafK0wC9d+/kIFXoGr3p3BaK6AIi272PBS4uEJJRD2WqyL6WDGMZc1URMdpkZGWgDljdd36y/2///0vVq5ciezsbBw9ehQxMTGYOXMm9Ho9MjIyoFQqXZXr1eeOo2L7h6gtKYBkq4SgVEMZm4SIK6YhfFjzxQ6yLKNoxQLYTx9GxBVTETvlAQDNV657Q0VFRYurjKdPn0bDryGNRoMBAwY0u8qYkpLS5oqtL/ii4Ku7WO0ODHt6I2S0v5K914z/h8LFs6GM7oO+9y5x21++9X1U7VmFfg+vgEJb/w+AiwuPKHAxUBIRwT97NsqyjAMHDsBgMMBgMODEiRPo1asXpsz6DXZGTgIA1BQcgPXn7dD0S4cyIg5SXQ2sh7fC9vN2RI2fg+hrspqMW/X9F/UrRRaTW6Bs0NFLrbIs49y5cy2GxrKyMtex0dHRzVZNDxo0CElJSV3eV7ErBPOtBw1V3h2pZC9Z9zJsR75F0n3/gjK6N4ALPyPv/RFybTWS7n/X9XpWeQcP/pOAiAj1j2n0t8IAQRAwcuRIjBw5Es8//zy+//57GAwGrDwuA+FOCKICIf1HNKki1qZehXOVRbDs39gkUDoqilCx7UPET52PkjVNn22uEAV8tMfYpDdgXV0dCgoKWmy1U11d7Tq2b9++GDRoEC699FLcfPPNbquN3dGLs7v5suCrq7kKjzpQyR49fg5qjn+HopVPIfra2yFqtDDv34S64pOIn/6E6zhPCo/I/zBQEhEFAEEQcOWVV+LKK69E7ss5TZ4jfTFFaCQka0WT7WUbliAk5TJoh4xr9nVOScb6H04g7Oh/3EKj0WiE0+kEUN9QOyUlBYMGDcKkSZNw7733ugLjgAEDOlypHgweyhiM+HCN3xR8eUtbT5pqjiqmD3rPeQkVWz9A2YYlgOSEKmEAet32JLSpV7mOc0oy5oz1r/dLncdASUQUQCx2BwqbCZOyLAGyDKnGAtuRnag++QNir7/f7Rjz/o2wnzuGvvPebPUcJTUCnv3nSxik64dBgwbhtttuc1tlTE5OhkKh8Or7Cgb+VPDlLW0VHjX36FIAUPdKQcKsp1oc118Kj8h7GCiJiAJIQZkVzUUU08Y3YflpQ/0nCiViJ/8PIi6/0bW/4dnLMRm/hTIirtVzCIKAPQePIz3Jv24BCATJsVosv3eM3xZ8dcaiGcMxefE2r1ayK0UBi2YM99p45HsMlEREAaS2hWriqKtnI3zkDZBsFbDl74Vp89uQ6moQNWYmAMC04Q2oEwYgfOQN7TuP03tVyz3R4N4ReDozHU8j3S8LvjrCF0+aosATOD/RRETUYs8+ZVSCq+9f6KDRAICKbR8gfPh1qDEeRPXJH5A452+Q7e4rnLLTAanGAkEVAkHxy68E9gb0Hn8s+OqoYC48Iu9goCQiCiApcWEQgGYvezem6ZMGy49fwVFxHnUlBYDkxPkPH2lynGX/Rlj2b0SvmX+GNu1qAPWXZFPimn8kIPVcwVp4RN7BQElEFCDMZjPWrl0LlV1Erab1Fa+aggOAIEIZnYjw4ZMRomt6v1rRJwsROngsIq/MhKpXf9d2XZw2oC7JUvcJxsIj8g7+jUFE5MdsNhvWr18Pg8GA//znP6ipqcGlc59CXdKVkCGg7KvXIWq0UPdJgyIsGk5bFWxHd8L28w5EjplZ/0QSbZSrwfTFlBFxbj0E2RuQ2hKMhUfkOT4ph4jIz9TU1GDDhg0wGAz44osvYLVaceWVV0Kv12P27Nmwa2Jw/WvbAQCWA5thObAFdWWFkOxWiKoQqC4U37T06MUGBS9O88qTcogCvfCIPMdASUTkB2pra7FlyxYYDAZ8/vnnqKqqwogRI6DX66HX6zFo0CC34xue5e3NVi5d9SxvIgp+DJRERD7icDiwdetWGAwGrF69GiaTCUOHDkVWVhb0ej2GDh3a4msLTTZMXrwN9hbaCHWGRiliy/yJvM+NiDqMgZKIqBtJkoSdO3fCYDBg1apVKC4uxsCBA10hcvjw4RAEoV1jZe8zerU34Eszh7MCl4g6hTc4EBF1MVmWkZubC4PBgE8//RRnzpxBcnIy7rrrLuj1eowaNardIbIx9gYkIn/BFUoioi4gyzJ+/PFHGAwGrFy5EqdOnUJiYiJmz54NvV6PsWPHQhS90zw8e5+RvQGJyKcYKImIvOjQoUMwGAwwGAzIy8tDfHw8brvtNuj1eowfPx4KhaJLzltosnW4N+D41Hj2BiQir2CgJCLy0LFjx1wh8vDhw4iOjsbMmTOh1+vxq1/9Ckpl991dxN6AROQLDJRERJ1w8uRJrFy5EgaDAT/++CPCw8Mxffp06PV6TJkyBWq12tdTZG9AIuo2DJRERO10+vRpfPrppzAYDMjNzUVoaChuvvlm6PV63HjjjQgNDfX1FImIfIKBkoioFUVFRVi1ahUMBgN27NgBtVqNm266CXq9HtOmTUN4eLivp0hE5HMMlEREFykrK8Pq1auRnZ2NrVu3QhRFTJkyBXq9HrfccguioqJ8PUUiIr/CQElEBKCyshKff/45DAYDNm/eDEmS8Ktf/Qp6vR4zZsxAXFycr6dIROS3GCiJqMeyWCz44osvYDAY8NVXX6Gurg7jx4+HXq/Hrbfeit69e/t6ikREAYGBkoh6lOrqavznP/+BwWDA+vXrUV1djbFjx0Kv12PWrFlISkry9RSJiAIOAyURBT273Y5NmzbBYDBg7dq1sFgsuOKKK6DX6zF79mykpKT4eopERAGNgZKIglJdXR2++eYbZGdnY82aNaisrMSwYcOg1+uh1+sxePBgX0+RiChoMFASUdBwOp3Yvn07srOz8dlnn6GsrAxpaWmuEJmenu7rKRIRBSUGSiIKaJIkYffu3cjOzsaqVatw/vx5pKSkuELkZZddBkEQfD1NIqKgxkBJRAFHlmV89913rudnnz59GklJSZg9ezaysrIwevRohkgiom7Eh7oSUUCQZRkHDhxwhcgTJ04gISEBs2bNgl6vxzXXXANRFH09TSKiHokrlETk137++WcYDAZkZ2fj6NGjiI2Nxa233oqsrCxMnDgRCoXC11MkIurxGCiJyO/k5+e7ViIPHjyIyMhIzJgxA3q9HpMnT4ZKpfL1FImIqBEGSiLyCwUFBVi5ciUMBgO+//57hIWFITMzE3q9HjfccANCQkJ8PUUiImoBAyUR+czZs2fx6aefwmAwYPfu3QgJCcHUqVORlZWFm266CVqt1tdTJCKidmCgJKJuVVJSglWrVsFgMGD79u1QKpX49a9/Db1ej8zMTERERPh6ikRE1EEMlETU5UwmE9asWQODwYBvvvkGADB58mRkZWVh+vTpiI6O9u0EiYjIIwyUREHGanfgVJkVtQ4JaqWIlLgwhGm6v0NYVVUV1q5dC4PBgE2bNsHhcGDSpEnIysrCzJkzER8f3+1zIiKirsFASRQE8orMWJFrRM7RYhhNNjT+n1oAoIvVImNIAu4co8Pg3l13SdlqteLLL79EdnY2/vOf/8But+Oaa65BVlYWbrvtNiQmJnbZuYmIyHcYKIkCWKHJhoVrDmJHfikUogCn1PL/zg37x6fGY9GM4UiO9U7BS01NDb766isYDAZ88cUXsNlsGD16NLKysjBr1iwkJyd75TxEROS/GCiJAlT2PiOeWncYDkluNUheTCEKUIoCnslMR9ZoXafOXVtbi82bN8NgMODzzz+H2WzGyJEjkZWVhdmzZ2PgwIGdGpeIiAITAyVRAFqSk4dXNh3zeJxHp6ThoYzB7TrW4XAgJycHBoMBq1evRnl5OS655BJkZWVBr9djyJAhHs+HiIgCEwMlUYDJ3mfEgtUHvTbeSzOHQ9/CSqUkSdixYwcMBgNWrVqFkpISDBo0yBUihw0bBkEQvDYXIiIKTAyURAGk0GTD5MXbYHdITfbZzx5FxY6PYD9zBJBlqPsMRvSEuQjpd2mrY2qUIrbMn+i6p1KWZezZswcGgwGffvopzp49C51OB71eD71ejyuuuIIhkoiI3DBQEgWQuctysetEWZN7Ju3njuH8R09A0ycNkVdNB2QZlbmfobboBBLvWARN0iUtjqkQBVw9MA7zL1PCYDBg5cqVKCgoQJ8+fTB79mzo9XqMHTuWIZKIiFrEQEkUIPKKzLj+te3N7isy/AW1xSeQdP9SiKr6Z15LdhvOvD0PqtgkJM59uc3xz7x7P6KFGsyaNQt6vR7XXnstFAqFV98DEREFp+7vdkxEnbIi19hiayD7mZ8ROuhKV5gEAFGjRUjyMNiO7YLDYoIyPLbFsQVZwtxn3sE7/3M9lEr+tUBERB0j+noCRNQ+OUeLW2wPJDvrIChUTXco67fVlZxqdWxZEHG8OoRhkoiIOoWBkigAWOwOGE22Fver4nSwnz0KWf6lWEeWnKg9exQAIFWb2zyHscwGq93h+WSJiKjHYaAkCgAFZVa0drNz5KhpcJjOwLTpbTjMpXBUlcC04Q04KovrD2hHQY0M4FSZ1SvzJSKinoXXt4gCQG0zbYIaCx85Bc7qKlTuMsDy438AAJqkoYgcMxNVe1ZBER7nlfMQERE1h4GSKAColW1fTIgaexsir7wFdeVnIKq1UEYloGzDEgiqEKgTU712HiIioosxUBIFgJS4MAhAq5e9AUBQqqDulQIAcFQWw/rzDoSPvAGiStPmOYQL5yEiIuooBkqiABCmUUIXq0VBC4U5tSWnYDu6C+rEwRCUKtQVnUDlnlVQxfRF9IQ57TqHLk6LMA3/SiAioo7jbw+iAFBdXY1Y+zmcksIhiE2bjQsKFWoKDsD83ReQ6qqhjOyFiMtvROTYWRDVIc2M6E4hCshIS+iKqRMRUQ/AJ+UQ+TGr1Yq3334bL7/8MsqlEPS5940uO9eW+ROQmhDRZeMTEVHw4h34RH6oqqoKL7zwAlJSUrBgwQJMnToVh3dtwfjUeChE7z5TWyEKGJ8azzBJRESdxhVKIj9SXl6Of/7zn/jHP/4Bq9WKe+65B0888QRSUlIAAIUmGyYv3ga7F9v7aJQitsyfiORYrdfGJCKinoUrlER+oLS0FH/+85+RkpKCF198EXPnzsXx48fx1ltvucIkACTHavFMZrpXz/1sZjrDJBEReYRFOUQ+dP78efz973/HW2+9BQB44IEH8MgjjyAxMbHF12SN1qHUYscrm455fP7HpgyBfrTO43GIiKhn4yVvIh84c+YM/va3v+Gdd96BSqXCH/7wB8yfPx/x8fHtHiN7nxFPrTsMhyTDKbX/f2OFKEApCng2M51hkoiIvIKBkqgbFRQU4MUXX8S///1vaLVa/OlPf8LDDz+MmJiYTo1XaLJh4ZqD2JFfCoUotBosG/aPT43HohnDeZmbiIi8hoGSqBvk5+fjhRdewIcffojo6Gj87//+Lx588EFERkZ6Zfy8IjNW5BqRc6wYxjKb2xN1BNQ3Lc9IS8CcsTpWcxMRkdcxUBJ1oSNHjuD555/Hxx9/jF69euGxxx7D/fffj7CwrnvEodXuwKkyK2odEtRKESlxYXwCDhERdSkGSqIucPDgQfz1r3/Fp59+ir59++KJJ57AvHnzEBoa6uupEREReR3bBhF50Q8//ICZM2dixIgRyM3NxZtvvonjx4/jD3/4A8MkEREFLQZKIi/Izc3FtGnTMGrUKBw4cADLli1DXl4e7r//fmg0Gl9Pj4iIqEsxUBJ5YMeOHZgyZQrGjh2L48ePY/ny5Thy5AjuueceqFQqX0+PiIioWzBQEnWQLMv4+uuvMWnSJEyYMAHnz5+HwWDAoUOHMGfOHCiVLIAhIqKehYGSqJ1kWcZXX32Fa665BpMnT4bZbMaaNWvw008/Yfbs2VAoFL6eIhERkU8wUBK1QZZlrF27FldddRVuuukmyLKML7/8Et999x2mT58OUeT/RkRE1LPxNyFRCyRJwqefforLLrsM06dPh1arxebNm7Fr1y7cdNNNEATB11MkIiLyCwyURBdxOBxYsWIFhg0bhtmzZyMhIQHbtm3Dtm3bMHnyZAZJIiKiizBQEl1QV1eH9957D5dccgnmzJmDAQMGYNeuXdi8eTMmTJjg6+kRERH5LQZK6vHsdjv+9a9/IS0tDffccw+GDRuG7777Dl9++SWuvvpqX0+PiIjI77G/CfVY1dXVWLp0Kf72t7/hzJkzmDVrFtauXYsRI0b4empEREQBhYGSehyr1Yq3334br7zyCoqLi3H77bfjz3/+My655BJfT42IiCggMVBSj1FVVYU33ngDr776KioqKnDXXXdhwYIFGDx4sK+nRkREFNAYKCnolZeX45///Cf+8Y9/wGq14p577sETTzyBlJQUX0+NiIgoKDBQUtAqLS3Fa6+9htdffx21tbX43e9+h8cffxz9+vXz9dSIiIiCCgMlBZ2ioiL8/e9/x5tvvglZlvHAAw/g0UcfRWJioq+nRkREFJQEWZZlX0+CyBvOnDmDl19+Ge+88w6USiUeeughzJ8/H7169fL11IiIiIIaAyUFvIKCArz00ktYtmwZtFot/vjHP+Lhhx9GbGysr6dGRETUIzBQUsA6fvw4XnjhBXzwwQeIiorC//7v/+LBBx9EVFSUr6dGRETUozBQUsA5cuQIFi1ahI8//hjx8fF49NFHcf/99yM8PNzXUyMiIuqRWJRDAePQoUP461//ipUrV6Jv37549dVX8bvf/Q6hoaG+nhoREVGPxmd5k9/78ccfMXPmTAwfPhx79uzBm2++iePHj+Phhx9mmCQiIvIDXKGkTrHaHThVZkWtQ4JaKSIlLgxhGu/+OOXm5uK5557Dl19+iUGDBmHZsmWYO3cuVCqVV89DREREnmGgpHbLKzJjRa4ROUeLYTTZ0PjmWwGALlaLjCEJuHOMDoN7R3T6PDt27MBzzz2HzZs3Y+jQoVi+fDmysrKgVPLHlYiIyB+xKIfaVGiyYeGag9iRXwqFKMAptfwj07B/fGo8Fs0YjuRYbbvOIcsycnJy8Nxzz2Hr1q0YNmwYnnzySdx6661QKBTeeitERETUBXgPJbUqe58Rkxdvw64TZQDQaphsvH/XiTJMXrwN2fuMrR4vyzI2bNiAa6+9Ftdddx0qKyuxevVq7N+/H7Nnz2aYJCIiCgAMlNSiJTl5WLD6IOwOqc0geTGnJMPukLBg9UEsyclrsl+WZaxbtw5XXXUVbrzxRjidTqxfvx7ff/89ZsyYAVHkjyYREVGg4G9talb2PiNe2XTMK2O9sukYDBdWKiVJwqpVq3D55ZfjlltuQWhoKDZt2oTdu3dj6tSpEATBK+ckIiKi7sN7KKmJQpMNkxdvQ0X+j7AezoH9zBE4zSUQNWFQJw5G1LW3Q5OY6jq+6rt1sP53Gxzl5yDV2qAIi4EmaSiixmVB3as/AECjFPGHQRV46+W/4r///S+uu+46PPnkk5g4caKv3iYRERF5CQMlNTF3WS52nSjD+c8WwVltRtjQa6GKT4bTVomqvWtQez4fCbOfRWjKSABAxY4VgCBAnTAAYkg4HBXnUblnFZzmMvS5+zWo4voBkhPVBfsxqmo3nnzySVx99dU+fpdERETkLQyU5CavyIzrX9sOAHBaK6AIi3bbL9VW48y/fgd1fH/0vv35FsepKy3E2aUPIGpcFqInzHFt3zJ/AlITOt9SiIiIiPwP76EkNytyjVCI9fcxXhwmAUBUh0IVp4PDXNrqOKI28sIffqnSVogCPtrTetU3ERERBR4GSnKTc7S41YpuqcaK2qLjUMXrmuyTJSdkRx3qygpR9tXrELXRCB8x2bXfKcnIOVbcJfMmIiIi3+GjR8jFYnfAaLK1eoxp81uQ62oQNU7fZJ/x77cBzjoAgDI2CYl3vABlZC/3Y8pssNodXn9MIxEREfkOVyjJpaDMitZuqK3YvhzWw1sRc908tyrvBolzX0bi3FcQd/MjENWhOP/J/0NtSYHbMTKAU2VW706ciIiIfIqBklxqHVKL+yp2fozKXQZET7gLkaNubvYYTWIqNElDEZ6egd63LwJkoGLbhx06DxEREQUeBkpyUSub/3Go2PkxKnd+jKhr70DUuNntGkvUaKGK64e68jPtPg8REREFJv5mJ5eUuDBc/Jyaim8/qQ+T4/SIvvaOdo/ltFWiruQUVNF93LYLF85DREREwYOVEeQSplFCF6tFwYXCnKrc1ajcsQIhA0chdNBo2M8ccTtekzQUUo0VRYb/Q9ilk6CM6QtBqYaj/AyqvlsH2VmHqItCqC5Oy4IcIiKiIMPf7OQmY0gClucWwCnJsOXvBQDUnPge50983+TY/gvWQ1CqoU4YAMtPG+Awl0J21EIRFoMQ3XBEzlgIdaP2QgpRQEZaQre9FyIiIuoefFIOuWn8pJyuwCflEBERBR/eQ0luBveOwPjUeNfTcrxFIQoYnxrPMElERBSEGCipiUUzhkPp5UCpFAUsmjHcq2MSERGRf2CgpCaSY7V4JjPdq2M+m5mO5FitV8ckIiIi/8BASc3SnP4eFduXe2Wsx6YMgX5002d/ExERUXBgoKQmvv32W9xxxx24oZ+MF6YPg0YpdvieSoUoQKMU8dLM4Xgwo+ljGomIiCh4sMqb3Pz888+45pprMHLkSGzYsAEajQaFJhsWrjmIHfmlUIgCnFLLPzIN+8enxmPRjOG8zE1ERNQDMFCSy9mzZ3H11VcjKioK27dvR3R0tNv+vCIzVuQakXOsGMYyGxr/4Aiob1qekZaAOWN1rOYmIiLqQRgoCQBQWVmJCRMmoLy8HLt370ZSUlKrx1vtDpwqs6LWIUGtFJESF8Yn4BAREfVQDJQEu92OG2+8ET/++CN27tyJ9HTvVngTERFRcOOSUg8nSRLuvvtu7Nq1C5s3b2aYJCIiog5joOzhHnvsMRgMBqxatQrjx4/39XSIiIgoADFQ9mCvvvoqXn31Vbz++uuYOXOmr6dDREREAYr3UPZQ2dnZuP3227FgwQK88MILvp4OERERBTAGyh7om2++wa9//WvcfvvteP/99yEI3n1uNxEREfUsDJQ9zP79+zFhwgSMHTsW69evh0ql8vWUiIiIKMAxUPYgBQUFuPrqq9GnTx9s3boVERFsPk5ERESeY6AMYB1pLm4ymXDttdeipqYGu3btQmJiYjfPloiIiIIVq7wDjOvxh0eLYTQ18/jDWC0yhiTgzjE6DO5dvwJZXV2NzMxMlJSU4Ntvv2WYJCIiIq/iCmWAKDTZsHDNQezIL4VCFOCUWv62NewfnxqP5zIvxfz77sLGjRuRk5ODMWPGdOOsiYiIqCdgoAwA2fuMeGrdYTgkudUgeTGFKEB2OlC68S18/NyDmDZtWhfOkoiIiHoqXvL2c0ty8vDKpmOdeq1TkiFDROyvH8KpsDQvz4yIiIionujrCVDLsvcZOx0mGzT0mHxl0zEY9hm9MS0iIiIiN7zk7acKTTZMXrwN5jP5qNj+IWpLCiDZKiEo1VDGJiHiimkIH5bhOr7qu3Ww/ncbHOXnINXaoAiLgSZpKKLGZUHdqz8AQKMUsWX+RCTHan31toiIiCgIMVD6qbnLcrHrRBmsJ/fD+vN2aPqlQxkRB6muBtbDW2H7eTuixs9B9DVZAICKHSsAQYA6YQDEkHA4Ks6jcs8qOM1l6HP3a1DF9YNCFDBuYByW38vCHCIiIvIeBko/lFdkxvWvbW/1mHMfPgKnxYR+v3+vxWPqSgtxdukDiBqXhegJc1zbt8yfgNQENjUnIiIi7+A9lH5oRa4RCrH152srQiMhCK1/+0Rt5IU/KH55nSjgoz28l5KIiIi8h4HSD+UcLW7SHkiWJciSE05bJcw/fInqkz8gcuxtTV4rS07IjjrUlRWi7KvXIWqjET5ismu/U5KRc6y4y98DERER9RxsG+RnLHYHjCZbk+2mjW/C8tOG+k8USsRO/h9EXH5jk+OMf78NcNYBAJSxSUi84wUoI3u5H1Nmg9XuaPExjUREREQdwXso/czhs5WY+vrOJtsdlcVw2ioh2Spgy98Ly08bET3pbkSNmel2nP18PuB0oK7iHMz71sJhLkXvrOddld4NvvzDtUjvG9Wl74WIiIh6Bi5R+Zlah9TsdmVUApRRCQCA0EGjAQAV2z5A+PDroND+Egw1ian1/00aCm3qGJz5132o2PYhEm57sl3nISIiIuoo3kPpZ9TK9n1LNH3SAMkJR8X5Fo8RNVqo4vqhrvxMp89DRERE1BamCj+TEheG1uu769UUHAAEEcroxBaPcdoqUVdyCqroPm7bhQvnISIiIvIGXvL2M2EaJXSxWhRcKMwp++p1iBot1H3SoAiLhtNWBdvRnbD9vAORY2ZCoY2CVGNFkeH/EHbpJChj+kJQquEoP4Oq79ZBdtYh6to73M6hi9OyIIeIiIi8hqnCD2UMScDy3AI4JRmapKGwHNgCy8GvIdmtEFUhUCUMQNy0R1yPXhSUaqgTBsDy0wY4zKWQHbVQhMUgRDcckTMWQh2vc42tEAVkpCX46q0RERFREGKVtx9qz5NyPMEn5RAREZE38R5KPzS4dwTGp8a3+bScjlKIAsanxjNMEhERkVcxUPqpRTOGQ+nlQKkUBSyaMdyrYxIRERExUPqp5FgtnslM9+qYz2amIzlW69UxiYiIiBgo/VjWaB0enZLmlbEemzIE+tG6tg8kIiIi6iAW5QSA7H1GPLXuMBySDKfU/m+XQhSgFAU8m5nOMElERERdhoEyQBSabFi45iB25JdCIQqtBsuG/eNT47FoxnBe5iYiIqIuxUAZYPKKzFiRa0TOsWIYy2xw/+bJ6B8Xhoy0BMwZq2M1NxEREXULBsoAZrU7cKrMilqHhNtmTEfmdVdj8d9e9PW0iIiIqIfhk3ICWJhGifS+UQCA/lEKnDOe8u2EiIiIqEdilXeQ0Ol0MBqNvp4GERER9UAMlEGCgZKIiIh8hYEySOh0Opw9exZ1dXW+ngoRERH1MAyUQUKn00GWZZw5c8bXUyEiIqIehoEySPTv3x8AUFBQ4OOZEBERUU/DQBkkkpOTAYD3URIREVG3Y6AMElqtFvHx8QyURERE1O0YKIMIK72JiIjIFxgogwgDJREREfkCA2UQYaAkIiIiX2CgDCINgZKPZyciIqLuxEAZRHQ6HSwWCyoqKnw9FSIiIupBGCiDCHtREhERkS8wUAYRnU4HgL0oiYiIqHsxUAaRhIQEqNVqBkoiIiLqVgyUQUQURSQnJzNQEhERUbdioAwybB1ERERE3U3p6wmQd+l0OhzJP4nDZytR65CgVopIiQtDmIbfaiIiIuoagsymhUEhr8iMFblGrNp9BGZZA0EQXPsEALpYLTKGJODOMToM7h3hu4kSERFR0GGgDHCFJhsWrjmIHfmlUIgCnFLL386G/eNT47FoxnAkx2q7caZEREQUrBgoA1j2PiOeWncYDkluNUheTCEKUIoCnslMR9ZoXRfOkIiIiHoCBsoAtSQnD69sOubxOI9OScNDGYO9MCMiIiLqqVjlHYCy9xm9EiYB4JVNx2DYx6pwIiIi6jyuUAaYQpMN4x7+B8oPfAP7mSNwmksgasKgThyMqGtvhyYx1XVs6frFsB76uskYyth+SLrvbdfnGqWILfMn8p5KIiIi6hT2kgkwC9ccRPn3X8JhMyPyykyo4pPhtFWiau8anP/wESTMfhahKSNdxwtKDXrf/rzbGIJS7fa5Q5KxcM1BLL93TLe8ByIiIgouDJQBJK/IjB35pYi5/gEowqLd9oUOHIUz//odqnavdAuUEARokoa2Oq5TkrEjvxT5xWakJrClkL+x2h04VWZlX1EiIvJb/K0UQFbkGqEQBeCiMAkAojoUqjgdHObSTo2tEAV8tMeIpzPTPZwleUNDX9Gco8UwmmxofF8K+4oSEZG/YaAMIDlHi1tsDyTVWFFbdBwh/Ue4bZcdtSh8fQ4kWxUU4THQDh6LqPFzoAh1DyFOSUbOsWI8DQZKX2pPX1EZQIHJhuW5BXh/9yn2FSUiIp9joAwQFrsDRpOtxf2mzW9BrqtB1Di9a5s6YQDUCQOg6tUfAFBTeAjmfZ+jumA/+vxmMUR1qNsYxjIbrHYHL6f6SOO+ogDa7C3asH/XiTJMXryNfUWJiMhnmBwCREGZFS3Fi4rty2E9vBUx1/+PW5V35FXT3Y4LHXA51AkDUfr5C7D8tLHJfhnAqTIr0vtGeXXu1DZP+oo6LzS2X7D6IEotdvYVJSKibsdAGSBqHVKz2yt2fozKXQZET7gLkaNubnMc7ZCrIahCYD97tNn9v5o8BTFSJeLi4tr9odFoPHpvPZ23+4r2CtdAz5XKdmHBExGRd/BvzgChVjbtQV+x82NU7vwYUdfegahxszswmgwIQrN77szSAxWnUVZWhrKyMhw5csT15/LycjTXtjQsLKxJyIyNjW01hEZFRUEU2Ve/0GTDU+sOe3XMv6w7jHGD4nlPZQtY8ERE5H0MlAEiJS4MAuD65Vfx7Sf1YXKcHtHX3tHucWxHvoVcZ4em75Am+wQAzz3xcIsrNE6nExUVFa6AefGHyWRCWVkZSkpK3IJodXV1k7FEUWwzdPaE1dCFaw7CIcmoLTqBiu0forakAJKtEoJSDWVsEiKumIbwYRlur5GdDpi//wKWg1vgKD8HKFRQxycjOuNehPS7hH1FW8CCJyKirsNAGSDCNEroYrUoMNlQlbsalTtWIGTgKIQOGg37mSNux2qShsJRWYzSdS9De8kEqGL6AIKAGuNBmL9bB1W8DuEjb2hyDl2cttXLfQqFwhXsOqK6urrFENr4o72roR0Nov66GtrQVxQApBoLFBHxiL5kIpQRcZDqamA9vBVl6/8OR2URoq/JAgDIkhMlq59Hzen/ImrMrdAkDYVUZ0ft+XzIdTUA2Fe0OSx4IiLqWgyUASRjSAKW5xbAlr8XAFBz4nucP/F9k+P6L1gPQaOFGBaNqn2fQ7JWQJadUEYmIGLUzYi6ejZEdYjbaxSigIy0hC6Zd2hoKPr164d+/fq1+zXNrYY2rIA2/igtLcXRo0e9uhracGxISEgzM/Oehr6iTklGSP8RTVo+aVOvwrnKIlj2b3QFSvP361F94nskzvmbW8N6bepot9eyr+gvWPBERNT1GCgDyJ1jdHh/9ykk3vlim8cqQsKRMPPP7R7bKcmYM9Z/VmD8YTVUq9V2+JJ8R1ZDW+sr6vo6hEZCsla4Pq/6bh00yentevoR+4qy4ImIqLswUAaQwb0jMD41HjvzSyCj+aKazlCIAsYNjAuKy6OdXQ2trKxsM4R6YzW0YZ82KhYFpqavl2UJkGVINRbYjuxE9ckfEHv9/QAAR1UJnJVF0KZehfJtH8CyfzOk6iqo4vohcsytCB9+ndtYPb2vKAueiIi6T8/8TROgnE4nQg6tgVNxFUSlusVK7Y5SigIWzRjulbECkUKhQGxsLGJjYzF4cPsvaXZkNbThcr3JZIIsy1AlDEDfe15vMqZp45uw/LThwsSUiJ38P4i4/EYAgNNcBgCwHPoayoh4xE65H6JGC8tPG1H25WLIzjpEXPZr11g9va9oQ8FTc8z7N8L01esQVCHQPbLKtV2WZVj2b4Tlx69QV34WgqiEqld/RI65FdrU0Sx4IiJqgSA3d62P/I7FYsEdd9yBL7/8EvNeeA8byzt2Kbg1L80czst43USSJFRUVGDnz6fx8PrCJvsdlcVw2ioh2Spgy98Ly08bET3pbkSNmYma0z+j6KPHAIUSSfe9A2VU/T2vsizj/Pt/gtNWiX4Pvu823rW23UiJACIjI1v8iIqKQmRkJFQqVXd8CbpFXpEZ17+2vdl9DnMpzi59EKJKA8lucwuUFds/QuWubIRffiO0aeMgO2ph/n49ak79iF4zFkI7ZBwAYMv8CUGxoh8M2EuUyD/w/7oAcPbsWdx88804duwYvvjiC9x0000eFRo09tiUIQyT3ajhkvigAQoATQOlMirBFRRDB9UX2lRs+wDhw69zPX9dFdvPdQwACIKAkIFXoGr3p3BaK6AIi3bt+/H7fdh5+iiqqqpQVVUFp9PZ4txCQkJaDJtthdHGH/7Q2qlxwdPFTBveQEhyOsSQCNiOfuu2z3JwMzT9LkXcDQ+6toUOuByFr8+F5eDX0A4Zx4InP8BeokT+h4HSzx04cABTp06FLMvYuXMnRo4cCQB4KGMw4sM1rlYobRV3NKYQBShFAc9mpjNM+sjFfUVboumTBsuPX8FRcR7qxFQIqhbCWsOFhka3QQgA9n79pWu1RpZlVFdXu8JlZWWl68/NfTTsP3XqVJPtDoejxTmr1epOh9HGHyEhIRA6eVtHSwVPlkM5qCk8hL7z3kLF9uVNXygqIWrC3DYJSjUEpQqCUg2ABU++xF6iRP6LgdKPbdiwAbNmzcLgwYPxxRdfICkpyW1/1mgdrhkU3+ZfsA0a9o8bGMe/YH2scV/R1tQUHAAEEcroRAiiAqGDx8J25Fs4KoqgjO4N4EJQPPEDlNF9oND+cr/kxX1FBUGAVquFVqtFYmJip+cuyzJqamraFUYbf5w+fdptX2VlJerq6lo8j0ql6lQYVYWGN/t1dVorUP71u4iZdDeUkfHNnjPyykyUf7MM5v2boB1Sf8m7KvczyHYbIq785dGmPb3gyRfYS5TIv/FvQz/11ltv4Q9/+ANuvPFGfPLJJwgPD2/2uORYLZbfO+aXS0DHimEsa+YSUJwWGWkJmDNWx3u//ERDX1GnJKPsq9charRQ90mDIiwaTlsVbEd3wvbzDkSOmekKitHj56Dm+HcoWvkUoq+9HaJGC/P+TagrPon46U+4xu7KvqKCICA0NBShoaHo3bu3R2PZ7fZ2hdHGH+fOncPRo0fdttXU1LjGbLHgadObUMUmIfzym1qcT+ToWyAo1TBtegumr/4JABBDItDrtr8gpN+lruN6esFTd2MvUSL/x6IcPyNJEh5//HH8/e9/x8MPP4xXX30VCoWiQ2PwJvXA0LhwxHJgMywHtqCurBCS3QpRFQJVwgCEj7yhyaMXa0tOoWLrB6gpPARITqgSBiBqnB7a1KvcjutJhSO1tbUwm82orKzEdydL8fiWErf91iPfovSLl9Hnt/+EOr5+lap0/WLYjn7rVpRjObAZZRvfROSoaQgZOApwOmA59A2q83PRa8ZChA4c5Tp2bNUODIgUEBsbi5iYGLf/Nvw5LCys05ftqV72PiMWrD7otfFYhEjUNRgo/YjNZsOcOXOwdu1aLF68GA8//LCvp0RdbO6yXOw6Udahe2Db0tBXtKe2tjl8thJTX9/p+lyqrcaZf/0O4ekZiBqnd20v2/QWqvNy6yvjRSVkyYEzS36D8BGTETvlAbcxz69YAEdVCfo9sMy1LWr3mzAbf4bJZEJlZWWzc1GpVC2GzdaCaExMTFBV3XdWocmGyYu3odpqQeWubNQWnURt0XFI1VWIuuZ2RI+/0+34tto+AYBGKWLL/IlBe8sPFxTIV/hT5ifOnz+PzMxM/Pe//8XatWsxbdo0X0+JusGiGcMxefE2rwbKnt5X9OKCJ8lWBclagaq9a1C1d02T4wtfy0Lo4LGIGnsbZIcd6j5pTY7R9BkMe+EhSLXVENWhEADs2rjW9Yu64VGhJpMJ5eXlzf634c8nT57EDz/84Npmt9ubfR/h4eHtDqKNj4mIiAiaVdGGXqJStRnmnzZCnTAA2rSxsOzf1OzxlTtWuNo+RU+629X2qWTVM662T8HYS5RV7+QPGCj9wKFDhzB16lTU1dVh+/btuOKKK3w9JeomybFaPJOZ7tVLes9mpgft6kt7XFzwpAiPQe/bFzU5rnLPKtgLDyFh1tMQtZEQ1fVfM/vZI25PHZJlGfazRyGGhENQ1T/f/eKCp84+KhSob5DfniBqMplgNBpdn7f0yFCFQuFa5ezoyqg/tHxqkFdkxo78UgCAIioByX/KhiAIcNoqWwyU7Wn75JRk7MgvRX6xOeBvCWHVO/kTBkof27x5M2677TakpKRg/fr1SE5O9vWUqJtljdah1GJnX1EvalzwJCjVCOk/oskxloNfA4Lotk+bNg6WnzZCUKgQOuhKyI46WA99Dfvp/yJq/BwIguD1gqfQ0FAkJSU16eLQFkmSUFVV1SR0NhdIz5w5g4MHD7q22WzNdxfQarWdCqIdeYZ9ezXuJdruFdd2tH0CEBS9RFn1Tv6GgdKH3n33XTzwwAOYMmUKDAYDIiIC+1/L1HnsK+pdd47R4f3dpzr8uvjMR1H1/XpYD30Dy4HNEEQllLF9EXfzIwi7dBKA+l/Mc8b6/ussiiKio6MRHR2NgQMHdui1dru9xQB68cro4cOH3fY11xxfEARER0d36hJ9aGhos3NsqZdoa9rb9inQe4my6p38EYtyfECSJCxcuBAvvfQSHnjgAfzzn/+EUslsT+27hNWgYT8vYTWPBU/eJ8syzGZzs5fjW7tUX15eDrPZ3OyYGo2mSdiMjEvAzt7T3Rr1N3DaKnH6n3c2W5QDAOYfv4JpyzuAs77HqRgSgfjpTyA05TK34wQAh56+IeAKVlj1Tv4qsP5PCgLV1dW466678Nlnn+HVV1/Fn/70p6C5gZ48x76i3sOCJ+8TBMHVQL5///4dem1dXZ2rcKmtIHr0TBmQ2PG/Fy0HNsO05Z0mbZ9KPnuuSdunQOwlWmiy4al1h7065l/WHca4QfH8Byl5jCuUndDZtgzFxcXIzMzEgQMH8PHHH2P69OldP1kKeGwD0nlczQlMPxrLMeOtXc3ua2mF0llj6VDbJwCYVLsPlyTUN+jv3bs3EhISXP9Vq9XwN62tulef2g/r4RzYzxyB01wCURMGdeJgRF17OzSJqS2O2dNX3cl7+FupnTxty/Dzzz9j6tSpqK6uxvbt23HllVd229wpsIVplAG1iuJPWPAUmNTKjhf4OMpOt7vtU4NdO7dj7bEfUVFR0eQ1MTExbiHz4o/G27Xarl/da1z13hzLj/+Bs9qMyCszoYpPhtNWiaq9a3D+w0eQMPtZhKaMbPZ1wVT1Tr7FQNkGb7Rl+OabbzBz5kwkJycjJyenw5eKiKjzWPAUeC7uJdoeivBYAO1r+wTULwT8uH0TwjRK1NbWori4GEVFRa6Piz//+eefUVRUhNLS0ibtmsLDw5sNms1ti4yM7NRtTo2r3psTO+UBKMKi3baFDhyFM//6Hap2r2wxUALBUfVOvsdA2QpvtGWoPvQ17rvvPvzqV7/CypUrERXFlSai7pY1WodrBsV3uOBp3MA4Fjz5wMW9RAGg+vh3kOpqINdWAwDqygphPVL/RKTQQVdCGZXQrrZPDRr3ElWr1ejXrx/69evX5tycTidKS0vdwubFAXTv3r2ubQ6Hw+31Go2m3eEzNjbW1Y6prar3i8MkAIjqUKjidHCYW17ZBAK/6p38A++hbIEnbRkaK9/2IbKGR2PJkiV8lBqRH2DBU2B4et1hVy9RADj95j1wVhU3e2zS/cugjO5d3yboQtsnR2WRq+1TxKhpCLt0kitQKkQBc8f07/IVOVmWUV5e3mr4bPxRU1Pj9nqlUolevXohoW8yKiY/2WzVe2ukGitOv3UPQvqPQMLMP7d6bKBWvZP/YKBshrdv5H9x5nA2kCXyQyx48l95RWZc/9r2Lht/y/wJfvUPBlmWYbFYmg2feSXV2BkxvsNjln7xCqw/70DiXX9vtTCnwZd/uJb3a1On8W/OizS0ZagtOoGK7R+itqQAkq0SglINZWwSIq6YhvBhGa7jS9cvhvXQ103GUcb2Q9J9bwMAnlp3GNewLQOR32HBk/8a3DsC41Pju6yXqD+FSaC+JVNERAQiIiKQmuoe/n40lmNnC1XvLanYvhzWw1sRc/3/tCtMAkCtQ+rQOYgaY6C8yMI1B+GQZEg1Figi4hF9yUQoI+Ig1dXAengrytb/HY7KIkRfk+V6jaDUoPftz7uN0/gxXw5JxsI1B9mWgYioA9hLtF5Hq94rdn6Myl0GRE+4C5Gjbm77BZ08D1FjDJSNNG7LENJ/RJPn/2pTr8K5yiJY9m90C5QQBGiShrY4LtsyEBF1XHKsFs9kpnv1FqRnM9MD7mpRR6reK3Z+jMqdHyPq2jsQNW52u88hXDgPUWfxnyONNLRlaI0iNBKC0PEvW0NbBiIiar+s0To8OqVpb8nOCNReog1V722p+PaT+jA5To/oa+/o0DkaV70TdQZ/ehppri2DLEuAXH8J3HZkJ6pP/oDY6+93P8ZRi8LX50CyVUERHgPt4LGIGj8HitBfViPZloGIqHPYSxTIGJLgVvV+sarc1ajcsQIhA0chdNBo2M8ccdvf2lU0hSggIy3Bq/OlnoeB8gKL3QFjo55nDUwb34Tlpw31nyiUiJ38P4i4/EbXfnXCAKgTBkDVq75ZeU3hIZj3fY7qgv3o85vFbk9lMJbZYLU7+K9AIqIO6um9RO8co8P7u0+1uN+WvxcAUHPie5w/8X2T/f0XrG/xtU5JxpyxgRu2yT+wbdAFh89WYurrO5tsd1QWw2mrhGSrgC1/Lyw/bUT0pLsRNWZmi2NZj3yL0s9fQMyv5iHyqulu+9iWgYjIMz21l2hrz/LuLD7Lm7yFS2UXtNQuQRmVAGVU/aWA0EGjAQAV2z5A+PDroNA2Hwy1Q66GoAqB/ezRdp+HiIjaZ3DvCDydmY6nkd6jeomy6p38GYtyLmhvuwRNnzRAcsJRcb6NI+Vmn2rAtgxERN7T0Ev0cl0M0vtGBW2YBH6pevemQKx6J//EdHNBQ1uGttQUHAAEEcroxBaPsR35FnKdHZq+Q9y2sy0DERF5glXv5K+C959yHdTQlqHgQmFO2VevQ9Rooe6TBkVYNJy2KtiO7oTt5x2IHDMTCm0UHJXFKF33MrSXTIAqpg8gCKgxHoT5u3VQxesQPvIGt3OwLQMREXmKVe/kj5huGmnclkGTNBSWA1tgOfg1JLsVoioEqoQBiJv2iOvRi4JGCzEsGlX7PodkrYAsO6GMTEDEqJsRdfVsiOoQ19hsy0BERN7S06veyf+wyruRvCIzrn9te5eNv2X+hKCqOCQiIt/rqVXv5F8YKC/CtgxERBSoelLVO/kXBsqLFJpsmLx4G+xebO+jUYrYMn8iLzEQERFRUGKV90XYloGIiIioYxgom8G2DERERETtx0vercjeZ2RbBiIiIqI2MFC2odBkc7VlECBDbqX9eUNbhvGp8WzLQERERD0GA2U75RWZccujr8DZeyhqVZFsy0BERER0AQNlO1ksFsTExOD111/H3N/OY1sGIiIioguYgtpp165dcDgcmDhxIsI0SqT3jfL1lIiIiIj8Aqu822nbtm1ISEjA0KFDfT0VIiIiIr/CQNlOW7duxcSJEyEILRflEBEREfVEDJTtYLVasXfvXkyaNMnXUyEiIiLyOwyU7bB7927X/ZNERERE5I6Bsh22bt2K+Ph4XHrppb6eChEREZHfYaBsh23btvH+SSIiIqIWMFC2wWazITc3l/dPEhEREbWAgbINe/bsQV1dHe+fJCIiImoBA2Ubtm7diri4OKSnp/t6KkRERER+iYGyDdu2bcOECRMgivxSERERETWHKakV1dXV2LNnD++fJCIiImoFA2UrcnNzUVtby/sniYiIiFrBQNmKrVu3IiYmBsOHD/f1VIiIiIj8FgNlK3j/JBEREVHbmJRaUFNTg927d/P+SSIiIqI2MFC2YO/evbDb7QyURERERG1goGzB1q1bER0dzfsniYiIiNrAQNmChvsnFQqFr6dCRERE5NcYKJtht9uxa9cutgsiIiIiagcGymbs27cPNTU1vH+SiIiIqB0YKJuxdetWREVFYeTIkb6eChEREZHfY6BsxrZt2zB+/HjeP0lERETUDgyUF6mtrcW3337L+yeJiIiI2omB8iLfffcdqquref8kERERUTsxUF5k69atiIiIwGWXXebrqRAREREFBAbKizTcP6lUKn09FSIiIqKA0ONTk9XuwKkyK2odEkRI+Hbvd/jL/3vC19MiIiIiChiCLMuyryfR3fKKzFiRa0TO0WIYTTY0/gLIsozEcAVuHKHDnWN0GNw7wmfzJCIiIgoEPSpQFppsWLjmIHbkl0IhCnBKLb/1hv3jU+OxaMZwJMdqu3GmRERERIGjxwTK7H1GPLXuMByS3GqQvJhCFKAUBTyTmY6s0bounCERERFRYOoRgXJJTh5e2XTM43EenZKGhzIGe2FGRERERMEj6Ku8s/cZvRImAeCVTcdg2Gf0ylhEREREwSKoA2WhyYan1h326ph/WXcYhSabV8ckIiIiCmRBfcl77rJc7DpR1uSeSfP+jTB99ToEVQh0j6wCAMiSE+bv1qH65I+oKy2AVG2BIqoXtIPHImrsbRBDwgHU31M5bmAclt87ptvfDxEREZE/Cto+lHlFZuzIL22y3WEuRfk3/4YiPBaS/ZeVRtlRi4qdHyPs0okIHzkFitBI1BYdR+UuA6rz9yLxN4shqjRwSjJ25Jciv9iM1AS2FCIiIiIK2kC5ItfYbGsg04Y3EJKcDjEkAraj37q2C0o1kh5YBkVopGtbSP8RUET2QunnL8J2dBfCh2UAqF+l/GiPEU9npnfPmyEiIiLyY0F7D2XO0eImYdJyKAc1hYcQO+X3TY4XRIVbmGyg6ZMGAHCaS1zbnJKMnGPFXp4xERERUWAKykBpsTtgvKhwxmmtQPnX7yJm0t1QRsa3e6yaggMAAFV8f7ftxjIbrHaH55MlIiIiCnBBGSgLyqy4uNLItOlNqGKTEH75Te0ex2EuRfm296FOHIzQ1NFu+2QAp8qsnk+WiIiIKMAFZaCsdUhun1uPfAtb/l7E3vgHCILQrjGc1WYUr3wakIH46U9AEJp+qS4+DxEREVFPFJRFOWrlL+FPqq2GafNbiBx1M5ThsZBqLAAAWaq/XC3VWABRCVEd4nqNs8aC4uz/g9NiQu/bn4cqOrHN8xARERH1VEEZKFPiwiCg/rK0ZKuCZK1A1d41qNq7psmxha9lIXTwWCTc+n8ALoTJT/4MR2URemc9D3XCgGbPIVw4DxEREVFPF5SBMkyjhC5WiwKTDYrwGPS+fVGTYyr3rIK98BASZj0NUVtf3e0KkxXnkZD1V6gTB7V4Dl2cFmGaoPzyEREREXVI0CaijCEJWJ5bAKdSjZD+I5rstxz8GhBE1z6pzo5iw5OoLTqBmMm/AyQn7GeOuI4XtVFQxfQBUN+HMiMtoXveCBEREZGfC9pAeecYHd7ffardx0vWCtSeywMAlG95p8n+sGHXIX7afAD1fSjnjNV5ZZ5EREREga5HPsvbE3yWNxEREZG7oC5TXjRjOJRi+9oEtZdSFLBoxnCvjklEREQUyII6UCbHavGMl5+3/WxmOpJjtV4dk4iIiCiQBXWgBICs0To8OiXNK2M9NmUI9KN57yQRERFRY0F9D2Vj2fuMeGrdYTgkuUP3VCpEAUpRwLOZ6QyTRERERM3oMYESAApNNixccxA78kuhEIVWg6UIGRIEXBIDvDMvg5e5iYiIiFrQowJlg7wiM1bkGpFzrBjGMhsafwEE1Dctn5TWCxuX/B9qik/hhx9+gCgG/d0BRERERJ3SIwNlY1a7A6fKrKh1SFArRaTEhbmegLNr1y5cc801+Oijj3DnnXf6eKZERERE/qnHB8q2zJgxAz/99BOOHDkCjUbj6+kQERER+R1ex23DokWLYDQa8a9//cvXUyEiIiLyS1yhbId58+Zh7dq1OH78OCIjI309HSIiIiK/whXKdnj66adhsVjwyiuv+HoqRERERH6HgbId+vXrhz/+8Y949dVXcf78eV9Ph4iIiMivMFC20xNPPAGVSoXnnnvO11MhIiIi8isMlO0UExODhQsX4p133kFeXp6vp0NERETkN1iU0wHV1dVIS0vDuHHjYDAYfD0dIiIiIr/AQNlB7733Hu655x7s27cPV155pWt7aw3SiYiIiIIZA2UHOZ1OjBgxAomJiXhrxRp8vLcQOUeLYTQ18wjHWC0yhiTgzjE6DO4d4aspExEREXUpBspOeG/lOixYvR+hA66AQhTglFr+EjbsH58aj0UzhiM5VtuNMyUiIiLqegyUHZS9z4in1h2GvbYOEBXtfp1CFKAUBTyTmY6s0bounCERERFR92Kg7IAlOXl4ZdMxj8d5dEoaHsoY7IUZEREREfke2wa1U/Y+o1fCJAC8sukYDPuMXhmLiIiIyNe4QtkOhSYbJi/ehmqrBZW7slFbdBK1RcchVVch6prbET3+TrfjC16c1uJYyth+SLrvbWiUIrbMn8h7KomIiCjgsa9NOyxccxAOSYZUbYb5p41QJwyANm0sLPs3NXt84tymz/y2nz2K8q/fhTZtLADAIclYuOYglt87pkvnTkRERNTVGCjbkFdkxo78UgCAIioByX/KhiAIcNoqWwyUmqShTbaZf9oAQED4yCkAAKckY0d+KfKLzUhNYEshIiIiCly8h7INK3KNUIgCAEAQBAiC0OExJLsNtiM7odENgyqmr2u7QhTw0R7eS0lERESBjYGyDTlHi1vtM9ke1p+3Q66rca1ONnBKMnKOFXs0NhEREZGvMVC2wmJ3wGiyeT7Ogc0QNWEIG3JNk33GMhusdofH5yAiIiLyFQbKVhSUWeFpCXxtSQFqzx5FWPokCEp1k/0ygFNlVg/PQkREROQ7DJStqHVIHo9hOVBfuBM+8oYuPQ8RERGRrzBQtkKt9OzLIzvrYD2UA3ViKtS9B3bZeYiIiIh8iUmmFSlxYeh4TfcvbHm5kKqrED5iSovHCBfOQ0RERBSo2IeyFWEaJXSxWhQ0KsypPv4dpLoayLXVAIC6skJYj+wEAIQOuhKiKsR1rOXAZghKDcIundjiOXRxWoRp+G0gIiKiwMUk04aMIQlYnlvgah1UtvFNOKt+afVjO7ITtguBMun+ZRCj6wOlo6oENSd/RFj6JIghza9AKkQBGWkJXfwOiIiIiLoWn+XdhrwiM65/bXuXjb9l/gQ+KYeIiIgCGu+hbMPg3hEYnxrvelqOtyhEAeNT4xkmiYiIKOAxULbDohnDofRyoFSKAhbNGO7VMYmIiIh8gYGyHZJjtXgmM92rYz6bmY7kWK1XxyQiIiLyBQbKdsoarcOjU9K8MtZjU4ZAP1rnlbGIiIiIfI1FOR2Uvc+Ip9YdhkOSXZXf7aEQAKVCxLOZ6QyTREREFFQYKDuh0GTDwjUHsSO/FApRaDVYNuyPqTmPdU/ezsvcREREFHQYKD2QV2TGilwjco4Vw1hmQ+MvpID6puUZaQmQ87bj+Sf+gCNHjiA1NdVX0yUiIiLqEgyUXmK1O3CqzIpahwS1UkRKXJjrCTjV1dUYNGgQpkyZgvfff9+3EyUiIiLyMgbKbvKPf/wDjzzyCFcpiYiIKOgwUHaT6upqDBw4EL/+9a/x3nvv+Xo6RERERF7DtkHdJDQ0FE888QSWL1+O48eP+3o6RERERF7DFcpu1LBKedNNN2HZsmW+ng4RERGRV3CFshuFhobi8ccfxwcffIATJ074ejpEREREXsEVym5ms9kwcOBATJs2DUuXLvX1dIiIiIg8xhXKbqbVal2rlCdPnvT1dIiIiIg8xhVKH7DZbBgwYAAyMzPx7rvv+no6RERERB7hCqUPNKxSvv/++zh16pSvp0NERETkEa5Q+ojVasXAgQNxyy234J133vH1dIiIiIg6jSuUPhIWFobHHnsM7733XpNVSqvdgcNnK/GjsRyHz1bCanf4ZpJERERE7cAVSh+yWq0YMGAAZsyYgUeffQUrco3IOVoMo8mGxt8UAYAuVouMIQm4c4wOg3tH+GrKRERERE0wUPrYky++hnf3WxCScjkUogCn1PK3o2H/+NR4LJoxHMmx2m6cKREREVHzGCh9KHufEU+tO4ya2joIoqLdr1OIApSigGcy05E1WteFMyQiIiJqGwOljyzJycMrm455PM6jU9LwUMZgL8yIiIiIqHNYlOMD2fuMXgmTAPDKpmMw7DN6ZSwiIiKizmCg7GaFJhueWnfYq2P+Zd1hFJpsXh2TiIiIqL14ybubzV2Wi6+/+QZVB7+B/cwROM0lEDVhUCcORtS1t0OTmOp2vP18Pipy3oP97FFAVCCk/wjE/OpeqKITXccoRAHjBsZh+b1juvvtEBEREXGFsjvlFZmxI78Uld9/CUdlMSKvzETCrKcRM/k+OG0VOP/hI6g+td91fF1ZIYo+/n+QnQ70mv4E4m/6IxymMyj66HE4bZWu45ySjB35pcgvNvvibREREVEPp/T1BHqSFblGKEQBsVMegCIs2m1f6MBROPOv36Fq90qEpowEAFTsWAFBoULCrKcgaupbBKkTU3HmX/ehKnc1YjJ+63q9QhTw0R4jns5M77b3Q0RERARwhbJb5RwthlOSm4RJABDVoVDF6eAwlwIAZMmJ6vx90A4Z5wqTAKCMSkBI/+GwHdvt9nqnJCPnWHGXzp+IiIioOQyU3cRid8DYSuGMVGNFbdFxqOLr+0o6ys9BdtihThjQ5Fh1rwEX9te6bTeW2fiYRiIiIup2DJTdpKDMitaqn0yb34JcV4OocXoAgLO6/n5IMSS8ybFiaDgAGc4ai9t2GcCpMquXZkxERETUPgyU3aTWIbW4r2L7clgPb0XMdfOaVHlDEFp8nYCm+1o7DxEREVFXYKDsJmpl81/qip0fo3KXAdET7kLkqJtd2xWhEQAAqbpp5bZUbQEgQAwJa/d5iIiIiLoK00c3SYkLa7KeWLHzY1Tu/BhR196BqHGz3fYpY/pAUGpQW3KqyVi1Jacu7Fe7bRcunIeIiIioOzFQdpMwjRK62F+qtSu+/aQ+TI7TI/raO5ocL4gKhKZeBdvR3ZDsvxTzOCqLUWM8AO2QcU1eo4vTIkzDTlBERETUvRgou1HGkAQoRAFVuatRuWMFQgaOQuig0bCfOeL20SB6/B2QHXYUr3oW1ce/g+3oLhSvegaK0EhEXjXDbWyFKCAjLaG73xIRERERH73YnfKKzLj+te04v2IB7IWHWjyu/4L1rj//8ujFI/WPXtRdePRiTJ8mr9syfwJSEyK6ZO5ERERELWGg7GZzl+Vi14kyOCXvfdn5LG8iIiLyJV7y7maLZgyHUmy5FVBnKEUBi2YM9+qYRERERO3FQNnNkmO1eMbLz9t+NjMdyY0KfoiIiIi6EwOlD2SN1uHRKWleGeuxKUOgH63zylhEREREncF7KH0oe58RT607DIckd+ieSoUoQCkKeDYznWGSiIiIfI6B0scKTTYsXHMQO/JLoRCFVoNlw/7xqfFYNGM4L3MTERGRX2Cg9BN5RWasyDUi51gxjGU2NP6mCKhvWp6RloA5Y3VsDURERER+hYHSD1ntDpwqs6LWIUGtFJESF8Yn4BAREZHfYqAkIiIiIo+wypuIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgjDJRERERE5BEGSiIiIiLyCAMlEREREXmEgZKIiIiIPMJASUREREQeYaAkIiIiIo8wUBIRERGRRxgoiYiIiMgj/x+TBTrbbpqr3gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Select 10 random vertices from the networkx graph as terminal nodes:\n", "\n", "terminals = np.random.choice(random_graph.nodes, size=10, replace=False)\n", "stree = steiner_tree(random_graph, terminals)\n", "total_weight = sum(random_graph[u][v][\"weight\"] for u, v in stree.edges)\n", "print(\"Cost:\", total_weight)\n", "nx.draw(stree, with_labels=True)" ] }, { "cell_type": "markdown", "id": "caf651ab", "metadata": {}, "source": [ "We need to check if the solution is the optimal one or not. **Since we've used an heuristic (but fast) method, we don't have any proof about optimality**. We will see now how CORNETO can be used to get a better solution (which is also the optimal). Due to CORNETO's exact modeling of the problem, and the use of advanced mathematical solvers, we can not only achieve a solution superior to that of NetworkX, but also have the guarantees that the solution is optimal and no better network exists. " ] }, { "cell_type": "markdown", "id": "8d9119a2", "metadata": {}, "source": [ "## Exact Steiner Trees with CORNETO\n", "\n", "Now, we will do the same using CORNETO to convert the random graph to a corneto graph and then we will create an ILP steiner tree problem, which is exact, and we will solve it using the SCIP solver (open source)." ] }, { "cell_type": "code", "execution_count": 3, "id": "a1684e94", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", "
\n", " \n", " \n", " \n", " \n", "
Installed version:v0.9.2-beta.1 (latest: v0.9.1-alpha.6)
Available backends:CVXPY v1.4.1, PICOS v2.4.17
Default backend (corneto.K):CVXPY
Installed solvers:CBC, CLARABEL, COPT, CPLEX, CVXOPT, DIFFCP, ECOS, ECOS_BB, GLPK, GLPK_MI, GUROBI, MOSEK, OSQP, PROXQP, SCIP, SCIPY, SCS
Graphviz version:v0.20.1
Repository:https://github.com/saezlab/corneto
\n", "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import corneto as cn\n", "\n", "cn.info()" ] }, { "cell_type": "markdown", "id": "c24bc2c6", "metadata": {}, "source": [ "We need some code to convert the `networkx` graph into a `corneto` Graph:" ] }, { "cell_type": "code", "execution_count": 4, "id": "c7d5f9bc", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(50, 145)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from corneto._graph import EdgeType\n", "\n", "\n", "def to_graph(nx_graph, weights=None, directed=False):\n", " G = cn.Graph()\n", " if weights is None:\n", " weights = [e[2].get(\"weight\", 0) for e in nx_graph.edges(data=True)]\n", " for e, w in zip(nx_graph.edges, weights):\n", " etype = EdgeType.DIRECTED if directed else EdgeType.UNDIRECTED\n", " G.add_edge(e[0], e[1], weight=w, type=etype)\n", " return G\n", "\n", "\n", "G = to_graph(random_graph)\n", "G.shape" ] }, { "cell_type": "markdown", "id": "d8444117", "metadata": {}, "source": [ "We will use the `exact_steiner_tree` method. This method converts the original graph into a network flow problem equivalent to the steiner tree problem. The method returns the optimization problem and a new graph, which is an augmented version of the original graph, containing few extra edges required for the network flow problem" ] }, { "cell_type": "code", "execution_count": 5, "id": "501e5630", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'_flow': Variable((155,), _flow),\n", " '_flow_i': Variable((155,), _flow_i, boolean=True)}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from corneto.methods.steiner import exact_steiner_tree\n", "\n", "P, G_steiner = exact_steiner_tree(G, terminals)\n", "\n", "# The problem contains two vector variables, _flow and _flow_i.\n", "# _flow is the normal _flow variable to define a network flow problem on the graph.\n", "# _flow_i is an indicator binary vector. If _flow_i[i] = 1, then _flow[i] = 0.\n", "# The _flow_i variable contains the edges that are selected in the Steiner tree.\n", "P.symbols" ] }, { "cell_type": "code", "execution_count": 6, "id": "5c347f91", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "===============================================================================\n", " CVXPY \n", " v1.4.1 \n", "===============================================================================\n", "(CVXPY) Dec 24 03:28:24 PM: Your problem has 310 variables, 7 constraints, and 0 parameters.\n", "(CVXPY) Dec 24 03:28:24 PM: It is compliant with the following grammars: DCP, DQCP\n", "(CVXPY) Dec 24 03:28:24 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n", "(CVXPY) Dec 24 03:28:24 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n", "(CVXPY) Dec 24 03:28:24 PM: Your problem is compiled with the CPP canonicalization backend.\n", "-------------------------------------------------------------------------------\n", " Compilation \n", "-------------------------------------------------------------------------------\n", "(CVXPY) Dec 24 03:28:24 PM: Compiling problem (target solver=SCIPY).\n", "(CVXPY) Dec 24 03:28:24 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCIPY\n", "(CVXPY) Dec 24 03:28:24 PM: Applying reduction Dcp2Cone\n", "(CVXPY) Dec 24 03:28:24 PM: Applying reduction CvxAttr2Constr\n", "(CVXPY) Dec 24 03:28:24 PM: Applying reduction ConeMatrixStuffing\n", "(CVXPY) Dec 24 03:28:24 PM: Applying reduction SCIPY\n", "(CVXPY) Dec 24 03:28:24 PM: Finished problem compilation (took 3.863e-02 seconds).\n", "-------------------------------------------------------------------------------\n", " Numerical solver \n", "-------------------------------------------------------------------------------\n", "(CVXPY) Dec 24 03:28:24 PM: Invoking solver SCIPY to obtain a solution.\n", "Solver terminated with message: Optimization terminated successfully. (HiGHS Status 7: Optimal)\n", "-------------------------------------------------------------------------------\n", " Summary \n", "-------------------------------------------------------------------------------\n", "(CVXPY) Dec 24 03:29:09 PM: Problem status: optimal\n", "(CVXPY) Dec 24 03:29:09 PM: Optimal value: 4.100e+01\n", "(CVXPY) Dec 24 03:29:09 PM: Compilation took 3.863e-02 seconds\n", "(CVXPY) Dec 24 03:29:09 PM: Solver (including time spent in interface) took 4.489e+01 seconds\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "2\n", "\n", "2\n", "\n", "\n", "\n", "30\n", "\n", "30\n", "\n", "\n", "\n", "2->30\n", "\n", "\n", "\n", "\n", "46\n", "\n", "46\n", "\n", "\n", "\n", "47\n", "\n", "47\n", "\n", "\n", "\n", "46->47\n", "\n", "\n", "\n", "\n", "6\n", "\n", "6\n", "\n", "\n", "\n", "13\n", "\n", "13\n", "\n", "\n", "\n", "6->13\n", "\n", "\n", "\n", "\n", "17\n", "\n", "17\n", "\n", "\n", "\n", "6->17\n", "\n", "\n", "\n", "\n", "20\n", "\n", "20\n", "\n", "\n", "\n", "17->20\n", "\n", "\n", "\n", "\n", "8\n", "\n", "8\n", "\n", "\n", "\n", "8->47\n", "\n", "\n", "\n", "\n", "8->13\n", "\n", "\n", "\n", "\n", "14\n", "\n", "14\n", "\n", "\n", "\n", "8->14\n", "\n", "\n", "\n", "\n", "9\n", "\n", "9\n", "\n", "\n", "\n", "23\n", "\n", "23\n", "\n", "\n", "\n", "9->23\n", "\n", "\n", "\n", "\n", "34\n", "\n", "34\n", "\n", "\n", "\n", "9->34\n", "\n", "\n", "\n", "\n", "23->30\n", "\n", "\n", "\n", "\n", "33\n", "\n", "33\n", "\n", "\n", "\n", "23->33\n", "\n", "\n", "\n", "\n", "34->47\n", "\n", "\n", "\n", "\n", "41\n", "\n", "41\n", "\n", "\n", "\n", "34->41\n", "\n", "\n", "\n", "\n", "18\n", "\n", "18\n", "\n", "\n", "\n", "18->30\n", "\n", "\n", "\n", "\n", "48\n", "\n", "48\n", "\n", "\n", "\n", "18->48\n", "\n", "\n", "\n", "\n", "25\n", "\n", "25\n", "\n", "\n", "\n", "25->33\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P.solve(solver=\"SCIPY\", verbosity=1)\n", "G_steiner.edge_subgraph(np.where(P.symbols[\"_flow_i\"].value > 0.5)[0]).plot(\n", " orphan_edges=False\n", ")" ] }, { "cell_type": "code", "execution_count": 7, "id": "e3c48b24", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "41.0" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The cost is lower than the solution we obtained with networkx\n", "P.objectives[0].value" ] }, { "cell_type": "markdown", "id": "276e90a2", "metadata": {}, "source": [ "The Steiner problem is only guaranteed to yield a Steiner tree when we allow the solver to run until optimality is proven. However, we can enforce the tree structure more strictly by setting `strict_acyclic=True`. In this scenario, instead of using a standard network flow problem, it employs an acyclic network flow. While this approach introduces more variables, it also adds more constraints. These additional constraints exclude many structures that aren't trees from the feasible space, potentially aiding solvers in finding the optimal solution more quickly. Moreover, we can stop the solver early, even before optimality is confirmed, and still obtain a valid Steiner tree." ] }, { "cell_type": "code", "execution_count": 8, "id": "6e7589c9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "===============================================================================\n", " CVXPY \n", " v1.4.1 \n", "===============================================================================\n", "(CVXPY) Dec 24 03:29:11 PM: Your problem has 670 variables, 17 constraints, and 0 parameters.\n", "(CVXPY) Dec 24 03:29:11 PM: It is compliant with the following grammars: DCP, DQCP\n", "(CVXPY) Dec 24 03:29:11 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n", "(CVXPY) Dec 24 03:29:11 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n", "(CVXPY) Dec 24 03:29:11 PM: Your problem is compiled with the CPP canonicalization backend.\n", "-------------------------------------------------------------------------------\n", " Compilation \n", "-------------------------------------------------------------------------------\n", "(CVXPY) Dec 24 03:29:11 PM: Compiling problem (target solver=SCIPY).\n", "(CVXPY) Dec 24 03:29:11 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCIPY\n", "(CVXPY) Dec 24 03:29:11 PM: Applying reduction Dcp2Cone\n", "(CVXPY) Dec 24 03:29:11 PM: Applying reduction CvxAttr2Constr\n", "(CVXPY) Dec 24 03:29:11 PM: Applying reduction ConeMatrixStuffing\n", "(CVXPY) Dec 24 03:29:11 PM: Applying reduction SCIPY\n", "(CVXPY) Dec 24 03:29:11 PM: Finished problem compilation (took 9.297e-02 seconds).\n", "-------------------------------------------------------------------------------\n", " Numerical solver \n", "-------------------------------------------------------------------------------\n", "(CVXPY) Dec 24 03:29:11 PM: Invoking solver SCIPY to obtain a solution.\n", "Solver terminated with message: Optimization terminated successfully. (HiGHS Status 7: Optimal)\n", "-------------------------------------------------------------------------------\n", " Summary \n", "-------------------------------------------------------------------------------\n", "(CVXPY) Dec 24 03:29:45 PM: Problem status: optimal\n", "(CVXPY) Dec 24 03:29:45 PM: Optimal value: 4.100e+01\n", "(CVXPY) Dec 24 03:29:45 PM: Compilation took 9.297e-02 seconds\n", "(CVXPY) Dec 24 03:29:45 PM: Solver (including time spent in interface) took 3.385e+01 seconds\n", "Optimal value: 40.999999999997144\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "2\n", "\n", "2\n", "\n", "\n", "\n", "30\n", "\n", "30\n", "\n", "\n", "\n", "2->30\n", "\n", "\n", "\n", "\n", "46\n", "\n", "46\n", "\n", "\n", "\n", "47\n", "\n", "47\n", "\n", "\n", "\n", "46->47\n", "\n", "\n", "\n", "\n", "9\n", "\n", "9\n", "\n", "\n", "\n", "23\n", "\n", "23\n", "\n", "\n", "\n", "9->23\n", "\n", "\n", "\n", "\n", "34\n", "\n", "34\n", "\n", "\n", "\n", "9->34\n", "\n", "\n", "\n", "\n", "23->30\n", "\n", "\n", "\n", "\n", "33\n", "\n", "33\n", "\n", "\n", "\n", "23->33\n", "\n", "\n", "\n", "\n", "34->47\n", "\n", "\n", "\n", "\n", "41\n", "\n", "41\n", "\n", "\n", "\n", "34->41\n", "\n", "\n", "\n", "\n", "14\n", "\n", "14\n", "\n", "\n", "\n", "14->33\n", "\n", "\n", "\n", "\n", "17\n", "\n", "17\n", "\n", "\n", "\n", "20\n", "\n", "20\n", "\n", "\n", "\n", "17->20\n", "\n", "\n", "\n", "\n", "42\n", "\n", "42\n", "\n", "\n", "\n", "17->42\n", "\n", "\n", "\n", "\n", "18\n", "\n", "18\n", "\n", "\n", "\n", "18->30\n", "\n", "\n", "\n", "\n", "48\n", "\n", "48\n", "\n", "\n", "\n", "18->48\n", "\n", "\n", "\n", "\n", "25\n", "\n", "25\n", "\n", "\n", "\n", "25->33\n", "\n", "\n", "\n", "\n", "35\n", "\n", "35\n", "\n", "\n", "\n", "25->35\n", "\n", "\n", "\n", "\n", "35->42\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P, G_steiner = exact_steiner_tree(G, terminals, strict_acyclic=True)\n", "P.solve(solver=\"SCIPY\", verbosity=1)\n", "print(\"Optimal value:\", P.objectives[0].value)\n", "G_steiner.edge_subgraph(np.where(P.symbols[\"_flow_i\"].value > 0.5)[0]).plot(\n", " orphan_edges=False\n", ")" ] } ], "metadata": { "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.10.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }