diff --git a/.gitignore b/.gitignore index 2e1371e..96b44c4 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,5 @@ out.mp4 district_assignment.csv maup_concated.csv + +examples/.ipynb_checkpoints \ No newline at end of file diff --git a/examples/ensemble_analysis.ipynb b/examples/ensemble_analysis.ipynb new file mode 100644 index 0000000..e1d747f --- /dev/null +++ b/examples/ensemble_analysis.ipynb @@ -0,0 +1,290 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Ensemble Analysis for New Hampshire" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import networkx as nx\n", + "from gerrychain import (Partition, Graph, MarkovChain, proposals, updaters, constraints, accept)\n", + "from gerrychain.proposals import recom\n", + "from gerrychain.tree import recursive_tree_part, bipartition_tree\n", + "from gerrychain.random import random\n", + "from functools import partial\n", + "import pandas as pd\n", + "\n", + "from rba import constants\n", + "from rba.util import (get_num_vra_districts, get_gerrymandering_score,\n", + " get_district_gerrymandering_scores, get_county_weighted_random_spanning_tree)\n", + "from rba.visualization import visualize_partition_geopandas" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "# CONSTANTS\n", + "\n", + "random.seed(2023)\n", + "GEODATA_FILE = \"../rba/data/2010/new_hampshire_geodata_merged.json\"\n", + "COMMUNITY_OUTPUT_FILE = \"../rba/data/2010/new_hampshire_communities.json\"\n", + "VRA_CONFIG_FILE = \"vra_nh.json\"\n", + "NUM_DISTRICTS = 2" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "305\n" + ] + } + ], + "source": [ + "# LOADING DATA\n", + "\n", + "with open(GEODATA_FILE, \"r\") as f:\n", + " data = json.load(f)\n", + "nx_graph = nx.readwrite.json_graph.adjacency_graph(data)\n", + "graph = Graph.from_networkx(nx_graph)\n", + "del nx_graph\n", + "\n", + "with open(COMMUNITY_OUTPUT_FILE, \"r\") as f:\n", + " community_data = json.load(f)\n", + "\n", + "edge_lifetimes = {}\n", + "for edge, lifetime in community_data[\"edge_lifetimes\"].items():\n", + " u = edge.split(\",\")[0][2:-1]\n", + " v = edge.split(\",\")[1][2:-2]\n", + " edge_lifetimes[frozenset((u, v))] = lifetime\n", + "\n", + "with open(VRA_CONFIG_FILE, \"r\") as f:\n", + " vra_config = json.load(f)\n", + "\n", + "vra_threshold = vra_config[\"opportunity_threshold\"]\n", + "num_combined_vra_districts = vra_config[\"combined\"]\n", + "del vra_config[\"opportunity_threshold\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [], + "source": [ + "# UPDATERS\n", + "\n", + "rba_updaters = {\n", + " \"population\": updaters.Tally(\"total_pop\", alias=\"population\"),\n", + " # \"gerry_score\": partial(get_gerrymandering_score, edge_lifetimes=edge_lifetimes),\n", + " # \"district_gerry_scores\": partial(get_district_gerrymandering_scores, edge_lifetimes=edge_lifetimes)\n", + " \"gerry_score\": lambda p: 0.5,\n", + " \"district_gerry_scores\": lambda p: [0.5] * NUM_DISTRICTS\n", + "}\n", + "\n", + "vra_updaters = {f\"num_{minority}_vra_districts\": partial(get_num_vra_districts,\n", + " label=f\"total_{minority}\",\n", + " threshold=vra_threshold)\n", + " for minority in vra_config.keys()}\n", + "\n", + "rba_updaters.update(vra_updaters)" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAGvCAYAAACQBNO0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACgGUlEQVR4nOydd3gc1bn/vzOzRXXVq61mddlyt2VhbJJrBwHOvTiQhOIkJjhw4QIJmNBCQiAFEpPcAIHgHze5kAIJkJsQQnFi7BiDMS5ylWRVq1q9rtruTjm/P87sasvM7qwsyZI9n+fRA9o9U1bed845b/m+DCGEQEdH55KBvdA3oKOjM7PoRq+jc4mhG72OziWGbvQ6OpcYutHr6Fxi6Eavo3OJoRu9js4lhm70OjqXGIYLfQOzGUmS0N7ejsjISDAMc6FvR0fHL4QQDA8PIzU1FSyrPp/rRu+H9vZ2pKWlXejb0NEJitbWVsyfP1/1fd3o/RAZGQmA/hEtFssFvhsdHf9YrVakpaW5vrdq6EbvB+eS3mKx6EavM2cItBXVHXk6OpcYutHr6Fxi6Eavo3OJoRu9js4lhm70OjqXGLrR6+hcYuhGr6NziaEbvY7OJYZu9Do6lxi60evoXGLoRq+jc4mhG72OziWGbvQ6OpcYutFfQhD7hyBC04W+DZ0LjG70lwiEOEAG7wcZ/9uFvhWdC4xu9JcKYidArIDtnQt9JzoXGN3oLxUYE8DGA2IziOPkhb4bnQuIbvRzECINTeIgHiDj9H/12f6SRjf6OYZkfRyk/xZI/V8DIaL2A8UWgIzS/7e9F9yxOhcVutHPIQgRgLE/A0Il4PgUsL0fxNH8xP9KPfR4nUsS3ejnCIRIIIN3AXBMvDb6IghfqW3W5qs8z2f75xTfoc5cQTf6OQIZfhqw/8vzRaEOpO96kP6vgYjddBwZn/h/sQfEtheE8CD8Cc9j7XtACJmBO9eZbegS2HMAYj8A2P6u8q4E8EdArN8HTCtBbLsAoRGI/A7I+GsAf4p67aVer8O66XumJdN+/zqzC93oZzlk/B2QoYfgsSdXwr4HxL5n4jjrdwBI9Bdvg3eOse8Foxv9JYe+vJ/thGwCDPmTOFAKPMT+wSTOqzPX0Y1+1sMj4Cw/WYQ6PRf/EkQ3+lkO6d8KCDXTdwHbruk7t86sJCijF0UR3/ve95CVlYXQ0FBkZ2fjhz/8oYcX+JZbbgHDMB4/V111lcd5+vv7sWXLFlgsFkRHR2Pbtm0YGRnxGHPq1CmsW7cOISEhSEtLw44dO3zu580330RBQQFCQkJQXFyM9957z+N9Qggee+wxpKSkIDQ0FBs3bkRdXV0wH/mCQKQBSL3Xgdj+BQjV03utoGL9OhcDQRn9T3/6U7z44ot4/vnncebMGfz0pz/Fjh078Mtf/tJj3FVXXYWOjg7Xzx//+EeP97ds2YLKykrs3r0b77zzDvbv34/bb7/d9b7VasWVV16JjIwMlJeX4+mnn8bjjz+Ol156yTXmk08+wU033YRt27bh+PHj2Lx5MzZv3oyKigrXmB07duC5557Dzp07cejQIYSHh6OsrAw2my2oP9JMQWx7QEZfBrH+GBAqQAb/cyKLbroQzoAILdN7DZ3ZBQmCTZs2kVtvvdXjteuuu45s2bLF9fvWrVvJtddeq3qOqqoqAoAcOXLE9dr7779PGIYh586dI4QQ8qtf/YrExMQQu93uGvPQQw+R/Px81+9f/vKXyaZNmzzOXVJSQv7zP/+TEEKIJEkkOTmZPP300673BwcHidlsJn/84x81fd6hoSECgAwNDWkaf75IkkSk8feJ2PN5InbkztiPNPz/ZuTz6UwvWr+vQc30l112Gfbs2YPa2loAwMmTJ/Hxxx/j6quv9hi3b98+JCYmIj8/H3feeSf6+vpc7x08eBDR0dFYuXKl67WNGzeCZVkcOnTINWb9+vUwmUyuMWVlZaipqcHAwIBrzMaNGz2uW1ZWhoMHDwIAGhsb0dnZ6TEmKioKJSUlrjHe2O12WK1Wj5+ZhGEYMCFXgYn7G5ion07Sax88xK5n511KBBWnf/jhh2G1WlFQUACO4yCKIn784x9jy5YtrjFXXXUVrrvuOmRlZaGhoQHf+c53cPXVV+PgwYPgOA6dnZ1ITEz0vAmDAbGxsejs7AQAdHZ2Iisry2NMUlKS672YmBh0dna6XnMf434O9+OUxnjz1FNP4YknngjmTzItMAwLhH4BMOSA9F0//RfkT4GI7WC41Om/ls4FJyijf+ONN/Dqq6/itddew8KFC3HixAnce++9SE1NxdatWwEAN954o2t8cXExFi9ejOzsbOzbtw8bNmyY2rufYh555BFs377d9bvVakVaWtqFuyGxfeauZdsNhG+duevpXDCCWt4/8MADePjhh3HjjTeiuLgYX/3qV3HffffhqaeeUj1mwYIFiI+PR319PQAgOTkZ3d3dHmMEQUB/fz+Sk5NdY7q6ujzGOH8PNMb9fffjlMZ4YzabYbFYPH4uKNJI4DFThF6Ac+kQlNGPjY2BZT0P4TgOkqSe/dXW1oa+vj6kpKQAAEpLSzE4OIjy8nLXmL1790KSJJSUlLjG7N+/Hzw/kZSye/du5OfnIyYmxjVmz549cGf37t0oLS0FAGRlZSE5OdljjNVqxaFDh1xjZjvE8dHMXYwvBxGV03V1LjKC8Q5u3bqVzJs3j7zzzjuksbGR/OUvfyHx8fHkwQcfJIQQMjw8TL797W+TgwcPksbGRvLBBx+Q5cuXk9zcXGKz2Vznueqqq8iyZcvIoUOHyMcff0xyc3PJTTfd5Hp/cHCQJCUlka9+9aukoqKC/OlPfyJhYWHk//2/CS/zgQMHiMFgID/72c/ImTNnyPe//31iNBrJ6dOnXWN+8pOfkOjoaPK3v/2NnDp1ilx77bUkKyuLjI+Pa/q8M+2990Yc/N7MevFHtUU1dGYnWr+vQRm91Wol3/rWt0h6ejoJCQkhCxYsII8++qgrtDY2NkauvPJKkpCQQIxGI8nIyCC33XYb6ezs9DhPX18fuemmm0hERASxWCzk61//OhkeHvYYc/LkSXL55ZcTs9lM5s2bR37yk5/43M8bb7xB8vLyiMlkIgsXLiTvvvuux/uSJJHvfe97JCkpiZjNZrJhwwZSU1Oj+fNecKPv3TKjRi/2ff2CfE6dqUHr95UhRC+qVsNqtSIqKgpDQ0MXZH8vDdwJ2PcEHjhlGMAkHgTDRs3gNXWmCq3fVz33fjZDHIHHTCkCYN87w9fUmWl0o5+lECICjiMzf13di3/Roxv9LIVhODCW78z8he0fg0hjM39dnRlDN/pZDBN2I8Aq5xRMH3bAsX+Gr6kzk+hGP4shQj0g9QUeONXXte2e8WvqzBy60c9SCCEgwzswbao5/rD/C2TGnYg6M4Vu9LMUhmEAofnCXJyMAI5PLsy1daYd3ehnKcRxHBAbL9z1dS/+RYtu9LMUcqGVau179X53Fym60c9CiNgOOI5ewDswAWwyXW3oXHTozS5mIWR4B8BPwuDYOIBbIP/CuL0hyq2qHQDGAWkcMKQDvNuDhZ0HcKkAGQOEekCoAuz/BMwroXNxoRv9LIMIjSqzPAMYcun/im3UOLlcgI0CIFGjBg/wNQC0yHxlAcblAAggDQHiWUA65znEvheEPEKdijoXDbrRzzbYBNpnzhsuGxBq5V8sdGYWmwBxkiE9YgeEE/7HiC2AUAcY8yZ3DZ1Zib6nn3VIABPt+RJjAST32dsqz8rnEcMn2lJtid4M46JDN/rZBhkBzOsBcBOvGRYARGH2Px+kAW3j9GYYFx260c8yGC4VbPTPADaavsCmAt695acC0guPB4saYgMIXxt4nM6cQTf6WQiRxiZy7tnpEu8gABurbaQ+219U6EY/G2FCafgNAJiwabyORoUc2y7oAksXD7r3fhbCMAwQ+ShgyAMZeW4KzsgCCAeYEIAxAYyZ/pcI2g4XG3Qv/kWEbvSzEGI/COI4DIABpF7AsARgOFDjZeQfAkCSfwhAJACCnIDD05CcIR1wnAJgAzAMkGF62GTuyfY+GN3oLwp0o59lEPt+kMF7qRd//E/ndzIpHtTgpwDb+yAR39QTdS4C9D39LIOM/Yka/FTAaPDOa0U865YcpDOX0Y1+tkHsU3gy4xSeS0/UuVjQjX6WwUT/nKbcTglT/M+rh+4uCnSjn2UwbDSY2N8DXFbgwQGZ4np48SwIXze159SZcXSjn4UwXDyY2N+e/4mmQ+fO/o+pP6fOjKJ772crbDRN0iHjXm+EAEwkwIYDMAGMkcbcYYBnDb1E4/FTDLG9Dybi7ik/r87MoRv9LIVhQkDM1wD8YYCIALEBZBSAjf6/2KPhJIlTf2NCHYjQBMaQOfXn1pkR9OX9LIYJWQ+IrYDUDpB+AEF69kn/tNyX7tCb2+hGP5sxbwDC/+s8TiDQWvwpRlfKndvoRj+LYRgT2Mh7wUTcD8/9ehBMR9tpoRJEaJr68+rMCLrRzwXM6wCYJnlwyFTeCYXLBbHr/e7mKrrRzwHIyPMIej8PAMYVgDhFcXVDIT0fG0fPOf5/U3NenRlH997PBcSOSRxkAfhj53FRC2DIAhgDFeAUzni+LZwBEVrAGNLP4xo6FwLd6Gc5ktBJC3AM+bKghvOfzCsmDwKqby/IPw4vMU0NcFkTarxiMyCc9D/e9g8g4rbgrqFzwdGNfpbD2P4KIk6ykSUTBZAhP++HAYYcOk6oob3zguifR+z/AKMb/ZxDN/pZDBE7QWzvTP4EbBQgehk9EwkY8miCj1AD8KcARGBSdff8aRCxAwyXMvl71JlxdEfeLIaM7KQyVZOFiZj4r3EFYFhE03r5ckCoBOCUyxqhD4Lg7xCw7Z78/elcEHSjn8UwhrTzPIEcruMyZEOvwIShe4+NmNQl9ESduYdu9LMYorUhhSqycg4TGnioU3I7WPijINI0pfvqTAu60c9imPDbp+ZEYpuGMQ0Ao00H3xMJsO2bxHE6F4qgjF4URXzve99DVlYWQkNDkZ2djR/+8IcemuiEEDz22GNISUlBaGgoNm7ciLo6z31pf38/tmzZAovFgujoaGzbtg0jI566cKdOncK6desQEhKCtLQ07Nixw+d+3nzzTRQUFCAkJATFxcV47733PN7Xci+zGvuH53kC+d9F6gRYDVsFw+SEO4hdX+LPJYIy+p/+9Kd48cUX8fzzz+PMmTP46U9/ih07duCXv/yla8yOHTvw3HPPYefOnTh06BDCw8NRVlYGm23CO7xlyxZUVlZi9+7deOedd7B//37cfvvErGa1WnHllVciIyMD5eXlePrpp/H444/jpZdeco355JNPcNNNN2Hbtm04fvw4Nm/ejM2bN6OioiKoe5mtECKBOA5M7mCuADAupvt412vJGi6qUQffG/sBEGl0csfqzDgMCaJ1yec//3kkJSXhN7/5jeu166+/HqGhofjDH/4AQghSU1Nx//3349vf/jYAYGhoCElJSXjllVdw44034syZMygqKsKRI0ewcuVKAMCuXbtwzTXXoK2tDampqXjxxRfx6KOPorOzEyYTzTl/+OGH8dZbb6G6uhoAcMMNN2B0dBTvvDMR0lqzZg2WLl2KnTt3arqXQFitVkRFRWFoaAgWy3S1l1KG8FUgfZuDOMIEGIupTr5SXN9QAAjVAc4RBtoJN/huuEzUM2BCrwn6OJ2pQ+v3NaiZ/rLLLsOePXtQW0ulkE+ePImPP/4YV199NQCgsbERnZ2d2Lhxo+uYqKgolJSU4ODBgwCAgwcPIjo62mXwALBx40awLItDhw65xqxfv95l8ABQVlaGmpoaDAwMuMa4X8c5xnkdLffijd1uh9Vq9fi5YBgKgfA7A49jEwDjOuqp58uVDR6QQ3+RAU42NsnQnb7En0sElZzz8MMPw2q1oqCgABzHQRRF/PjHP8aWLVsAAJ2dnQCApKQkj+OSkpJc73V2diIx0VPRxWAwIDY21mNMVlaWzzmc78XExKCzszPgdQLdizdPPfUUnnjiCQ1/iemHYRiAy1BvSMNlAWwMTa6RDtJEHL9rNhEwZgfugKvF0+8NmwxIgyDEDmYaJLp0ppagZvo33ngDr776Kl577TUcO3YMv/3tb/Gzn/0Mv/3tFIg4zgIeeeQRDA0NuX5aW1sv7A0ZFni9wAKGxXTPLjbKBTUC/eG8xyqhoSZf0PCZmXDqMzCuoJJcUifg+ASwT9IHoTOjBDXTP/DAA3j44Ydd++Hi4mI0NzfjqaeewtatW5GcTJ1FXV1dSEmZSM3s6urC0qVLAQDJycno7u72OK8gCOjv73cdn5ycjK6uLo8xzt8DjXF/P9C9eGM2m2E2z46ZikiDIKO/kn8LAYyLaLWdcEr5AL4WtObejwKuoCGvnnQBXCatrHOHS5Nn9GFaWsv73gexfwAm5N8CX0PnghLUTD82NgaW9TyE4zhIkgQAyMrKQnJyMvbs2eN632q14tChQygtLQUAlJaWYnBwEOXlE57lvXv3QpIklJSUuMbs378fPD/hUNq9ezfy8/MRExPjGuN+HecY53W03MusZvyv1EiNqwAYAP4oIJ3zc8AQdeT5gwwCXE7ga7OJtBjHuAQwltDZXGwF+COAWA1VPX3bHhAyxVr7OlNOUEb/7//+7/jxj3+Md999F01NTfjrX/+K//7v/8YXvvAFAHQfeu+99+JHP/oR3n77bZw+fRpf+9rXkJqais2bNwMACgsLcdVVV+G2227D4cOHceDAAdx999248cYbkZqaCgC4+eabYTKZsG3bNlRWVuL111/Hs88+i+3bt7vu5Vvf+hZ27dqFn//856iursbjjz+Oo0eP4u6779Z8L7MZwldRpxx/BIDG3nZaSmlZPwk4bBJgXAlII7Q0lz9JFXhJt/ox7pCB86zh15kRSBBYrVbyrW99i6Snp5OQkBCyYMEC8uijjxK73e4aI0kS+d73vkeSkpKI2WwmGzZsIDU1NR7n6evrIzfddBOJiIggFouFfP3rXyfDw8MeY06ePEkuv/xyYjabybx588hPfvITn/t54403SF5eHjGZTGThwoXk3Xff9Xhfy734Y2hoiAAgQ0NDmo+ZCiRHNRE7CojYkRv8T/dVAd4v8/r9s0Ts3aJ+XO/NwV1/6KkZ/VvpTKD1+xpUnP5S40LF6cnor0GGfTMQNWFcBvDH/Y8xLKYhPrEDkAI47rh8QKzRfn0uA0z8P/WW1heAaYnT60w/xHEMZPiZyZ+APx04h57haBONQAYPyKIaQfh7xWaax68za9GNfpbBmJYDodedxxkEwBCo620ws7Aj+Jx8vcZ+VqMb/SyEMZUAOI/QodAAv7OzcDbIG4oOajix7wk8SOeCoRv9bCRkIxjLdzDpfx7S7z98pzV050TqDe76/CkQUTnrUefCoxv9LIRhzGDCbgITce/kTxIofOcvdOeN2Ehz/IPBvi+48Tozhm70sxgm4g4g/K7JHSw2+BbPMKET6bMkyFJYLiOo4cT+r+DOrzNj6Gq4sxwm/Osg429SLfqgDw6j2XVcOhXEFGon0meDbWxJ1DQIGPpwYSwA+ImCHvtBEGID49Tp05k16DP9LIdhLYD5s5M7mD9Ne9vzR2X1W7c6eWKlhTtaEWox0RcvRE7RXQYgUpbSPgLadMOJDXB8Orn71plWdKOfAzChXwBMQdYLMEmAcXmAZbyf4hyf80UCplVURhsiTdHljwNw8x14+RGIbW8wd6wzQ+jL+zkAY1oOhN8O4lAW//CAy6BNJvlTAN/lf6x4FmBTAEmlVx47D+BSAWkQEOsBx0f+z+e9BbHvAyFEz86bZehGP1cwFgGGXJo6S9wKcAwL5f8Jpc0m+U/V1XOU4NLcjN5Ar8FEAlIP9dr7rexzR6GFltRJl/7GILYROtOObvRzBIaNAWL/BAi1IMPPARgHpCF5ry5jXBb8iaU+wLAcYERAqPftTqsVLt63hRZAQ3e60c8q9D39HIJhI8GYVgDmK6iX3LvZJF+PoJ/jYgNdujvLaCcLq6y/p4fuZh+60c9BmPCvA+HfgG8O/bCGvHsFyPgU3JVJ+WX+uJ6dN8vQjX4OwjAMGDYBivn5wcbfwWIykte++NHM1/vdzSp0o5+rmNYAsPu+HnQSzxQlz/jZGhD7B1NzDZ0pQTf6uQobR8Nt3ojNVMBSK8ZFU3M//ppYOo6CSApOPp0Lgm70cxU2klbimT/j+x6nscW1YYm25paKWGi0wLgcYGIAvx12hSnoy6czVehGP0dhmFAwIVeCiXrGdx9PxgKfwFAICFUAY9R+US6DqvMa8gGM0ow8/hgt1fW3pwdAbPoSf7agx+nnMITYQEZfpHn07gg1tNhGzfi5bEBsAcADfjvSMFQjj0uiDwixWSXxxwwgQFNQ+z4QaQwMG+Z/nM60o8/0cxpOxbAFmlmnBDufimK4HG/eobZIwLhUXrZHUZ17x0cAm+rnPlTCdR7YAqfx6swI+kw/ByHSIK16I2PA2J9URin80zKJAOxe6bIG2hKLjafLdKFBod+dRNNywWKiks6tmw5jCtBHT75v2wdgQsoCD9SZVnSjn4OQ4R3A+J/9DxKaQJN3nNbI0tx8qd1zHGOimnliAN08qR0wltK0XWlAPtcQwIbS7QKvQVLLvg+ECGAY/Wt3IdGX93MQxvxvNGTnD9InO92Wgj7bJSiX0gbR9oA/CICnWwepCzAupAIdksaMOzIEOI5qv57OtKAb/RyECdkIJmE/EPk9ecmuAn8YEE4Ahhz6AFAUuAyy14nYCPCf0LAgsQPiALTt6eWr2XV57AuNbvRzFIYx0v0x6aazvsv5FqMwOFJWtlFikg2OxGYaIZCa6FI/IGbAUAwIk80L0Jkq9M3VXMZxmP5X6qM/TrhMql5Lxmgc3m9TScnPe4FwxuZVvkZssrwiGKdhROE0ANqckzEWncd1dc4H3ejnKETsBhn9f8pvik2+/eVVTzTZYpsQGuM3rpRj/gDA0io/JprWAIjNyvt9+x4qCqJzQdCX93MUhksEE/0r9Xi8ZhSKdtRgE6mRGxYBkADHASq6KXVTnwEiAKFO7mOvrt5DbHoHnAuJPtPPZYgVYCLO8xz+MukYOYYfR731jIUauSISPEQy/SFUgYgdYDiFgiGdaUef6ecy9v2B21IHwiejz0hncuNK2v1WbKBRALE5QMquGNx19dn+gqEb/RyFkHEQoT64g9hEugznFky8Jo1S775xuRzTNwJCBZ3RSZ/n8UILVBHqAHCab0Vvcnnh0I1+jkGIHcS2C6RrKWD7u7aDuExaBiv10f02GzXxnnE5ne35YzSmDz8VeqSbJuMovjcanH/BcRhEGgk8TmfK0Y1+rjH2KiA0QlN83bAIMBRRTz5/HK4lOH96Ip7PHwYQqv36bJL6e4yyOKYyvF6Ac4HQjX7OYQQZed7/EEMRjY8LFbQk1geBZuk5/98YxAztrxuuFKC5hhfEtiuo8TpTg270cw02En4Tagx5gNgDiK3+zyM0YCJ4E0RWnlgHMOEq77XQdlpase8HIUG01tKZEnSjn2OQ4Weg6innMgChAzBkaThRP8DNk1N3Q6DdCSf5l9k2aJTqAqgfQEurLp0pRTf6uYYhU/l1NgmQRgAMy7O8hv5xYjN14Dk+BQzBdKHxI7GlSUNfzsM3rgCx7w/iujpTgW70cwgijQEIg88/GxMNwDARYpM66L4+GJggZKyEJj/v1YLeozusrK+3AjAuBsDQPHy+HLD9A4ScT/6/TrDoGXlzCDJwK52ZDYsBoRqAA2BCATZWQQRDe8wcgFv+vJYb6aNhQMX8fp7KahM73fuTEbkWQEVfT+qmzsapkuLWCYg+088RiDRMFW4AQDgFcFl0hueylFVvhApaaacVqYueSyv+zs0fp4bMH6ENMQMt+fV+dzOKbvRzBGJ9QpaalhFraKxdtU+cRGfjYGDj1d/jZCEO4wr59Brz7DVAbLrRzyRBGX1mZibto+b1c9dddwEAPvOZz/i8d8cdd3ico6WlBZs2bUJYWBgSExPxwAMPQBA8NdP37duH5cuXw2w2IycnB6+88orPvbzwwgvIzMxESEgISkpKcPjwYY/3bTYb7rrrLsTFxSEiIgLXX389urqCiyPPFsj43wHb275vSI0AGwYwKtJZQiOC+if2UNYJpdsI4wr54VJPZ27+OH1dCmI7EAihAkRsDzxOZ0oIyuiPHDmCjo4O18/u3VT66Etf+pJrzG233eYxZseOHa73RFHEpk2b4HA48Mknn+C3v/0tXnnlFTz22GOuMY2Njdi0aRM++9nP4sSJE7j33nvxjW98A//4xz9cY15//XVs374d3//+93Hs2DEsWbIEZWVl6O6e6ON233334e9//zvefPNNfPjhh2hvb8d1110X/F/oAkJGfwOpcyHI0P3qg8Q2qkuveIJeuQw2EFGyFn4jYLqMetYh0m0EXw4Q9+41En19SjrdOjGC2PXsvJmCIYRMUi8JuPfee/HOO++grq4ODMPgM5/5DJYuXYpnnnlGcfz777+Pz3/+82hvb0dSEv2i7ty5Ew899BB6enpgMpnw0EMP4d1330VFRYXruBtvvBGDg4PYtYtmcJWUlGDVqlV4/nmamSZJEtLS0nDPPffg4YcfxtDQEBISEvDaa6/hi1/8IgCguroahYWFOHjwINasWaPp81mtVkRFRWFoaAgWS7DdYCcHISLI8M8ACDSUJtQEPsi4TL3ajiug2vVKsAl0H8+fArh4+gAxrvIjrTXVMFTmi68BTEvBxr4yQ9e9ONH6fZ30nt7hcOAPf/gDbr31VjDMREz41VdfRXx8PBYtWoRHHnkEY2MTBRwHDx5EcXGxy+ABoKysDFarFZWVla4xGzdu9LhWWVkZDh486LpueXm5xxiWZbFx40bXmPLycvA87zGmoKAA6enprjFK2O12WK1Wj5+ZhDjKQYYeAsZ+A/BntBm8YbH/8lqxxbe3HZtKl+3SAKhCrm3C4L275UwrRFbuGQMcR0Ak9c63OlPHpEN2b731FgYHB3HLLbe4Xrv55puRkZGB1NRUnDp1Cg899BBqamrwl7/8BQDQ2dnpYfAAXL93dnb6HWO1WjE+Po6BgQGIoqg4prq62nUOk8mE6OhonzHO6yjx1FNP4YknntD+R5hCCHGADN4zsa/mDwPG1XJBjBom9dbUTDRNyRWqaeKO2Er/y82nDwmn/r17GuyMzfDu9+lU0pULcEKumvl7uMSYtNH/5je/wdVXX43U1Il2R7fffrvr/4uLi5GSkoINGzagoaEB2dl+UjdnCY888gi2b9/u+t1qtSItLYi00vOFS/d0pvGH/S+3jUt832PiaRouf3rigcGfkbcAp3yLYoQqKmCpVbt+qnGL3RPbXjC60U87k1reNzc344MPPsA3vvENv+NKSkoAAPX1VOwhOTnZx4Pu/D05OdnvGIvFgtDQUMTHx4PjOMUx7udwOBwYHBxUHaOE2WyGxWLx+JkpyMizyqq1/BFZf84LxkL3wt5wsfKDwF0Ga9SztNbnGJUa+ZlA6gNY+cFq/1DPzpsBJmX0L7/8MhITE7Fp0ya/406cOAEASEmhWmilpaU4ffq0h5d99+7dsFgsKCoqco3Zs8dTVWX37t0oLS0FAJhMJqxYscJjjCRJ2LNnj2vMihUrYDQaPcbU1NSgpaXFNWa2wZjWq7+pZPiGfChq0gVV0y4jNOCCpmw4ow9kIIBct85UEPTyXpIkvPzyy9i6dSsMhonDGxoa8Nprr+Gaa65BXFwcTp06hfvuuw/r16/H4sWLAQBXXnklioqK8NWvfhU7duxAZ2cnvvvd7+Kuu+6C2Uz11+644w48//zzePDBB3Hrrbdi7969eOONN/Duu++6rrV9+3Zs3boVK1euxOrVq/HMM89gdHQUX//61wEAUVFR2LZtG7Zv347Y2FhYLBbcc889KC0t1ey5n3EMC+S0VRVnFn9kYo/Ppig0mXSOO0W7zXo0qQwA6ZMbUZwO+rYnDWOhDy4ySu/ZeSv2f4ExrZy5+7gECdroP/jgA7S0tODWW2/1eN1kMuGDDz5wGWBaWhquv/56fPe733WN4TgO77zzDu68806UlpYiPDwcW7duxQ9+8APXmKysLLz77ru477778Oyzz2L+/Pn49a9/jbKyiW6nN9xwA3p6evDYY4+hs7MTS5cuxa5duzyce7/4xS/Asiyuv/562O12lJWV4Ve/+lWwH3fmYCNoyaqbAfjAHwaMJQAEgEsFwNDlsdiKicYTPHXgBe2Um3TkVhtMHJ3RmXBa3OM4qHyPtg+AyAem914ucc4rTn+xM1NxeiKeA+n/Kg2bBULRo2+m+3I2GoBAPfJCFYIzZFaWuu4J4hglIgEuBWDlvxcZAcR2z1CgoZDm5KvAxL8HxqXso6MVrd9XvcruAkOIDWTwfm0Gz85T2fPaqaKNu5/OsIgW3WgmHGBjJmf0xqU03i510eiDOOxfEVtoAGCGaqMN224gQjf66UIvuJlBiOM4pN7/gDTwTZDR34HY9oL0fFaj84qjZbQQAo6cFMYCuRZ+MjCAUKnSFVcJh1/1HWLfN8n70NGCPtPPIMT6fZosI1SD2IMUhTQuD26fLlTSzDtJYyELf4QKXPjzKahBgmx0Aajr7AEAfxJEGgTDRgd/Xp2A6EY/jRChHnCcABl/i5bFTnomjaJGHNzVqQaet9GzaXQZL3ZQ6S2hnv6XiADfQFN2A4lq+jCJJphSv783AfvHQOjngz+vTkB0o58mJKEHGPw27TUn9WjvIquEcTLeeMjOslAA43JPOgudyUk/dfY5z8nLVXTGFbQ+31/oUAkyiS2H2AjAArX+d8S+D4xu9NOCbvTTACF2YPi7nprzTDRgLKbdZMR2qmOnBXbe5PvVkRG5iGacOvWcq3DVXICTVCyDjZFXJVrlqYPofOtCounCwkmVU+4HISIYJkjZL52A6EY/hRBpGLDvBRn7E61D93hzkOrGOV9n46hYJBhAOAcQldx3NhaQzk3+poJaIQj0YcSE0bAbEbRdm0zG6AEw/lR1B6mD06SQgqxzXuhGP5WQUZChR+DjYWcSAYbzzJKT+ugPQJ10goNWwDEhdCYWWgHDvJnNkgPoPXLptPpNGtZ4zCQFNQKo5RDbHjC60U85eshuCmG4ZJpO640hg86gavXxYgvdZwunaOKNUAnASrPYjKtpht1M/lMxIXRFItZq08P3aXetEandf0ccvbPttKDP9DOBWs07QGdVNflpqZsaHiD7BJbIsyoPSEM0GSYYh5tW+Aq3/P1Ae2oO55U7YJgP8CrahWIziNAAxl9HHZ2g0Y1+CiGjv/ENyzFxynrvTtgkFaNnPENnZJAavHcKLhNNZa/YCNB/ToEuy6U+L227YBinGX38Ebrq4BYoy2wDoI0tNG4DFAlQSmvf57+Nlk7Q6EY/RRChGWR4h+8bhgyA7/NzoMp+mE1U6AKrYCBkEBAHVdJeQwEugc7aTIj8moPKV0s91LuvBl8DIARUSqtzQl7L2/iZUIBMwuiZMHnbYgDduigbP7H9C0z4tuDPr6OKbvRTheMogq9UM9DkGCXYOF+jJ8EmwYwrryK4TNngLfJDIRxgzKD97xy0Jx4ZALh82fs/NhF14DI8Vy6M2fdju3e/MRQCQjMtBiJ2OmuTcZqZ6CwP5vKpjr8S/DEQaQgMGxXkZ9dRQzf6qYJh5G4zjZ6v+yukMeSqV5sp9paborx74lTVsQKiihAml0/j9oYl9FnAy/F0JhqAu9HLGndcFt1miO2ywYdRhV2hEYANkGRnn9Kqh7X4KdARAPuHQOh/aP54Ov7RvfdTBBN6HRDiqeILQ5F/J55flRuFrrOTjYf7nEeDt11sARgDTZ7hTwKGhfR14aScX+Ac1wMw8+nDjj8MSM6H3Ji8yrB5n9kXn22M1+3qba+mFN3opxAm7CtUEcb1QoCFlL/4tpKBEw0GpAVNIbZxz8637qm2TLTbOCvAnKckg9hCxTmV4OYD/GmQoLc2OmroRj+FMFwKmJhfYWLXFEbTYA0LQZ1i7lj8F+AoFaRMSXguBJq3CfxRmicA0D23IZ/+v/dsPxXVcNx8+d6MNDfAuEoOZ7YBTCSIozzQGXQ0ou/ppxjGtBqEjaHed/5Tt3cM1KnFRFIvOBtFjUoRTlmSeiqMngkPbsXAH3ZztLnNEWzihEOP8X6gBQmXSVdFhhzq2BS8OvIIlYD9H4B5luobzjF0o58WRFq6ChMmilYET6cds4x66NkU2XPO0eOIHTQfvwWAu7ilAdoLYPzAhFEhzGCQegCE0PtnE+jvHisRBf+DPwx5NIwo9QBit9y/vsn/MbYPQCIf8+impDM5dKOfDsxX0l7y/rrTCKdp0otfSSsDzcuX+uXkGxMAHhB75UKYSWjET2ZWNiygzTMAugyXemi8nomQQ39BPoyYqOBLhaUu+rcyFgd3nI4PutFPB0KthnJYge6F/YrOCAActE20zzgTFclgYmm4UBqQl9sB9uuuNlIaMSyWl9tOx6LzQUPodgVQl+NWQ2yEa2UTBMT2LzC60Z83utFPB4ZMbbp34lkE/PKrdnxx0H0/cc8LMAFc7oQqrtjpW7fvr5zVG8MS33p3aWTinvkjAFjaSov4CU16ECLH9NMBIcjGFo4PAXwzuGN0fNCNfjrgqwKPAaiQpGEpIJzwN0j9LSbCK+zn8FXFZSxUAsvlwAtiplcKOYoNXv31JIAN0T5ps7F0ZeDu/dcKXwEidoPhEoM/VseFbvTTABNyNciISk94HwIl3PiJgTOhgU9PrJPQ1wNV7PH2ojvhj9NcfLD0/sQgml9K7TSEKVT6pvQGhAD2vUDYjUEco+ONHqefBpiIOwHzFdoGC2cCzHr+ZnoNRj8Z2HlylZ5aiFCgufj8ERp2lPykGivhTDFm5wV9a8S+N+hjdDzRjX6aYCIf8czO8wcb7+dNf0Z/nvFxNcgoNKXPTgYmlvoVDIvog4ONC+54+6cgU5WZeImiG/00wRgWgIl7nSaeBII/7ZXa6obf1s3TtDvjEqbuXGwqDTsaV9D/J/2A4xM5VGkHuGBr5W2A/eDU3d8liL6nn0YYQzYIlwqIA4AxC4BBjnG3wnMGd1CHnmJc35/RT9cze7LnZelDjo2jZcBiC93D+2u4wVfQFRFRqfZTgNj3gAn57CTvUUc3+mmGMX8GxPGJVyw7EjCkyx71ETkjTS12PYnuMecLEwGaZaehkIYJncjJ5+toGFJVZUeJMcCg1JTTD/a9IEQCw+gL1cmg/9Wmm9DrANNarxeHqfeaP0y18YmdLu9N6+TKNvew2gUwer488LaES6dxfCLSBxp/AsAkawOEs/A//3i9J/VOrv2WDgDd6GcAQYOclEjj62RQbpDB0LCWcVWAmvtJpOFqggksa03sVL13KuoBSC/to+eNcTnN03euJDwuryvlThZ9eT/NkNHfTWJWsnvG1tlUmnJLxuXec9my5z4M/lpDTR5D4AeV1AUYiqdOl5+IgHElbbdFhuT4vUEuPzYACIfHSsK+D4i8f2qufYmhG/00w7DR2pXz1Dz1Hs4wE3V88U4PNks94GwcNVShDur592FUmoo4qBddEZZuMdTaTTnvwbgEwWsCKsDlALB5Xs+4ihq9NCi/INBCG/d6BqEGRGgFY0g7/3u4xNCX99NN6HWAaR2YmFc0DNayXHfQ4hr3Y8QG2T9wZiIGblzlm/RjLJDr9EU3oYpcrzEK+fYuTHQ2ZixyYk65uuKNVliLr46gUzVIrAOYGOeL8n85gMsDjKupg1QnaPSZfpphWAsQ82uNpaRanXZ+9tHOZpVO3HvmwSCr7PYBglutviGPrjKkYZUOtGa65xYavIQ/CHXoKQl+aCJSuU5BcCtEMiyQs/+qAeNSur0Ra+mfijEBYTdM8tqXLvpMPwMwDEN16RQcUh4QjcvlYAQypT5a8UeGAf4QNVJvhFq6WmBDvPboJtmZGE4fWkpbAlel4CQw5kM5828EMC4DTJe55fXbaITAXavfcRjEqbKroxnd6GcAIllBhv9bvX2VC43adWQyHnO5pJY/TqW3fU8qO88Y2WteBMCsbuxOpF66nZgMzgaeSvBHaVWg3665juDi+zoAdKOfERjWQnu+B+zuqlHxdTK55+519ESAosSVYbGsPntMDh1q7VwzCT1+Lte3R4A3/DG5C446xP5x8Ne+xNGNfgYgRAJ4laYWHgO1Gs9klrRuS3Cxkc7mTgyF1AiFU57987QiVFGtPy0wSXJobhiaXEqKTT/c0KvugkZ35M0Eo7+iiTcB0ZpsI8JTdFMLXvtuoWFC117QKPqhCqFCHd4qPe6w8wAuieYsOLvUGpcFlhXjG+E3JViygvC1YIz+VwQ6E+gz/TRDiAjCaxSxCKahQ9BltV5GQwZpjv15G7yMUAePOcS4gnbs5TJkcc9OWULMbTXjisP7Y0g9Jdi4HAAH6B1wgiIoo8/MzATDMD4/d911FwDAZrPhrrvuQlxcHCIiInD99dejq8uzZVFLSws2bdqEsLAwJCYm4oEHHoAgeC5r9+3bh+XLl8NsNiMnJwevvPKKz7288MILyMzMREhICEpKSnD4sKdDR8u9zAhkFHBoVH4NxkEXjNYdAMVVBH9MbjIxBZABwOjm0CM2Ku0ttsjGrhCOFBsnxDX9wSbA5YgE6KrBUChHJQb0lNwgCcrojxw5go6ODtfP7t27AQBf+tKXAAD33Xcf/v73v+PNN9/Ehx9+iPb2dlx33XWu40VRxKZNm+BwOPDJJ5/gt7/9LV555RU89thjrjGNjY3YtGkTPvvZz+LEiRO499578Y1vfAP/+Mc/XGNef/11bN++Hd///vdx7NgxLFmyBGVlZejunhBnDHQvMwXDWoCQf4eH44yJowkmPs60YBxiLAADlZPWgqK/QFCv4w+KSNoJhxAaS+eyaBqx1I7AWXsaHl78YXpemOl1pB7PHgL8SRBRqzCnDkOI1uCwL/feey/eeecd1NXVwWq1IiEhAa+99hq++MUvAgCqq6tRWFiIgwcPYs2aNXj//ffx+c9/Hu3t7UhKSgIA7Ny5Ew899BB6enpgMpnw0EMP4d1330VFxUSCyY033ojBwUHs2rULAFBSUoJVq1bh+eefBwBIkoS0tDTcc889ePjhhzE0NBTwXrRgtVoRFRWFoaEhWCwaVXC8II7jINbvUq16QxYACZCGQLvYdNOuLkQEhBr6Hhw0442NkltIh4IumxnZcO109SAOgC590zWEAiG3jFZxJhoWy8UzQcLEyl1pqibi51wOlezWDEuz+lRr7hkAHN0uiM2qiUCM5SkwYdcHc/cXHVq/r5Pe0zscDvzhD3/ArbfeCoZhUF5eDp7nsXHjROfWgoICpKen4+BBmid+8OBBFBcXuwweAMrKymC1WlFZWeka434O5xjnORwOB8rLyz3GsCyLjRs3usZouRcl7HY7rFarx8/5QAgPMvLfdL9LBuhylD9BPeRkiApH8Mcm0l65TNBiFysdI1RTRxd/RE6zPUZnULEJru43RGMWn7+WWFIPglLJZVNkD/wovS/3hBk2Qvt56MVpMZEiEXIfQJYmFvnJ/COO/UFe99Jl0kb/1ltvYXBwELfccgsAoLOzEyaTCdHR0R7jkpKS0NnZ6RrjbvDO953v+RtjtVoxPj6O3t5eiKKoOMb9HIHuRYmnnnoKUVFRrp+0tPMs5hBqAb5G6Q3Z0Nyx0Tz0oGPeGj3+kj+j75CXzwHgsqnHXeqW03EVMgP50zQsFwz8adBqQTfYFICLkVOKNfg67AdpaFQnIJM2+t/85je4+uqrkZqaOpX3c0F55JFHMDQ05PppbZ1EzNodLgVBhdUm039eEmlGHBtA144M+X+fr1AR6OToA8FQKBf2HIf/GgGRqgIFhY0WA7kumUcTmYLJGSCDAO+vMlDHyaTi9M3Nzfjggw/wl7/8xfVacnIyHA4HBgcHPWbYrq4uJCcnu8Z4e9mdHnX3Md5e9q6uLlgsFoSGhoLjOHAcpzjG/RyB7kUJs9kMs9ms8a+ggbE/auwFLyPUATAjsBa+G4zoVmBjAQzz5IQWAzUEsRUA67kEV75Z2p1W6pV/DweMC+lWIti2VUItgv4cYgMAI40A8JWYjDgHse8HY1oW9HGXGpOa6V9++WUkJiZi06ZNrtdWrFgBo9GIPXsmwic1NTVoaWlBaWkpAKC0tBSnT5/28LLv3r0bFosFRUVFrjHu53COcZ7DZDJhxYoVHmMkScKePXtcY7Tcy4wQegPol18rNpW8eH+4L2mt1FnHl9Plt9hDZ0xG4z6bPw4YlslFNgzdr0uT8IqToeAbTUr9gLEUrv59k8Gu7+u1EPRML0kSXn75ZWzduhUGw8ThUVFR2LZtG7Zv347Y2FhYLBbcc889KC0tdXnLr7zyShQVFeGrX/0qduzYgc7OTnz3u9/FXXfd5Zph77jjDjz//PN48MEHceutt2Lv3r1444038O6777qutX37dmzduhUrV67E6tWr8cwzz2B0dBRf//rXNd/LTMBw8SCGnOA6zATbwEI1dVekCjv8keBKX4XTsjptoJVBAIJ6WJjkGX4/ACOtRhQUfCHGpbJ6bgQtBTZk0i2RWEM7AEudIGIXGC5In8IlRtBG/8EHH6ClpQW33nqrz3u/+MUvwLIsrr/+etjtdpSVleFXv/qV632O4/DOO+/gzjvvRGlpKcLDw7F161b84Ac/cI3JysrCu+++i/vuuw/PPvss5s+fj1//+tcoKytzjbnhhhvQ09ODxx57DJ2dnVi6dCl27drl4dwLdC8zBRP6eQDXgoz8UoNWHiaR++5nf60p480bga42+EOTONYNsUW5AaY3TDTAJro1/ORpqNJQIBu+W0RZGgdMKwC+FjAky2FGI00wElvovdt2AeFbz+/eL3LOK05/sTMVcXoAkBzlwOB2/7np7nBpQRh/KAA/1Xtsqn/d+YmBci4+CwidAPoxqeo5d9RmbNclE0DzFVRWIu65+UwkABPNdfAQ8vA9ho17fbJ3PKeZ9ji9jnYYaZAaPJNAv8jGlTRrTe3Pr7ViDUBAwwyUZstE0Sw3NpE6BIVTALqpE+98EWrkz6l2b1n+tx78cepfMK6i0QMyJDsJ/cAfBxE0JCtdwuhVdtMMERpAhh6Sf+kBePf4fDjdlzLhdOkvttMvtlKLaFUCFOlI7XQGZyLcBCcYakSMmcbIlYQoJiXUoQAbq143r2W74y4zZljkKQWmhu1dIOJObfd3CaIb/XTDOTvZKGX3jfo6+ZgkWQ56tfyCBLqvFeXsO0H+Ly9X5QlyWq/KEl9so5l+QhUNyYnNVHdOqIPfB4ZwJshthgp8BYAouDIInTCh8j0EAaMtEkJsu2jnYB1FdKOfZhjGCGJaCtg07udJF0CStOfCc1kAumj+OptIS27JGI23S90AJPq62ES93Ey09nJaNvn8jR52KqrpLQzK5QSpmc9SDQAtCGdAhBYwQScJXRroe/qZgE0MbnwwZbNsAgBC98bCKTlHv0LeKxtoKMtDKmtQ+7n5quBDiEqIzfAR8XDpAWjMYzAUBnfv9n9qH3uJoRv9NEOEemD8z0EcwQBCAO04d7w14z1wyA0lJytRPSoXvJwHbKLcyXYedSqy86ijUmyifgY21m0r4wcyMqH0owFi+2DSt3yxoy/vpxFCCFXB9Vfh5g2XS3XdNY3NkGfRAGjcCyviUxikBXbiYSFU+E/UMRTRpb+hyP+2w/k5uXxabesvFAgA/AkQsRcMp1RPcGmjz/TTie2vgD3IGYeNDmKsxsyzoFV23BCbab29putE0PAam0z368JpBBTRkHrpGLFDpeDH+35q5FCgcmNLtxMDjo+03fclhm700whxTEKTXQpC0kv0pwnvziQXdFwW1aETqv23r+Ky5Pp6UU771ZIMBPrQcvayJwNUlEPrvYq11PgNxaoaesS2W9u5LjF0o58miNCo3LLJH0yCtuU6QJf2fhtBuBPkP7NhkVxK2yinxzoUOuMYaKIRlyePOwq/mYFKeJ9TrKXnDAbhNPUPGJfS7EMnxlUAGYUknWcNwUWIbvTThdBIZ8hgCCbEpHVpD0CxsYXvxemszmXImXle0lqiLEXNxFDHGxNFM+a0+h+UkBRyF/gjwRs+QMt/pW65wWYK/QyOg2B0pVwfdKOfJpiQfwNCv0i/hMbVVOPNsEie3QI0cNCCq+5d083IQpxKWOQ03Bg6q6utNKQewPxv1IvOHwaIn5ZUWmDj1B8YfLVyz72ACHTFwaW6ugkR23uTv8eLFN17P02QkeeB8b9AVc6KiaRffMbi5l030JCWs520GkzCxF5YC1IPbftsXCHLXmdQpxkRaE2A1n5wUi+tfhPP+e9vpwUuy08vu3FaOsuEBSdC4sJNvMN+AISMg5mKfIOLBN3opwnCO9Vt1QYMA6JX7jmXJu/TTTSpho0GQGjNuNQ7UaVnyPDK4Q+ANEzPw5fTh4zYRH+4zOAch3wlgHBaG6A1D16NQGFMqU2unz8R3HmZaK9tlQ2w/QsIvSa481zE6Mv7aYIJv2USRzlnfAedmfkjdLkqnJYN3izn8kfK1WcrAeMSKljJxPg5r9uqwb0GQKuijgsBMObSpb1QQe9Bi269N2ycNn8HfwIwrQXA0SQfNkXeiqSoJ+oYcuFdeUjs/1Aee4miz/TTBGNaAcJly9pvmg8KMMBOxSLYJN9cdgCAmS7b2Rg5zZUBIAJqjSCkoYmadXYePS7Q7C25Fc7wR+Tcfz5AZqAX3AL/baq9r8fETZT+Olc7Uof8sAvxLFpSEg6x7wchDjAB/76XBvpMP52YlgQ3XnMSjZo33k63B0IFXSHwR0Cba6gYpNRK209xeRPHGZf7XwGIdfQB4fq9kTbyMK6C9jlEqziHGRDqAdKt/DASG6jBc3l0xcPOo/fnDRkFHJ9ovObFj27004kYbAqrVqMJRt/dRGdKNYQGgA2f+J0/RrcPXI76MT7NKWzyrD/fT5TA/ZqNcmPLADn3hkJ67kCItVT+2o/HX8/Fn0A3+mmECf9akEdo/OcIRh+fP0xbRKvioHtndwOUOuSEl9XwqY4DfMN6xpU0m07qpk01jKsA+OmqSwbpOfhj/htjMFryC9zxs4Kw74WuDEfRjX4aYcxXBNkgUuOXXIvijMdpA7W1JgBfDyp24USQHxjZsj6dEyONgztzDbhMGhUgIwBCAXKOzvpsjIaqOAZg1UJpxuBENrg8ebWhIiEu9QL8JPr1XYToRj+NEMfJ4GrAA3Z4lZEGtI3j8uXsOYu8XPdj/MYswKiwpBdraWyeTZLPFSGLVcrxcyZcvm+HZ8KO1EGr5owr6PW9YaLl3Psm5fsxFAQnw83KDyw/8XgSbPHTRYpu9NOJVvXbiQO0DSND0LT/Z0PobO3YRzvJcvOUjzOuprO1s5uuC4aWyBK5pTV/mBbGODEsDqzpz5eDdq5ZCo+VDBur7mAENEQy3HFT1RFOq4uB6n3sAehGP62Q0f8N8giNHWgBDdsGg69yrNhAvfPuGFdMZOSJ9VTairHIJbIp1KiF43JDCbeVAJcZuD+eE9JH/QaGggmlX79ltIbglvaGQrcMQaJeESjUgwhN2s97kaIb/XRiKAg8xh3VbjUKsFH+3zfkufLPPeDL6bIfoDO1d8ab2EkdhUolsqJc0GJcRSMCWisCAVlks5OG1YwlcgGPwqqDiaG5A4pCoip4+yz4U6oPRTL0AIgQRArzRYienDNNEGKj/dkMBfKX0ghXsgyx09ixNCzPlrIqbTCy04FyyZlwlTdE6mcwFMsrAa/VhdQnL72VJLascgktF5w2PxNKK+qIHRCOAggFTKtpDN6nPJhVSTxSw6SgouMADEuVawr4kyC91wBxb4MxaggvXoToRj9djP1ZuzgjE0HTR4ldXn4bQHPuHXTGEzvhU6seSALLr8CGUW4DpRQDd9DsN39NKIxL3NpQacApciGck685DohdNO/fuFruT+/8fAbQBahG/4Zx4UQXHHf8ds6VQEaeARMz823OZgO60U8TJJjZirEof3EnBlBPNxsvrxo4Of9+OdW+Jza6lCcO+v9sDCCpLL2ZBPrA8OdEU9oWuGDV03qdNfl8Bd1+sBb5s5XDx4jFapqOyx+WP1suVfPlEgFhCJqScgDlyABAVzPGlSotsFjAtF7b+S9CdKOfBoh4DnCUaz+Amx9AYorQWdFZEcfOkwUr3bYDxlUTy2IuX9no2TjQ2U+k/6+W/85Gq/sUjUvUH1DGZbIRx9GcfP4Q/Don2Rj6tvOzGZcCQge0G3yMbNQWAAo+ALETdEvlFQo1FCin614i6I68KYaQcZCB/wquVXOwoT0uFb493N1+J2OymKXbvp6JARBGZ3g2QVVXDoD/Cji1HAEue2JWlfoA/hPAdJn/jDu+wnOm5k8AGJYz+jR8NQ051DdiVHGYSm3Ud+ENYwbGfg/i8Le6unjRjX4KIWI3yOC3fKWm/MFlBd9FRumB4q5tL1TSpTLDTDTMNBR61s7zbisR4wo6gzvLcwkPwD0LT8awWD2Zho3CxIwaQmdtx3Fa+6+KnXr1ubwJbzsZoysWxkwluRinLr6RPgzcH1bOkKHYCMV0YTAKxUMhVJkHALFfmmq5+vJ+CiGjLwP2fcEdxMarN3hUgrEohMosymIYZMRzKW5cTZff/BFqQEINXeryh0GXweGyoKQDAPFtraUWlzcs9XzQGYvdPPABcg+ESlCDXuK5/ybjcDn3+MP04cCfBPXMLwbA0fwBgG513NtaA1Tkg4wD/KdyLwF5OW8smAhTBlMOfBGhz/RTBCEOYPyN4A9UmznVMGT6vsYozXIK8CcmymL5I3RidoW1CIAR+rpwyvecxqXKcXnjSkA4qe7801IcZFwUWCFHbIVrCyOc8i1D9lG9FWQtAwlg3TQJ3YtuzqcJyBxGN/qpgj8VfCEMlzOJDjJKX1RW3YvtgcMrqUclAYYJB3j3bD6O5hwowsLHUSbUApDzCPxGAmSkYWivsXce4+UHEevkYhuDnEnotnriTwPsfLlT7oS/gvHbLOPiRTf6KYKMvxX8QWxs8Mco7edJn3/HnDtCFZ2d/cHN98x9Ny6V4/oK8BXy+dxmXsMCuL5aWnLoxXoavtMKl6nsB2EtNPxHhuEZn5doebEh3/P1oIqhLh50o58CCBEmV7YZrANPcT/vRNLWFgqghu/tVedyqNOLsVDxD1dhjRkQ/KXbjtG9uEurngEIC0AWvhQ1ymIF8wBU6wLMnwBgV67O40/6LEiINAqiuoK5eNEdeVMAwxhwtn4ZQs3RCLeMIjJygBaA+JNv5rKCc+ABgCFLdmYpIFRQR500iIBLZTJGZ2NBdv5xaQB46i2XeuGx7Dcup++RVDkJaFi+hteynT9MIwRMON0rG4pkiW2NEtZ8JX3oaCmnVd0yiLL+nkJPAMbs2yBz7DcgUg+Y6J9pu8eLBN3op4h9f8vDn35aBYZh8V5XAlihH4CRhpzYKLqfZEygXnIetvEYmA1tYJx595oIsFTmD9PZ3vmlNyykDxalh49Q4SaKmSCLaHjv8cPoHpgoxeZDADZiwshhlmW5BulMz4RQhxubrFEmYBwwuCUY+UOopPes5A/ha+jf2vvBYMhXTh1WEtK8yNGNfoq49p6r8aefvgUA2PndGFzztTXIzOuiX1DRN5z26I2b0Fy1GJ+7ORqrNxJkFQ4gKqoODPzNdCHy/pnBRKaZSBtDQJBDbSLAsnJRj40avpohCU0ALNQY2DhACoNLHAOQvepqjTBsgGQD4D6rRgGw05AdjPKWIcR/yq87mhtbSHKzDCUn6JDyw0MtikCGQIgEhrl0droM0YXDVLFarYiKisLQ0BAsFv/ecYfNgf+wfA2iQOPSDAO8WWtCZPgJwGs2J4TBtTkrYR/3zKrjOGDNNdFYu8mEopXDSExpAsfISTdMhGwUGgtRuBzqIAMnl7U2qYyTtxnusWwAgAVgpODUawDqJScCLRQyrQYcH4OuUAIYNJdJr8WmyC2uA8Akyko9CnkATJzszJP/vkyMnGOg/LdjIh8FE7418DVnOVq/r/pMP0U0nm5xGTwAxKeaZCFGA7yNvrcrxcfgAUAUgQN/H8SBvztfSUL+igXYcEM4LrsKSIgPovWyK39eBDU6hRx0Np6OY/Kpc8zd6I35QZa4ykhtNDGGyQAch6jnnzh8E33c4fJoCI4Ma3dGkm65PFjhAUH6PIttDDn+P4spgCrvRcals6aZZkaHPGey7/w6ApaoLsW69/aWVJ/X1KgpH8WvHuzG8KDGIhSAXpN3y5ATa30Vc2CSRS2P08w8x2E55z2UzqJqDkMtSL10W2DMl9txdcnefYUkIsMSGpFw5jgI1f7lt7Uinpu4nqSSSThxE+d/vTlE0EZ/7tw5fOUrX0FcXBxCQ0NRXFyMo0cn0idvueUWMAzj8XPVVVd5nKO/vx9btmyBxWJBdHQ0tm3bhpERz2XkqVOnsG7dOoSEhCAtLQ07duzwuZc333wTBQUFCAkJQXFxMd57z7NDKSEEjz32GFJSUhAaGoqNGzeirm56qqsWX1EEo5nGqufnhKBwyUc0JMdwE0o1Mp2t/lpQKTM/M4jwnqEQrpCZE76OeucBuTVUnNeeWFbLMS4CDGnwLejRCkOLYEBoCE1spo0onQ5D4yq4lHSNq2k2n3fNO6sl0QiyQ08lfCd1UOkvNsV/O+2wW8AYVRR0L1KCMvqBgQGsXbsWRqMR77//PqqqqvDzn/8cMTGeX+KrrroKHR0drp8//vGPHu9v2bIFlZWV2L17N9555x3s378ft99+u+t9q9WKK6+8EhkZGSgvL8fTTz+Nxx9/HC+99JJrzCeffIKbbroJ27Ztw/Hjx7F582Zs3rwZFRUTnVB27NiB5557Djt37sShQ4cQHh6OsrIy2GxBzJoaeeGb/wveTpfxbfU2EGfBitQja9Otco3tPhdIktqTtLwQmAxBSFMRpXx3K13mGpfShBqpQznmzx91q9KbBEyiV9MJMhGalDrlYhoOMK1TdxLyp+h5AiLRphmqb/fJ4Ug/t+u3J8DFSVCOvIcffhgHDhzARx+pVyfdcsstGBwcxFtvvaX4/pkzZ1BUVIQjR45g5UqaGbZr1y5cc801aGtrQ2pqKl588UU8+uij6OzshMlkcl37rbfeQnU1TaO84YYbMDo6infeecd17jVr1mDp0qXYuXMnCCFITU3F/fffj29/+9sAgKGhISQlJeGVV17BjTfeGPDzBuPI+9GN/40P3zjo+r2oJAKP/e84YmKcYSID/YKKDfj5g5vxzz9oj9Hf9oMEfPEbWuWbDfDrODMsBYQTgU/DxKiE6vzAZQNg5Iq/AA5AZ/GP6vsaw3dMnOykU8lN4BZRbX2hAYrttY1LwMS+elH0udP6fQ1qpn/77bexcuVKfOlLX0JiYiKWLVuG//mf//EZt2/fPiQmJiI/Px933nkn+vomsrIOHjyI6Ohol8EDwMaNG8GyLA4dOuQas379epfBA0BZWRlqamowMDDgGrNx40aP65aVleHgQWp4jY2N6Ozs9BgTFRWFkpIS1xhv7HY7rFarx49W+to9DaTq0Ai+vgaw2Z06bIIsYgF0NgeXZ77q34JYmRjy4NdTTrxlrv2cx71nXSCMS+k+WmxD4K2BYUKyWg2+nPoB/LSqAiA77RYpv8dlAGKFXFw0Lp8v0+s6Jz3Lki8BgjL6s2fP4sUXX0Rubi7+8Y9/4M4778Q3v/lN/Pa3v3WNueqqq/C73/0Oe/bswU9/+lN8+OGHuPrqqyGKdMnZ2dmJxETPpZvBYEBsbCw6OztdY5KSPJddzt8DjXF/3/04pTHePPXUU4iKinL9pKX5Xxq6k70k0+e18WERQwO++/fO5uDCYKkZQeznVQUxZcRGgAsgamlcTVVvyKD/5TO9oDz+BABRNsAADzVDoWdjDEUk6gcQ22QnoB/1X0nloci6/9uPy+dros5Dd+NXqyu4SAnKbSlJElauXIknn3wSALBs2TJUVFRg586d2LqVxjndl83FxcVYvHgxsrOzsW/fPmzYsGEKb33qeeSRR7B9+3bX71arVbPh3/nMLTi25zRaq8/BYDIgJMyMkHAzmutNaKhYh5FhM8ZGojAylISUrAhERFvQ0diF8WH/s3hmYSiMXBCedDGACo9hoZ8GFcyEwQPUISfaaUotE0bfd9a0AwAiaKmva5nOQ5OoZVAlrZIcifDT8kuslotwmrxeVzFmQf57GhYBjAFk/H3AVApGa4nyHCcoo09JSUFRkWd/ssLCQvzf//2f6jELFixAfHw86uvrsWHDBiQnJ6O727NSTBAE9Pf3IzmZNilITk5GV5dnFpvz90Bj3N93vpaSkuIxZunSpYr3ajabYTZPrsa6tbodg11DYDkWgkPAiEPAyOAovvdl91GeGWQMwyCzOA1hEaHgDByGeq1oq+2AJE4YzbrNgT39BAYA4RCRDEOgvTQRQJf33s4+A/V2i+1edQGCXJm3akJ8gz8qz5S8QgtpDS4i0Z8eoALGwgDCoZDFSJomfjfkK0hje+G8d7EfGN0JRNwV3H3NUYIy+rVr16KmxvMPWVtbi4wM9SVgW1sb+vr6XIZXWlqKwcFBlJeXY8WKFQCAvXv3QpIklJSUuMY8+uij4HkeRiMNg+3evRv5+fmuSEFpaSn27NmDe++913Wt3bt3o7S0FACQlZWF5ORk7Nmzx2XkVqsVhw4dwp133hnMx9bE//v2bzE8ENyyPTE9Hk2nPZfuRrMRC9fm4/R+Gmc/U25BU8M6jI4YYO3nMNjLoK+DoK1BwNnTNvR38hgepMvp7/8+AZdd2Qy/s61YI3vwT7i9GCInsMhORyaOhvfclXKcmvz8Eaok6/gEist4EsD/wKYEEAFVgtDYPRsprzQUPhtfQRWCnfF+H5ksFbhsQGwAGXsdJGwrWFbjcXOYoIz+vvvuw2WXXYYnn3wSX/7yl3H48GG89NJLrlDayMgInnjiCVx//fVITk5GQ0MDHnzwQeTk5KCsrAwAXRlcddVVuO2227Bz507wPI+7774bN954I1JTadLKzTffjCeeeALbtm3DQw89hIqKCjz77LP4xS9+4bqXb33rW7jiiivw85//HJs2bcKf/vQnHD161HUvDMPg3nvvxY9+9CPk5uYiKysL3/ve95CamorNmzdPxd8OACAKIp7++gs4+o/gk1ni58eiq9lz9uftvMdCtqmyG/+5TptDsbikm2bVBfR8u3mqmSga63afsUmfvA0Yw0Q2oR3gCuhSmozDx+CNK6lzjgmhMttERRyEmx+kECgL8GcBWOnihE2izj2+Cp65CDaaoccfAdXB06hTyMbSkKo0CAz/DMTyfTBBt8ieWwTlyFu1ahX++te/4o9//CMWLVqEH/7wh3jmmWewZcsWAADHcTh16hT+4z/+A3l5edi2bRtWrFiBjz76yGPZ/Oqrr6KgoAAbNmzANddcg8svv9wjBh8VFYV//vOfaGxsxIoVK3D//ffjscce84jlX3bZZXjttdfw0ksvYcmSJfjzn/+Mt956C4sWTXhyH3zwQdxzzz24/fbbsWrVKoyMjGDXrl0ICQkuTu6P/330j9jz6uQEFjmD8h7SYRewcG0+8ldmo/ectnrv2GQTIiLkL7pQSQ1PFXlpzybRGVFJDlqolJtOOn+vpjOtcTWtr3efSblMWhpLBmRP+SB9CLAKjSS1KOm4w2XDo/pP6pJj/ZD77blFGMQ2AKysjquleMdAu+wAAGyAcAZk8C4QsRNk5AWQQP6ROYpecOMHLXHPO1c8iPrjQdbFy+Qsy/J77KJ1haj4SNuMdceTCfjCLW6xfLW8dABUiHI5/cL79aJHAYw4UXRjKKIS2FIHXeKLTXSlIDRDWXqLpRV3kp2uEJgoOVsviJBloHg+GHpf4Gh+v3Ep1cETNKy8vJ2azl4AztJcLhtM1BOAccWccPJNS5xex5fYlOhJH9ta0w6DUf3LJAlKmXXKlHzOK+1WOD2hZsMkwGMnx2XTL3vAsNmQbFAyhAekUaoRQEZppEBogqrWHiQ5Di7n0xtXQGNxvdspFAQxPCATkt9cOoBo3269ajBhnr9LAsDMpw9MNoXu9fu/Ergd9xxDN/rzZFXZssCDVLCP2TG/QLn4hmEYnKvTvrxMSlVIdhEaaWadIXNiqW5YAohngyiZZamRA3QbYMigTj8AdL+vQQyUy6OpvY69tGefcTUADQ4zJoreq1bEFno/XBwNx/nF5NnUg0kAwABcDNUhdOvlR4Z/BmLfr/0+Zjm60Z8nq65eqmlcmEW5y2xktPKXn2EAluOwcG1gxdZVn4sCxygJZg5Sx5lQTcNspvXykl9jMY1xFdWNN2S7nXOE5uxLo+rHeZ9DrJ/w2Ivt8p5ckvfkKr3kASoPFixkjO7thQp5xlbJKjQUeqoXM5y8rB+gqxjOraOt41OQ8bepzPlFgG7058m8nBRERAfIggPdv8fPi0VawTwUlU58odTCfJJEMNA1CMEReP975c1+8saF09RQjasBx35oE+EIc/OEAx5LcuKgpariuQDFLKFyaPCI5zWJlRbzODvZSN200QWnVOkW5D6aSfCMzQun5bLeVfBdWbj/XTn54TMKsKn0c0nn5I468oPHvhdk/K8gCipIc41Lq5B4Guht78fIoPKsx3IsEtPiEZVoQVttO/o7BgEA6YXzYA4zw2g2oLnKv5TU2HBgb3fuYk5dMw6QBSX8OcPcYKKpsIaHE9BtbpA66NKZjVff67IpNOtOrYGFR0aeNFG7b8gDECoLbhBACDI91pAF8N5/A0FeWUTLMlrlcojSTSTTuGxCcENspo5HZ+MPgD6UCAdYnwAxvQdE/QwM5y86MrvRjf48+b//fgdhllAkzI9DREwEOAMLh80Ba+8Iult60NnUjc4mz6V3y5lzyF+VDYPRgMpP/GeNWXv975k5DkhOOUDzz42rqBPLu/2UeFZb1Rw7D7QUtsnzdTJCw3tSFwBJ9pKvUm6lZVhE99b+YvFqy2SnA45NlRtianxQOfHXOIQMUiPmMumPw63+gj8qt/220zwDQwjgOApXnT9/knr6IQCOgyC9V4KEbQMTceec8Op7oxv9eRKdaMGYdTzgjO0OwzJoqmxFTFIUCtfkwmgygnfwsPaNoLetD2mF82DtHUZMcjQaTvgPB6YXhE2IafJHAJhpEg1roXtbqZ2GoQyFct93leU9l0eNWKlfnVDl2wDTJxeemUjR9buFYALr/UvtcsONCOp4FKoCOx7ZNG2S4mKTXHRTQHUHnDkKUi/dsrDxsvGHydsQXt7vRwNGE2hTDw7gK0AGt4MwJrDRTwe+7ixCN/rzZN31a/DxXw6h+nB94MEyienx6GrqQWcj/XGHYRi01bTDNmpHd0svkrMS0dmo3vY6o8C7VkCOiYug+3KnA004ox7zNhQDQh389oUXztBlsNQv97V3N+wQuTGkhpmZy9HWG17qpiFFvg9AhJwUVKe+WuFSACmIakTJSv82hkWyNl8cgEi6IiDyw809q9G9HoFLA5hUAN2A2ARp6PsAGw0m/D/BsF5hwFmIbvTnSWp2Mp458CM8sOEJV768GvHzYpGanQxRENHVpLwUJYTANjohHxUVH4nOxm4YjAxyl4Yja5EZyelGzM8WkZDKIznNT8aetz4fXz6RSuvEsETem2tImOGP0xmVTYZrpufy6CwcqAGlE5dgpx/kfPgJRuQHSogcDWj0jd8H24GWm0eNXqgAwNB9PulSyQsIpyXLhnw5jNgGMDa6GuAygXGqDEXAgIm8N7j7uADoRj8FcByHL9xzTUCjt8RH4tT+Kr9jvLnsGuCZt9rBMm6zq7vSq18ELwefSHPimViqImNYLjvNNGbIsXE0lijWA1IMYFpLFW+DybDT0srLmQ/vg02efY30byC20BWBS+47CDwKg4h83lB5RXHaM13YkC13yo0CuFD6sGBjQFt+udUsjP4KJHwrGDZ4DcSZRA/ZTRHz8wMr3Ha3aOzr5sZQrwiW6YHnclrjPxt/jIbXjFcAkOsNpD6ASwYMK+V2z1oNPpGew1mjTgao2KahQNvxAJ3B3ZJeVGFC4L+bD08felI/zfLjNMpmu1BT7hmXVxThclMReU50RhvIEJUaMy6jfw+f7UyI9uq+C4hu9FOEluy5kYERJKYH9wXtalVIW+WrtRmbcQmdnfkPPSWlhCqqG6/FAAEagiNGGrt2h3RDc/MNQFuTSiaONshgwj3j5IoIdFvh0VZbA4Yc+JcU66UPFTZJDtd5+Tr443Kn3SKvA20IOrfgAqAb/RShtegmOVOLyusEXa1K4S2r/x5s3HzqbedPToTO+KP0CwzQmUpsoTFpfxlxgBzG4wGDSmabcFYOZwWC8ewZr4YhAwCZqNgTmxQ0+93gcujKJRiYSG3jpHM0V8AZOXB/AJExuv83XU63BMblAJMAIgS5zbgA6EY/RVzxpVJN47ybYgSi95xKqErqpG2l3WEiZEdXB/W2e8M30hg475bP7i+rjsuie1upVy61VaozsNGwlpvEty8GmnNPAhXPAL4rByMgqkcvwIbJIcXFGvT8nJfwE8/3xrCAGr1hMfVHcAUTAptcJl1pSAO0vp+MzIm4vW70U0TmonS/FXNOGitaNKXtOrH2jkAUVc7rMliOZpvBKDuk1NzjVvlB4fYgUatv5/KocThlo8moHNZTQKyV03KzPWvwARpB4BZ4iVSqEe4rfmFcEqABJlVWgnCKrl6My/1r5rMJ6n391CA86N9UoJEP/ric2ce5RTIMgLEYRBwM7twXAN17P4UoKRMYTAakZichKj4ShAADXUPYcr+IRStaMTIcgRFrOKwD4RgdNmF81IjxUQ62URbjowzGR4GxYYLhoQFExyo0pmCj6NJSbAAYQZtOvVjtGa8X6kGdZm7bCK4AkFp8u8iSEXVVHmc+ALvaLfsvBGAz5LyBRmr8/qrmjErtpM3yg8QrGYiIoKsC9z86kY83y/dZA5+yXy49uJle7KCxe2FI/ru5JR85jgCGFYAgC4nyhwH7XsC8Qvv5LwC60U8RDMOgeH0hxobGEBIeAlEQMdhjRWdjN1rOeDrATn+ahQ3/0QrNvVWMKwBeweiJMGG8klfTRn/wR2ieu1ALKjPlJiZhWCir6ahsQ/gjNKFFqIVitZ7YQmPgwhhdGrty3Hn60GDi1Zf5hPd9TaibSJZRwuDtTANcLbqYCDnf3l3BV7tGAZg4Nz0/Qv/WhiKqk09G6XZCKJcfMMdB24UPaj//BUI3+ikkIjocJ/Z6q8P6cmK/9iYaFBXhCR+v8kkvJVs/55PG4JrhnWIShsWyL4D3/wARKtSz+6ROem4uR0Gkolt2hkUrxNUjPevbAVnoMkBVm7+e9kQummGTAC6Vtsvi/aw0vDGkyRmBbghVNM/BuGTi8/NHJv4eQhDnv0Doe/opJCpeW+PFjrN96GwL0LnFHTWtNsnbwcWDxt01yHhLbRP7b7Fdls+qgksIkz9BM9DUEE7LYhheMJHU0y9W0tnPO+QmNlKDNxTQbEDnst2YN3FtJ2rNKT0+x6CGMV3yvSwI3OjDA6Pyy6Rfjoashmve5I9pjGJceHSjn0IyihSEIFU48WlgcQwXUoeC0GWkcpxdbKX69Vrgj8l73HNyGq17oo4AiL3U8aUEkRNZjBPtycAm0LwAV+RAkD3qCgtKoZrq2HHzaaKQu6CFk0DdemCmIT5O499drKPLckOxNk8/I/sGjKvoQ9GQRwtyYAAtCT5MH2pcOgBBbtBZBSIE0Wz0AqAb/RSSnJWIxVcUwRIXgXl5KUjKUK+5PrY/yMIMw3z/v7vDH1XZ63ojwNU2WinJhvTJ7a1VZjyALpm5bBpJIAzNS+eyqTEYV1JHob84u9hKJ3ulyECgTjhcAt3SBMo1cMLG0RWOcFp+OK6gPgZFzDTFmD8iOy456seQekH/bhYq/MHG0IiIcQ0t4uGSQEZ/re1+LhC60U8hkbEROPVhFax9IzhX2wFR4FH2VeWMsvIPutVDcYo4x8olrH6XvkT2UGsIDYrVAYyyfiKpRxGHXMZrpRl6xmU0miA2yp1wUuBXDNNQLCfLTEKU2ZlkI3YGbEkNwKsZpkTDbk5VIXgl7BgWwHO74e20tNKVA3+Y+jj4Y7QjMBMPjP8Z0sj/AxGCLAKaIXSjn0J2/3afx++ckcXN9xzDf7/DouQqupyMiA7DgsXJSM6KQlNdEHtAqQu0ymyxrC0foL5c6qEhsEAYlij4BrxwikwoHp8vp8HKdfgeAhnO/gKCSiJRCQ3tOT51i/G7PQiJHX5xOgqlNrldVyCUvu42ariM/DB15dt7PQSkQElVDoBxAMJR6q8YewWk93OQhneAaNUTnCF07/0UYhvz/JJ2NfXjry8vR/1JHtt/fhpJ6Utw5ZfrkFskt8o2rvbxXSnD0bg0lz4hLaWlaQR/TKGFlfOUCwDGpE0fHqAZZ1wGTd11Yiii3mrjSlDjOUqTFYxLQDvTnHGTnFpFw47GDDrDMgkAJLpiEdvk1QEAJgkwpMsS3YGiHG4PCKmfOtKIHRNtvZxGztAffzMvscqe/hSa1quUoxAISf73d5cRG/01yNhfQaJ+AMa4BAwXXBr2dKAb/RRy2bWrsfe1jz1ee+tFGj772X2F+P7/HEB0nFtiiBbPM0DDTd4lqVo7vwpnabzZqXHv/FLzJxDcktpGHzzORhDGJfRBAF6OX+fIIhtDymEx/gxgkt/nsulemD8sN8KMhEtKm3QBfBd9TQxUlegec7f716dnUwKH/wDqNJU6favlJC2JT60Tfx93jIuAwbtAjCVA+K1gQj7reotI/WC0FCJNIfryfgpZu3kVvv2//4WtT9zg817Vp+fQ2+3l2BMbNJZiKrTh4o979lhXg1gBLgVDA/F46clr8d4fiydh8DJSG13mG5cBfAUmlikOaqBiu0rGXTiAEVmGyjSxFwYmknl8YBFQUz9Qs0x3FK+hNjZDIZpgAxAoJDvqq+prXAk4PpRXOuUgY69D6rsR0siLkAZuB+m+HMQRpBbgeaLP9FOIwWhA2S2fxe7ff6j4fsr8NrSczUHFkWzUnAhF7bExdDT2IibRiMwiM9LzjUjNYJEwT0RsIg9LzBjCIocRApUSWO99svIgNNdF44FrwzHU2wQAWPO55YiNLQ/y05noVkFNNIN0AWwBgD7qVGNj6D6bMQJ8pRxNMIDWCRTQkB0TIefCV/uej4sFRAW9Pne0rpQAes9ssrZyYu+2167Xo6iH3h+Mm0mxSRNpxfwRWpHnkCXBhSbqAIUAMvYGGJNCzsM0oRv9NHBs9ynEz4tFXEoMzOFmSKKEnrY+3FkWi66mfsDLiMdHRLSfteGTd3zPZTBG4PfloYiNV/iyCafc0ml94YUMvPh4Nt79X89l8oPXGfDSvjiw0CjqYVhEQ1WBNPDEahq64j/1zQp0puOaSgGhhz6wiBVgi5QzCMVO2gDTkCGX5DJw5doTiT5UxCDq6PkT2uP5ahoBbKTbWxF09cDIDzJIoN1y3WoHpC5P6S+hHjAW0nsh4wCbTR96juMg0hAYNkr75zkP9OX9NNBafQ695/pRc7QBpz6sAm/nkZSRIBt8cAi8iF1v+JsFfMN+BBzOnPoMvlQY62PwANBaa8MrTy+Fr6KtJ309yfj0w3+nISmtghv8p54JO+4w4dSZZsimDjwmCuoiHON0OyF2yYq+5yb226QbdGuhcYvCya202MDqRgDUMyCZENCmnLLqkFgz0UdPqKBlvky01zFuUQAujRo7lyY7ceUcANILMvxTkGC2K+eBbvTTwNXf2ODx+5h1HKc+DE4bz533fzeM4SEV3TXhjEfftjFbIX78X5fj3qsGMD6qrmrz+i96UV1xher7u/+6EbdfkYUff70DZ2uCkMQC1FN4DYVUsVZoo/F7Vp7tDYXq5yIq4S7vvH5/OJ2GWtR6mTgVzX4zXZ1w8+UohZdjj1sAOD4C4KA+D+MqGk0QGyaSh4RmuuIRG+DhhCTjwPifQfq+BCI0af9ck0Q3+mlg0+2fwwMv3+X63dqnocmjCjnLsmCJi8Lj2y5TH0TGQBCNf71zFa7LNuGjtwY1nfvBa4dgc3g6nhprC3D/lzbhZ3f1YGRwDA4bjyfvmA/beDAZhGopvAJNxjGk0NmUkZezStV1gKw7r7KHZhScm6qI1Ji1YHBL8uHS5BTc1QAYwLFvQiPQ8yA5T4DIRT7H6QwuVFKHIJcpb2f64DJ2sd33NEINSN91ILZ/BPHZgkc3+mkiNWciNdTaNwKDKTj3SUFJDtIL56P+eCPqjzeh8uA59HSqLE/Fs3juu5/Fc98eBGfQnuVnH5fwg1uTQRCG8bEIvPTktfivDaGoOOAZz26t6cGLP7gyqPunKbwWeAhc8ifkdlflNNNPqKD+CCZa3hJ43bvUoz7TBwNfLstwBYCJlO9lBRXiEFvlPAMBfnsCGApBKxNXQTELUjgj5yXIBs/O99UbdEJGQAbvgWR9EkTtYXie6I68acJ9dieEICEtDh0N/uPEnIFDwZpc9LcPoPqQZ+kpIQT73l2GL22jM0TXuTScOpyP059G4NQBKzrO0lBZ/qps1BxRUnpVpnyvFTt/tBn/eqPZ5d1XYtdvG7HyM1dgXZlyZEIRsUH2WMu5C2oCHIJcwstl03g46aflq/560ysV6PhDUokEsAm0TkAapA8ixz6lQQqvmWjvPCaczvJSh1wYFU7Tcd0TosgQHevcbXGpAdSAQDP6+JNA9LNggtUADIBu9NNEZIxn/N3aOwyDyaDYhTY0MgQ5y7LQcbYbHMei46zyw+G939rQcHoz6k6M4lx9D4jUA8BTBYYzav8nTc1NQUioCW/9yk9Sixu/uNeOvA/SkDQvsHb98FAsQqPyYXAcnNDPUzJ4d8QGOuNLEXSV4K4p7zM2QOqwO0yMumIPl64hKtFJl+hsHF2+E4FmJgo1AIxy5Z0MGaVLd8NiuaeA6yZo2JJIgdOenfDHQXqvBaJ/DsZ8ubZjNKAv76eJ0AjPPefo0BgKSzz3z/HzYlG8nlbDnd5/Br1tfWiubFXV2muv78W/3mhEW203UrOVdXcaTzUjOsF//D4k3IzF64sgiSI6GrW3Xh4dGsdTdxdDFNQfLIQw+OdfN2Lb+iLcdrmEqpPrQfjKwAbvOgFL981CJWgyTKjCoAht0mBO2Hioe/o1mIDURlcWfDk1dLEBrlwF4yJfx5/UDZ/tgDRK1Y3E6uA0+sgAyMCtkHquhjT0HRAShPKPCvpMP030nvMNz53+6AwKS3MxPmxDSEQIag7X+4wb6h1G4ZpcnPnUf7+32KRonKvzDaONj9hgiY9ErCHa1RrbnfxVOehp63V12im6LB9VATrnunPm0Dn89plNuPXbf/N5r6kuH7//xWJ8/BaNnw/1APddbcP8nGLc92wYFi47CiZQlh1rlHP8WzGhb2egjjvGTP+fTZYNiwEYFhOGS+T/N9Kmk0TEH54pxIdv2cAwmWAYBgYjC5ZjwNtFMCyDq75ixrVf0fDBuXS5h58XQgu9D++HiqsAyEAjDWSY3jcZD1hTwFh+CDARII5DNO4v1NMHjbQAUzFP60Y/TdQd811OpixIwpjVhuZK/8tjx3hgB46/CHVXUw8S0+MRmzJh+PHz4xCXGoMar0abVZ/UYOHafFQe0G74bz7TjBXrV2LJaroXt42H4Q/PfQ5/eaEVhWt8v5Rt9Tbcv8mGdZuX4zsvVIJl/OzVxRaF1leCvKyWi17YeP+tsN3oOZeGljPq+REjQ1oVjFgariNOMU6JZtyRcYCEUyFRd5x97sQmauTESnMM2ITAhURMCJjQTWBCNwEAiDQCkCEwwaQS+/8kOlPNYM8QPnzjE5/Xh3qt6G4OrMTacLJJdfnuZKDLf4pqd0svWI5DQno8Fl9RBGuv1cfgXdc70YSUBdqrvySJ4L/vi4RtPBynj6zA3ddchjefbYIoiKg50qC6vaChRD8JQew8bb3uGD+iHl5Y+/0vh0NCNS6X+RZq4FK77LTrkmW365SdhIzJdxnPxNEwoGEh/P0diFeGJcNGTJnBA7rRTwuvfPdPiv3qx6zjyF3ur03TBHHz/FdedTR0IiTcf6VduCUUIWEmnPqwCg6b+urBNmoHy7EwhWg3prFhHj+88xp8+1oBrTUTDzLeziPdj2xY+YHFIEoFREAQRTHaF6jWfv9tt8yqRm+hWgPGVbL4Rg+NLijeTqa2eyS9sl+gkjr6DAUTrbPMZUDIvwMh14CZ5n54+vJ+GvBnsO1nu8GwDIjkP4X07IlmmENNsI8rtbWis21qdjLOnlLWY0tMi0drTTskUVuvuXN1nSheV4DTHykUv3ixcG0+WqvP4eg/lXPf6483IjQyBOPDvrHt7365F+n5S7H92RAULD4Kxr3xhlavdhBz1ciAf3ENUwg1eoctBBXHl+DY/mS0NzJ4bOfbvloDfA1oxaPX5xIb5c5C5zwls/0hn5tJOACGU5dVmw70mX4auHrbvyEkTHkW7m3rQ/7qnIDnGLWOISfAqiAsUsmzTelu7cWiy4NLnz39UTUK1+Sqvp+Wn4rMhWmoPFADa5+6qMSYdRy5yxaovt9SM457rxrAzcsL8cmejRClFOrk0uzV1t4au7/Lv3+k5kQS7v/SJlxXsBCPXG/Dm8824ZO/N2FsVCkDccizEagTpwCH1CFnHC6Sm3EEwJA74wYP6EY/LcTPi8OyjcWq79u9FHbUGLMqq+PkrczGgsUZqDjgf1auPFCDjIXaFXoBoKXqHOK9VioR0eFY8pmFOFfXgaYATkgnrTXtAdt89Xc68MRXe3BtdjIqjin1yVNB0qBiIzMcYKavO2FDxYE28G75E4QQNNaqPPzUinHokVR0E0a5dDaFtgczLJSX88X0dy6dZuiFfknz55hKgjb6c+fO4Stf+Qri4uIQGhqK4uJiHD060RSBEILHHnsMKSkpCA0NxcaNG1FX5xl+6u/vx5YtW2CxWBAdHY1t27ZhZMTzH/LUqVNYt24dQkJCkJaWhh07dvjcy5tvvomCggKEhISguLgY7733nsf7Wu5luhjqVQ9NNVe2Ib0w8P618XQLUhZQh57ByGHhZfmYl5sMzsCpLuvdEQURsSkxPjkD/hi1jiEiJsKVzrtwbT5YAwtRlCAF2JK4M9A1iILV6qsGdxLTE8Hbg/gqKoXOFCAwKrYac0cUlGfks2dUHKnSOT9Kw+xEmzFDCp35xVq3SrzTct+/FoD0gDFeGJ38oIx+YGAAa9euhdFoxPvvv4+qqir8/Oc/R0zMRAXYjh078Nxzz2Hnzp04dOgQwsPDUVZWBpttYh+0ZcsWVFZWYvfu3XjnnXewf/9+3H777a73rVYrrrzySmRkZKC8vBxPP/00Hn/8cbz00kuuMZ988gluuukmbNu2DcePH8fmzZuxefNmVFRUBHUv08W8HN/UydDIECxeX4TY5GifjD01kjITsPiKIkTGRaLykxqcq+tUzOpTY3RgFAsWa+zmKtNU0YJlG4qxYHEGXcr3DqOntRcs678U15v+Tv8JNGkFqSgoyUV7fRf+snMYEx5tPysE40rqBTcUyjNosexwWypXt61AX/9nseftz+Hj3RvUzyOjZvSNVX4elD6ONgN9EBhX0Ow+YpW7/Prp5Gtcrl6CPM0whAR6Fk7w8MMP48CBA/joo48U3yeEIDU1Fffffz++/e1vAwCGhoaQlJSEV155BTfeeCPOnDmDoqIiHDlyBCtX0g+9a9cuXHPNNWhra0NqaipefPFFPProo+js7ITJZHJd+6233kJ1NV3S3nDDDRgdHcU770woT6xZswZLly7Fzp07Nd1LIKxWK6KiojA0NASLRVv3GidVn9biW5c9CgCISYrG/PxU1B8/63JuGUwGhFlCYVVZEcTNi0VKVhLqyht8nHmZC9M0L7MXLM7A2VPNKCrNQ9XBwKITBiPnStgReE+DKCjJRfWh4FZKuSsXoO6oZ85CzjLqq6g/7ime8fuTcUhMbgfYcIAvx/BIMTiDgLAQt062PumtvrS3ZuLrJVGITrRgsNt/TDy9cJ5Pr0EAKFozD7/4y3sKRwACUwKOaQfDxYMKgNbBlUhkKPRsE24ocuvn5wabDDZxv997Cxat39egZvq3334bK1euxJe+9CUkJiZi2bJl+J//+R/X+42Njejs7MTGjRtdr0VFRaGkpAQHD1IF2IMHDyI6Otpl8ACwceNGsCyLQ4cOucasX7/eZfAAUFZWhpqaGgwMDLjGuF/HOcZ5HS334o3dbofVavX4mSzDcsENy7GISrDg9P4qD2+24BCQudBXqz2tIBWrr1mOoR4rKj4+o+i9t41q8wmEW8Jc1X3NVW1ISFMvL42ICUfeymwkZyXi1IdVPgYPALYRDQq8XkgCjR4wDFBYmov0wnly5aCvWs6fX4DcCrocI+Nr8F8bQvHP1718EsIpXx06L+ITaKbiyOAYsorTkbs8CwUluVi4tgDF6wqRlp+K3OVZyF+Voz7TV/SBELryEAUDKsqX47e/+A9sv34T3ng+Gox0jjrr+HJ4dMZlvJyrQpVbFxyPvwyIX//A9BFUyO7s2bN48cUXsX37dnznO9/BkSNH8M1vfhMmkwlbt25FZyf9Yyclee6HkpKSXO91dnYiMdEzEcRgMCA2NtZjTFZWls85nO/FxMSgs7Mz4HUC3Ys3Tz31FJ544gltf4wAHP0nDclIooSwSOWlYltNO1iOBcsyyF+Vg9GhMTRVtqKntQ9GleIcABBUvqjRiVFIyU6EwWDAQNcg2us7YQ6lD87RoTHEpsSAM7AQBc8wXlxqDFJzklHx0Rm/e+CmyjYsWJKBsye1t21qONGElVcuQWtNO84c9L9K+NtLfbjl4QKEhVTj1CcGdLcO45Un+1B2Yx4MBjsYhgfL8AAJ9TtbmUJssMSFw9o3isbTSvXvE6h1IYpJisJrv7oWZ45KqDjQhfERG4BmrC6Lwk13HQPVuasBEAG4hx35ClohSNyyAMVGuaFnglxTLwGhm8EE1Vdv6gjK6CVJwsqVK/Hkk08CAJYtW4aKigrs3LkTW7dunZYbnEkeeeQRbN++3fW71WpFWpqGzikK3Prkzdj/50/R3zGA4X5lb3N/5yBKNi1HzZF6VLrlv9tG7Vh0eQEqPlb2zoeEmRAaEYKUBUkIjw6jGnytfehu6cVgt2d2mHukoLX6HBZfUeRS8QmzhCJnWRaqD9WhrvwsDCYDeLt/f0Egj7wSo0Nj6NKQiQgAb7+chhvvrMbiVRVITF+K7pYBbM4KRc6yItQfb0Tu8iwIvIi7fmxA8UrvXvYThEeFwNoXuBbfNmYHyzJIyU5GbHI0JFFCe0Mnetr68Yef9HjkOXAc8P3fnHXLLbADxmKv7r4OKgfuXbkn1oOq70TIefxBFAxNMUEt71NSUlBU5Om5LCwsREsLfZomJ1PnVVeXZ+VWV1eX673k5GR0d3smYQiCgP7+fo8xSudwv4baGPf3A92LN2azGRaLxeNnsvS09mGohy77WmvaXftYgNbNpyxIQlFpHoZ6hxX3nSMD6l9YS1wkbKN2nD3VjNP7z6DyQA26W5Tz2YcHPB84FR+dQc6yTBSvLwLLsa5sPduoHTFJgYUZa4+exfy84Gao6sP1SM3WVhP+uye7IUopiLAM4ls/E1G4JheFpXloONGInGWZaDjZjMbTLXjrf/0/jFd/TvnfLiEtDoUluVi8vgjF64uQlBGP0IhQxCZH4/RHZ1D5SQ0GuobA23mfzymKwPiYV+KVUh69d3SBTaFefeMSWkzExtKQ3gUiKKNfu3Ytamo8CzNqa2uRkUG9w1lZWUhOTsaePXtc71utVhw6dAilpaUAgNLSUgwODqK8fEKCee/evZAkCSUlJa4x+/fvB89PJFbs3r0b+fn5rkhBaWmpx3WcY5zX0XIv04nTgKISLMhclOYRNstclIaOs12oOlir6hhrqmxVNRTBIUCr/7W/03Pmz12ZDWvfME7vr/J5sCSmaUsUiQpQuuuOOcyM4vWFCFXZ4ngjisDRj+jEsuLygxjpb4XBwGHR5YXoau7FwrVUe++Tv7eg85x6sUxUvPIWKCYxCmcO1eHU/iqc3l+F2qNnMWodQ3uD75YvOsH3Ifjea141CkItrRnw+BBnAdNlcqQhiobunBp9QiWV/9acfTj1BGX09913Hz799FM8+eSTqK+vx2uvvYaXXnoJd91F9eAYhsG9996LH/3oR3j77bdx+vRpfO1rX0Nqaio2b94MgK4MrrrqKtx22204fPgwDhw4gLvvvhs33ngjUlOpHNTNN98Mk8mEbdu2obKyEq+//jqeffZZj6X3t771LezatQs///nPUV1djccffxxHjx7F3XffrfleppPImAjkLM+E0WxEU0UrTn804dENlDPvxDtJxkkw0lv2MTss8ZGYn5dCVXUO16O7RTnO3VbXrin//syndQFrAwwmA4rXF7py/4f7RpCYrtYh1pNffGsIIrMKjCEHd/00FpzRAIZlkJqdBLOc6SiJEv72inpCj9oz8eypZkQohEv72gcwL9dzZreN+oZ2f/vDHgiSl2wZJ//OLaCGblwJgJEFNL2LcWQRDUnbdmc6CMroV61ahb/+9a/44x//iEWLFuGHP/whnnnmGWzZssU15sEHH8Q999yD22+/HatWrcLIyAh27dqFkJCJJ/2rr76KgoICbNiwAddccw0uv/xyjxh8VFQU/vnPf6KxsRErVqzA/fffj8cee8wjln/ZZZe5HjpLlizBn//8Z7z11ltYtGhRUPcynUgiQW9bHxLTPL/sDSeaNBl+U2WLooGzXHA5VXkrstFW2xFQRmuw24rMRYFLTaPiI5GtEvtnWQYLL8tHdGIUTu8/40pS6m7t1Wz0A908mmqiACYSOYvH0d9BnWKjQ2MYs070mNv1+y6MjSrnO5hDlGsOBF5EpkpBUFyqp+JwZ5PvbCyKwKG97ltcltb5c/l0huePyim5I0D4bWAivgmYPwefqjrzxmnTwAtEUHH6S43zidMDwP/94h389bn3kJSRAEmSwLIsbGM2DHQNITE9XlMN+8LL8j2cfAA8nHFaKF5X6LHS8MeqsqU48o8TPq+HWUKxoDgD46M2NJ5qhinUBIPJ4LFFKFidi6E+q6oWoMHIIToxSlFgxJ3krESk58fh4efeRViEHff8+zUIi6S58C3V55CWl4rTH1WBEODOn6Ri89fe9znH335/FX71kHJITE1HMHNRGpoqPPMf0vJT0VrjqVwbk2jEs+9ZERVvQUhIp6/eXeiNYCzfB8NMOD0JXwsy8jxg/wdgXg825n8w1Wj9vupVdtPI9fd9HiWfX4Gv53/T5z2tGXk0VORFkI9pRmMmXebC+TixrxKL1hWAZamTL6s4HWGRoag/3uiR628btaNgURoc4zwYjoHIi6g+7D8kJ/AiUhYkqRp9dIIF8/JSUPVJLTobu7H37bXYsPkE6sobYQoxoqAkF2NDYzi1v8r14Pvb/9jxH19hwbKeM7vRrF5d2HCyWbEKsLmyFVHxkR4p1NGJUS6jD40IQfaSTAiCiNvWt+DL9ybhK3cf9TgHTKVgLI95GDwAMMY8MDHPgQj1wfXVmwb0md4P5zvTO9n9uw+x45bnfV5PykjQFMpKyU7ymD21Ztc58Rf+UxvHMAyWX7kY5f8I3Mpa6/kBwBxqginUiOF+ukLgDCxyVyyAJBI0nGjySJaJTYmCyIuImxeHsyebERUfCd4huAqRClbnoKu5B3nLk9FxtgUxydHoaRnAQC+P8WERnIGF0UwfFvXHmwAALMcAhCCtYD46m7phMHAe/wbOzxKdaEFyZiJsY3bX7F98eQFOu33OpMxYfO/X3cgtkoVFDQvBxP4eDDu99fBqTEtGns7k2PjV9YqltkmZ2rzlCfMmMulylmch1E9JrRKBYu9OBN6z0qy+/CzCLIGv1d85qPle7OMOZBVnIL1wHorXFyI8KhzVh+pRe7TBJzuuv2MIQ70jsPaNoKAkBxkL0zwqDzubupGxMA0VB1rRUjOOkx92wBSeiPFheh5RkGAbtWN0cAwjAyMYGRiBtXcY1r4RiIKAvnP96GruQXSCBfmrsrHo8gKMWccQFR+JwW6rywnrxLvgqKupH798JB0SksFEPgwm7k8XzOCDQV/ezwAH3z4Km0I5bUtVGzgDp5oK6uRcXQcWXV6A/s5B1B9r1PywcKLkhVaCtwswGDlXCm5oZCjS5sUGnMXb6zuRt3IBar1y7NVormrF+LBNMeddid62PowNjyEiijaSSM1OQmxyDKqP1OPEXk+Z7PAo34cUr5DZWHO4ASkLEtFxthuDPVYM9vjG271ToDmvxKTwqDBcdu0VEKN2wBCqLSIzG9Bn+hkgNSdZ0Vs/2GNF9hL1CriMovlYtK4A46M2VHxcjfZ6Gku29g5Tx9qSDDBM4P366NCY6nsMQ7vpLFicAYOB88i5Hx0a9cnwU8M7tdcfQz3DyFupIj2lwtjQOJKzEpG/OgftDV2oOFCtmKZce/SsR0guOjEKkTG+XWcIIYhLVQ87JqTHoe7YWSxaW6BYmpy/Khuv1D6Hm79zHcxzyOAB3ehnhMyFafjyt69VfM9g9lxsGc1GLLq8AGkF89Bc1YaKj6pdS1qGZVC4Jhcp2UkYs47j7MlmJGcFnvWV0oDNYSYUrytEUmYiqg/V46zskfc8bhQjg2OaUm8bTjQhrUC7g6qvvR8Mw4BhGEQoGKUSjaebEegZx9t55CzNhMFkQPaSTISEm3H2tHKtQN2xRkREKffoi0mMBpEIKg5UwxIficXri9DbNpHfcPOj1ysm78wFdKOfIdSUdOrKzyIiOhzz81KweH0RQiPMYDkWrdUTS9+QcDMWX1EES2wEznxa51HwEptCZyslw0xMi0dRaR7yVuUgMpbOfqnZyVi8vgghYWac/ugMOhsnYtHn6js8zsMZOCRnJaLAj4SWOxHR2owXoDLdK8uWIGVBEubnaWshPdw/quka/Z2DSMpIQEdjF0YGRpG9JFNxnH3MjgUq70VEh7nyIbqaenCuvgM9bRNRhxWfW6zpnmcj+p5+hjj0rnJxCG8XsOjyQhzfcxpttTSu3NnYDYZhXKm2qdnJqnF5g5FDRFQYFizNhMPGwxRixNjwODrPdqO7tRfdrTQn3xmGO3uyWTHlFABYlkVRaT5aa84hrWAe2mraUX2oDvNyPdOBFyzOwJh13CN5xWAywGg2IDE9TjXjz0lBSS6sfVYc2XUCANBxthOJ6fGq9QMAEJsSg3k5yeAdAlJzktBer96Zp+VMG1IWJGJ82OaqYixeV4jRoTEfxaHOxm4sXJuPkcExCA4BsSnREAUJrTXtyF2+AGdPNYG3C+hrH0BEdDgKS3IwOjQ+55b07uhGPwOMWsfwf794R/V9b+93d0uvRwnr2VPNLqeTN33tAxgZGguYrOOMu/uj91w/rP0jMJo4j7Hn6jqRvzIbgiBCEiU0V7UiKp6GhMxhJuStzEFrdRtOfViF4nWFikZvMHIoKMlFf8cAWJbxMFpCgMT0BEWjT85MRPy8WJw5VIv+jgEs/ewizM9N9Wv0APWjxCbHgOVYnPgXDamxHIuYpCiPngHdrb2wxEei82wX7OMOnKubSOjpae1D0Zo8VH1Kw6Mjg6M4tf8MYlM8M/fmGvryfgZ4+4V/gLerp1w2V7YiOcuzkCM0wtMLHT9fWQCjvaFTUzGLmsimN45xB5IyfBtfNJxqRsOJJjSeboEoSLD2j2DF55bAHGrC6f1VrkrB2qMNiIydWIIzLINFawsQnRRNnZENXYpbkdqj9bDETTjg5ueloKg0D13N3ag4UO1yFNpGbag4UB1wmd9c2YahXqtHJqIkSkiQ/45GswEZC9OQtzIbY8Pjqn/DlppzPv82RNLutJyN6EY/Axz954mAY7zz0htONHp4/J06+N5IooS0/MB7YrUlvRLh0b7OLW9PueAQcPqjKoRZPMc64/AMy2Dh2nwkpsWj4kC1hxNMFH3zwRw2HpmL0pFVnI68lbRWoOpgrUfhDGfg0N3SizHrOBb4iXoAdNXi7ZgEAEu8hdYGJEShreYcao82oL2+E0az0adGAqAlzt6FT1p7CcxW9OX9DHD/r+/EHcseUE6plWmv7/TYx9tG7R5596PWMSxcW4BKBdnrkHBtM713Zp8qbobGsNTDrvRFz1uRrSjD3Xi6BSnZSaq1BY5xdbkvf0o3mYvS0HCiCQDNXQjUNMS9L4BTz+CoW12Be2ZjT2sfFq0tcPlAAOBzX7sCtz/9VYwN22hyz+AYrL1W1QYkcwXd6GeA1OxkZBWnq6bOGowcQiNCkJAW57GvdU/oYVlGVXZLizquJT4S8amxAY3eFGLE2ZO03XV60Xy0nGnDvNwURQNW+/IP948gc2Ea2hW66gI0913JcReoH0C4W3itr30gYHfftpp2LLq8AANdQ4qafN75C7ZxO6ISLAi3hKKgJBcPvHwXGIaRQ3P+ewvOJXSjnyGyl2ahubIN83KTERIRAsEuYNQ6BmvfCIZ6hsA7BB8jaDjRhMT0eJhCjDAYORzZdcJV9RWdaKFhqbNdHs4nAIiMDUdawTxwHIehXiu6mnpg7R0G0VCpY4mLxLzcZFQdrHM585xOO2/8ZRK6z5jeSKKEpAxfx13t0QYkpMWhp1XZ++89q/vrzxcSbgYY+M0mbK1p91gttMltwJIzE1wGfzGiG/0M8Y2f3AyBF/DPV/YpGktnYzciosIw4jb7hISbkb00E4fePeZaXkfGRiB+fhyGeqxoONnsmuXdi3cyF6YrltIOyCo6EVFhSM1NRkh4CCRBogkocRHIKErDmU9rfargmqvakFaQitbqiRLT2JRov51zu5p6sOjyArAci3N1Hehrp5pwDMNg0eX5ihV5hAApWUmqRt9xtssjTZgmBHnelxPbqB1xqbF+JbBjkqPR5/ZZnSrD/3bTOhiMF69p6I68GSIsMgzbX7oD/37Hlapjkty8xDnLMhEeFYaDbx/16ElXdbAWvW19SM1JhtEtmy8yNgIry5YgzBKqWjvf0dCJ1JxkjFrHUXv0LE59WIX2s11YfEURBF7E6Y/OKEpfA0CUV/ZZclYiBroG/X7mio+rcerDKleXHoCmv/a2D6gWAfkTGOk91++jbmOJU68mi02O9nt/MYnKGXWfvWmt3+PmOrrRzyDv/c8HOPDWYdX3wyxhiJsXi9XXLEP98SbX7HjqQ1o/7k5UXKRHPXj98Uac/uiM3zp9SSJor+/00Nfr7xiAY9wRMKR39mQTzGET3vAzn9ZpVsJxX8ZHRIf7bbw5ah3zEBH1xplZ6KT2aAMWrSvA4vVFyCpOR86yLFc4b2zY/2dSu4++jgunVDsT6EY/g+z948foaVNeuhqMHEwhRowMjODwe8exaG2Bx6ymZPje2MccmirwTCFGZBWnY/EVRShckwuDObAu3ph13GWMIeFm5K/Mdqn9ehMeFYbF64uQvyobDMOgu6UXi68owtINxeAdvIdDTgl/fQDP1XV6SIjxdp4mHu2vQuPpFtQfb0REdDgS0+MDlgV7KwU7+fj/Dvk9bq6jG/0MsmHLOtX3CtfkofyfJ2Efox7xigPVOFfXgbT8VJcUc2fjhOddUkgQYVnGRz+P5Vik5iRj0eUFWHQ5fZAIvIjG0y049WEVznxah4qPtGWZOWRvfUR0OMzhZoRGhHhUoKUsSETxukIIvIBT+6tQc6QB2UtpPL2pohUj/SOwjznQXNmKebkpqg+o1upzqv33BroGUVSa5/FaZ2OXh9Ots6kbY8PjrpWSEuZQE1rOtCm+99FfPlU97mLg4vVWzEKaq5S/ZAAw1Os7a6YVpCI6MQpj1nEUrsn1WIK794ePiAlHVnE6Ohu7cWJvBVZ8bjFVmBkex7naDrTXd7rKctWYn5eC/gDL2oYTTVi2sRin9lW5nH3F6wthH3OAEAl15Y0+qcLmMNq0s6myxRU2G+odxlDvMDgDq6gBCMBvp9268rMwh5lcD0hzWAhiU2M8nHJjQ2N+BUTn5aaodv5tOXOO1h/kX1hZq+lCn+lnAEmS8P0v7MD7v96jOqblzDnkrljg8ZrDxkPkRfS09uLsKRrbdmblRSdSB9bi9YWwjdpxev8Zl9d7bNiGUx9Wof5Yo9+EIHfaatsDquxKEoHIix7Rh+6WXtQebUBduW8cHABqDtehqarV4yHlRBQk1JY3IKs4zWcfX1ve4LMNyCiaj4WX5WN8xIbCNXlYuDYfGUXz0Vp9zsPgnfeavyIbBatzFO8rUDnvvtc/8fv+XEbXyPPDVGjk9Z7rw8+/8SKOatCaU+tGGz8vFrEp0QgJD6Ex975hRMZGIDY5GqyBQ5XXTGkwcgiPCvO7N1Yib8UC1Jb7V7+JjA2HY5z3SMzJWZ6F+mPKRg8Ai9cX4dR+/8U+83JTfPINitcVYqBrEFHxFnBGDg3Hm5C7Igv2MQeMIUZ0Nnajt63PR8bKdd0ritBxtssjeuDEFGKCw+aAfczuUsaNTrBAFERX48v/d+Jnfu95tqGr4c4Cqj6txfb1jwWUw3LSVNmKnGVZPtljvef6Mdg9hPTC+ehq7kHuiiyERoRifNQGiL4JJAIvImNhWlAy2YA2Pf15OSmoPeopH20w+BfZ6GoJLP55rq4DaQXz0Fp9Dpa4CGQuSodt1I7EjAQIdgHNla0YtY6h4kANFixOx3ifDd0tvZifn4q2Gt84PUD7DvS09qnG/Z0s21CMf7t5Hcpu+YwrFVpJ3uxiQTf6aeSXd/1as8EbTAYkZyYgRiW2LPAiohIsKF5XCN7OY7h/BJb4SAyoiFK21QbOTXdiNBsxPz8FplAT4ufFqkpUh0eFISImwmdmrS0/i7jUGFXHWVdTj+LDzOf80WEoXldI2059WAWj2QBzqBkjg6PgDBwKS3JhDDGit60PDhtPHYohJkREh2Nk0Lf3n+DQ1kxiVdlSXPX1z7p+ZxgGoRrqGeYqutFPI1fd+m94/p7feLzGsgzi58chNjkaplATeAePwa4hdDX3oq22A221HchanI7GUy1UBz43BWAYdJztwvE9p32u4S1w4aS/YwDJmYnoPdfnk3ATER2O1OwkhESEYGRwFC1nzqHxFC10WXxFkY/Rh0aEIGf5AjRVNKPqoK/TTRIlzMtN8estN6q04nKKXAi8iI76TnAGzpVlyNsFMAyDvBULEBIRgoqPqyGJEvJWZkPgBUQlWGi7q9xknHXLTnSixZ8RFhmKz9+pnjB1MaLv6f0wFXv6h8p+CN7Og0gE1r5hdDX3uLzOaiTJM361n2ISJ6ERIZifn6LqSFvy2YVwyHtgh82BvvYBv8vdqPhIjFrHITgEVyvrsyebPWbS/JXZqPFa4kfGRsA2alPNtPP2MySmxyM5KxGt1eeQmpOMygM1yF2ehTo/vgFnZ5rE9Hj0tQ8gMjYcCfPjUHesEYvWFvhU/PnL43ey4Svr8PDvfJuRzEX0Pf0sIX5eLP75yr6gjulqolrsWhgfsfkIbjhhGAadZ7s194YHaDht8fpCgGFQf7xRs19guH8EC9fmK1bjpeYkIyYpCgaDASc/rETRZXk4c7DOlakXmxyDRZcX+DgkvXE2r0zOSkRMcjQYhsGYdQyhkSGKpURjQ4GFQ/JWBKfKezGgh+ymmf/4r6smdVwwGmzD/SNYvL4IlvhIj9cJIUjKCE4jH4BLfkstNVdSWRw6u9YA9IGTvzIb2Usz0V7ficoDNS5NvXO1HR5Ow4aTTejx44V34tS8G+odRmiEGdWH6tBy5hwWLM5E7dEGn9beo9bASr5LP7vI7/sXI7rRTzNWhaQbLYhBqLM0V7Wh/gSVc3auEFIWJGLh2nywHONRmKOFsyebkZqtXD+emB6PZoWwIkAFKbOXZqJ4XSEyF81HzdEGl+gFAHQ19yB7aSZYjkVCmqf8l5aHU2dTFxatLUBzZSvs4zzSC+dh6b8tgsPmQFSCRbG1d7gfWa2cZVmqmX8XM/ryfpoxK7Sz0sLokK83Wg1JlJBROB9nDtUhzBKK/NU5qDva4MqOK7osP+DS2Zu41Fi0KwhuOMYdPo0tYpKikZyViPCoMFR/WouRoTEsXq9cJxASbqZtpbzKclvOtIHlWL9SVNGJUa54/+jgGNUVSIhCXflZpBfOg5J7KjwqzKNGgGUZLFiSiWX/tgifvely9T/ARYxu9NNMoPJONQI5oLwxykUzY9Zx1Byu93jPpjErz52Gk7TE1Vlj7sTaN4xlGxbDYXdAEiR0Nfeg91w/BroGwXIsouIjgSHl3u4Arc6LTrD4KAAPdluRvzrH597dcV/+97T2ItwShrY6mknYcuYcTCHdPvecMD8O6QXzULyuEAUluchdnqXqA7lU0I1+GpEkCa/++P+CPi46waLYW80f4yPqTit/EtpqjFnHsejyAtQda8S83GRExETAPmrDudoOlO9Wzi6URAlpBfMw0DWE7pZexUw9SZQwPz9VsellIJ0a95Bc9rIsVH58BiODQGFpHs4crIXAi7jiCyWIiApDYWkeCtfkYb5X/b2ObvTTyun9Z/DB7/cHfVxCepzL6JX6qCvhLybNsgxSs5MDGn1McjQS5sfBHGYCbxcw0DUI3s7j7MlmFJXmudJV/dFa0+5qyqmmPtN4uhkGk8Enrl5/vBGWuAjFPH0A6DvXD5ZlIEkEFR+dweqrl4EAWFiaj699/8soKs3zW6uvQ9GNfhrRokfvJCYpCrGpsQiLCAHLMViwOAPdrX3gOAaOcT5gZl9kbCQAz9x1g8mAnGVZaKtpR92xs+AMrGs/zhlYzM9LhSU+EoJDQFdTD/o7B1Uz/LQW7gx0DroEK2uPNiA2JRr9HZ7nHO4f9Wgi4UTgRWQuTFfM02cYBqk5tCVX/upcLFqbr7kdlo4nutFPI1nF6YiMCcfwwIRTLi41BnEpMQiPDgdv5zEyMIru1l4MdA2pas4pJcMAtPglOTMJoZEhPnnzBiOH+bkp6GrqdiXWLL6iCCBUPCIyNiKo3PzG0y0uUc5AOJOPJFHC/LxUH6MHALtK84/+TprVFxJmRu6KBVh4WT4Wri1AUWkeLHGRisfoBIdu9NOI0WTEui+VouF4E4b7h9F7rh997QPoax9AdIIFwwOjmnLzGYNvZHV+Xgraajsw3E+r4pKzErFwbQHaatsx1GPFgiWZPoUx7kYeGhGi6KjzR1SCRZPRnz3VjNScZL9trhuON3rIYCdnJWLhZfkoXJOHgpIcZC/JvKjFKS8k+l91mslamI73XvrA5/XBHisKSnJRfShwqm3d0QZkL82Ew+YACBBmCfXZX3c2dqOzsRvLNhTDPmZHdwDv//iITTF11R/OllXuSThqJKbFY6jbqriaMIeakLcyG6uvWe5qXxWbPLf7w80ldKOfZjbfczWqj9Rhzx8+8nlPFAI3qaDjJDScaELhmjyc+VS5YYaT9oZOdDVpS7sdtY4FHuSGw8Yjf3UuTgeojQeAmqP1rgq/xPR4FK7JRcHqXBRdlo+8FQv0WfwCov/lZ4Crb92gaPT1xxoRPz/Oo8+bP4Z61JfLTrqaejSJYQB0n56clejRoz4Qfef832uYJRT5q3KQvyoHi9YVIG9FtqrUtM6FQTf6GWDJZxYiNTvJJ8ONECA1O0mz0UfGRgJaetEF0ZklMT0+KKNvb+hC5qJ0NFW0gGUZpBXMQ2EJncELSnKRUTQfLKtnd89mdKOfIbZ894t4+usv+LzecuacK64dGG1V0HXlZz063vijrbYDDANoKbCOjI3AwsvyseLKJUgvnI/8VdkIt/iXs9aZfQT1SH788cfBMIzHT0HBRPeVz3zmMz7v33HHHR7naGlpwaZNmxAWFobExEQ88MADELz2tvv27cPy5cthNpuRk5ODV155xedeXnjhBWRmZiIkJAQlJSU4fNiziYTNZsNdd92FuLg4RERE4Prrr0dXl4ZZcppY98U1KFyT6/P6YPcQcperN3dwp7b8rKa03mCq6/o7BpC91Pf6nIFD7vIsfP4/P4f7f30nfl35C/xfz//ih28/jM13X43lG4p1g5+jBD3TL1y4EB98MOGNNhg8T3HbbbfhBz/4gev3sLCJL4Yoiti0aROSk5PxySefoKOjA1/72tdgNBrx5JNPAgAaGxuxadMm3HHHHXj11VexZ88efOMb30BKSgrKysoAAK+//jq2b9+OnTt3oqSkBM888wzKyspQU1ODxETaGuq+++7Du+++izfffBNRUVG4++67cd111+HAgQPBfuQpITQ8BM98/CP8/cV/4uXv/tGjY6qooGGvBJGIagqrN40VLYpZb0qYw0yYn5fi2ovnr8pGzrIsmEJ8+7vrzH2CUs55/PHH8dZbb+HEiROK73/mM5/B0qVL8cwzzyi+//777+Pzn/882tvbkZRESzd37tyJhx56CD09PTCZTHjooYfw7rvvoqKiwnXcjTfeiMHBQezatQsAUFJSglWrVuH5558HQHPc09LScM899+Dhhx/G0NAQEhIS8Nprr+GLX/wiAKC6uhqFhYU4ePAg1qxZo+nzToVyjhLNVa24b933XEk7DAPEpsb6yDgrEZcag4HOwYC154Bn/3V3YpOjUVCSi7wV2chfTY3cXzssnbmB1u9r0B6Xuro6pKamYsGCBdiyZQtaWlo83n/11VcRHx+PRYsW4ZFHHsHY2MSMdvDgQRQXF7sMHgDKyspgtVpRWVnpGrNx40aPc5aVleHgwYMAAIfDgfLyco8xLMti48aNrjHl5eXged5jTEFBAdLT011jlLDb7bBarR4/00FGURpKr12Fz9xwGWKSokAIMC9HWevOm772ASxYmqlprOAQYDQbUVSah+u+tQmP/vFe/KHxV3i9/X/wxF8fxJbvXo+VVy7RDf4SI6jlfUlJCV555RXk5+ejo6MDTzzxBNatW4eKigpERkbi5ptvRkZGBlJTU3Hq1Ck89NBDqKmpwV/+8hcAQGdnp4fBA3D93tnZ6XeM1WrF+Pg4BgYGIIqi4pjq6mrXOUwmE6Kjo33GOK+jxFNPPYUnnngimD/JpGmpasN3/ngvJFHCT776S7TWtAesJ3diClHvPTcvlya75K3MRmFJLrKX6pltOp4E9W24+uqrXf+/ePFilJSUICMjA2+88Qa2bduG22+/3fV+cXExUlJSsGHDBjQ0NCA7e/ZrkT3yyCPYvn2763er1Yq0tLRpudYvP33K9f9bn7gBj1z1o4D15E7qjzUizBIKc6gJ+atyXMv0gpIcWGL1/HQd/5zXFBAdHY28vDzU1yt/UUtKSgAA9fX1yM7ORnJyso+X3elRT05Odv3X28ve1dUFi8WC0NBQcBwHjuMUx7ifw+FwYHBw0GO2dx+jhNlshtk8OaWb82HF5xbjszet9enw4k54VBhyVyxAntyqKX9VNhLTg9e/09E5L6MfGRlBQ0MDvvrVryq+73T4paRQIYPS0lL8+Mc/Rnd3t8vLvnv3blgsFhQVFbnGvPfeex7n2b17N0pLSwEAJpMJK1aswJ49e7B582YA1JG3Z88e3H333QCAFStWwGg0Ys+ePbj++usBADU1NWhpaXGdZzbBMAwe/v038eSWZ9HfMYjRIdqjPW/FAuStzEbeqhzMy0n26MyqozNpSBDcf//9ZN++faSxsZEcOHCAbNy4kcTHx5Pu7m5SX19PfvCDH5CjR4+SxsZG8re//Y0sWLCArF+/3nW8IAhk0aJF5MorryQnTpwgu3btIgkJCeSRRx5xjTl79iwJCwsjDzzwADlz5gx54YUXCMdxZNeuXa4xf/rTn4jZbCavvPIKqaqqIrfffjuJjo4mnZ2drjF33HEHSU9PJ3v37iVHjx4lpaWlpLS0NJiPS4aGhggAMjQ0FNRxk2VkcIQ0VbUSURRn5Ho6Fxdav69BGf0NN9xAUlJSiMlkIvPmzSM33HADqa+vJ4QQ0tLSQtavX09iY2OJ2WwmOTk55IEHHvC5gaamJnL11VeT0NBQEh8fT+6//37C87zHmH/9619k6dKlxGQykQULFpCXX37Z515++ctfkvT0dGIymcjq1avJp59+6vH++Pg4+a//+i8SExNDwsLCyBe+8AXS0dERzMedcaPX0TkftH5f9Q43fpiuOL2OznQwbXF6HR2duY1u9Do6lxi60evoXGLoRq+jc4mhG72OziWGbvQ6OpcYutHr6Fxi6Eavo3OJoRu9js4lhm70OjqXGLrR6+hcYuhGr6NziaEbvY7OJYYunuYHZwHidAlk6uhMJc7vaaDCWd3o/TA8PAwA06aTp6MzHQwPDyMqSr1/oF5P7wdJktDe3o7IyMgLKlXlFOhsbW29aOr6L7bPNBs+DyEEw8PDSE1N9dtPUJ/p/cCyLObPn3+hb8OFxWK5KAzEnYvtM13oz+NvhneiO/J0dC4xdKPX0bnE0I1+DmA2m/H973//gmjyTxcX22eaS59Hd+Tp6Fxi6DO9js4lhm70OjqXGLrR6+hcYuhGr6NziaEb/QVg//79+Pd//3ekpqaCYRi89dZbAY959dVXsWTJEoSFhSElJQW33nor+vr6PMa8+eabKCgoQEhICIqLi30agU4X0/F5XnnlFTAM4/ETEhIyjZ9igsl8nhdeeAGFhYUIDQ1Ffn4+fve73/mMuVD/Pt7oRn8BGB0dxZIlS/DCCy9oGn/gwAF87Wtfw7Zt21BZWYk333wThw8fxm233eYa88knn+Cmm27Ctm3bcPz4cWzevBmbN29GRUXFdH0MF9PxeQCa3dbR0eH6aW5uno7b9yHYz/Piiy/ikUceweOPP47Kyko88cQTuOuuu/D3v//dNeZC/vv4ML0t9XQCAYD89a9/9Tvm6aefJgsWLPB47bnnniPz5s1z/f7lL3+ZbNq0yWNMSUkJ+c///M8pu1ctTNXnefnll0lUVNQ03GFwaPk8paWl5Nvf/rbHa9u3bydr1651/T5b/n0IIUSf6ecApaWlaG1txXvvvQdCCLq6uvDnP/8Z11xzjWvMwYMHsXHjRo/jysrKcPDgwZm+3YBo+TwAMDIygoyMDKSlpeHaa69FZWXlBbpj/9jtdp+tR2hoKA4fPgye5wHMrn8f3ejnAGvXrsWrr76KG264Af+/nfsLaaqN4wD+he2sVktBke1kmCVMLBRHKS4b0lX/BEmKgrKFeCEUGFTQTUpQdDMluggMc2FdSEIQFZQi3iQSdBo4XM2aZii4/niTSFbb7714aS9Tqa0Xtsn5fsAb95zn/B4evuec/Xkek8kEm82G7OzsuMfP2dlZWK3WuOOsVitmZ2dTXe4fJTKe4uJidHd34+HDh7h37x6i0Sh27dqF6enpNFa+sr1796KrqwuapkFE8PLlS3R1deHHjx/4/PkzgMyaH4Z+FQgEAmhpaUFrays0TcPTp0/x/v17NDc3p7u0v5LIeJxOJ06ePIny8nLU1NTgwYMHyMvLQ2dnZxorX9mlS5ewf/9+VFVVQVEU1NXVwe12A8Bvl7imC5fWrgLXrl1DdXU1Lly4AAAoKyvD+vXr4XK5cOXKFaiqCpvNhnA4HHdcOByGzWZLR8m/lch4llIUBQ6HA+/evUt1uX9kNpvR3d2Nzs5OhMNhqKqKW7duYcOGDcjLywOAjJqfzLsM0TILCwvL7hgGgwHAf1sjOZ1ODA4OxrUZGBiA0+lMTZFJSGQ8S0UiEfj9/hUvCJlCURRs2rQJBoMBvb29qK2tjY0zo+Yn5R8dknz9+lV8Pp/4fD4BIB0dHeLz+WRqakpERC5evCgNDQ2x9l6vV4xGo9y8eVNCoZA8f/5cdu7cKZWVlbE2w8PDYjQaxePxyOvXr6WtrU0URRG/378qx3P58mV59uyZhEIh0TRNjh07JmvXrpWxsbGMG08wGJS7d+/K+Pi4vHjxQo4ePSo5OTkyOTkZa5PO+VmKoU+DoaEhAbDsz+12i4iI2+2WmpqauGNu3Lgh27ZtE7PZLKqqyvHjx2V6ejquzf3798Vut4vJZJLt27fLkydPVu14zp49KwUFBWIymcRqtcqBAwfk1atXGTmeQCAg5eXlYjabJSsrS+rq6uTNmzfL+k3X/CzFpbVEOsP39EQ6w9AT6QxDT6QzDD2RzjD0RDrD0BPpDENPpDMMPVGK/M2OPEuJCDweD+x2O9asWYP8/HxcvXo1qT644IYoRX7tyNPY2Ij6+vq/6qOlpQX9/f3weDwoLS3F3Nwc5ubmkuskLb8DJNI5rLAjz7dv3+TcuXOyceNGWbdunVRWVsrQ0FDs9UAgIEajccWf+CaDj/dEGeLMmTMYGRlBb28vRkdHceTIEezbtw9v374FADx69Ahbt27F48ePsWXLFhQWFqKpqSnpOz1DT5QBPnz4AK/Xi76+PrhcLhQVFeH8+fPYvXs3vF4vAGBiYgJTU1Po6+tDT08P7ty5A03TcPjw4aTOxff0RBnA7/cjEonAbrfH/X9xcRG5ubkAgGg0isXFRfT09MTa3b59Gzt27EAwGERxcXFC52LoiTLA/Pw8DAYDNE2LbSjyi8ViAQCoqgqj0Rh3YSgpKQHw75MCQ0+0ijgcDkQiEXz8+BEul2vFNtXV1fj58ydCoRCKiooAAOPj4wCAzZs3J3wurqcnSpH5+fnYHn8OhwMdHR3Ys2cPcnJyUFBQgBMnTmB4eBjt7e1wOBz49OkTBgcHUVZWhoMHDyIajaKiogIWiwXXr19HNBrF6dOnkZWVhf7+/sQL+V+f/RNRwv60I8/379+ltbVVCgsLRVEUUVVVDh06JKOjo7E+ZmZmpL6+XiwWi1itVjl16pR8+fIlqTp4pyfSGX5lR6QzDD2RzjD0RDrD0BPpDENPpDMMPZHOMPREOsPQE+kMQ0+kMww9kc4w9EQ6w9AT6cw/2KIWNjIKGxUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# INITIAL STATE\n", + "\n", + "state_population = 0\n", + "for node in graph:\n", + " state_population += graph.nodes[node][\"total_pop\"]\n", + "ideal_population = state_population / NUM_DISTRICTS\n", + "\n", + "initial_assignment = recursive_tree_part(\n", + " graph, range(NUM_DISTRICTS),\n", + " pop_target=ideal_population,\n", + " pop_col=\"total_pop\",\n", + " epsilon=constants.POP_EQUALITY_THRESHOLD)\n", + "\n", + "initial_partition = Partition(graph, initial_assignment, rba_updaters)\n", + "\n", + "visualize_partition_geopandas(initial_partition)" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "metadata": {}, + "outputs": [], + "source": [ + "# PROPOSAL METHOD (weighted to prefer county lines)\n", + "\n", + "weighted_recom_proposal = partial(\n", + " recom,\n", + " pop_col=\"total_pop\",\n", + " pop_target=ideal_population,\n", + " epsilon=constants.POP_EQUALITY_THRESHOLD,\n", + " node_repeats=2,\n", + " method=partial(\n", + " bipartition_tree,\n", + " spanning_tree_fn=get_county_weighted_random_spanning_tree)\n", + ")\n", + "\n", + "# CONSTRAINTS\n", + "\n", + "# NOTE: we said we wouldn't have a compactness constraint but GerryChain uses one in their example\n", + "# showing that maybe it's necessary even for ReCom. This keeps the proposals within 2x the number of\n", + "# cut edges in the starting one.\n", + "compactness_bound = constraints.UpperBound(\n", + " lambda p: len(p[\"cut_edges\"]),\n", + " 2 * len(initial_partition[\"cut_edges\"])\n", + ")\n", + "\n", + "pop_constraint = constraints.within_percent_of_ideal_population(initial_partition,\n", + " constants.POP_EQUALITY_THRESHOLD)\n", + "\n", + "vra_constraints = [\n", + " constraints.LowerBound(\n", + " lambda p: p[f\"num_{minority}_vra_districts\"],\n", + " num_districts\n", + " )\n", + " for minority, num_districts in vra_config.items()]" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [], + "source": [ + "chain = MarkovChain(\n", + " proposal=weighted_recom_proposal,\n", + " constraints=[\n", + " pop_constraint,\n", + " compactness_bound\n", + " ] + vra_constraints,\n", + " accept=accept.always_accept,\n", + " initial_state=initial_partition,\n", + " total_steps=1000\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d0cbfdcdcdca47ecaebd5f836f1668b2", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/1000 [00:00 2\u001b[0m [\u001b[39msorted\u001b[39;49m(partition[\u001b[39m\"\u001b[39;49m\u001b[39mdistrict_gerry_scores\u001b[39;49m\u001b[39m\"\u001b[39;49m]) \u001b[39m+\u001b[39;49m [partition[\u001b[39m\"\u001b[39;49m\u001b[39mgerry_score\u001b[39;49m\u001b[39m\"\u001b[39;49m]]\n\u001b[1;32m 3\u001b[0m \u001b[39mfor\u001b[39;49;00m partition \u001b[39min\u001b[39;49;00m chain\u001b[39m.\u001b[39;49mwith_progress_bar()],\n\u001b[1;32m 4\u001b[0m columns\u001b[39m=\u001b[39m[\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mdistrict\u001b[39m\u001b[39m{\u001b[39;00mi\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(\u001b[39m1\u001b[39m, NUM_DISTRICTS \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m)] \u001b[39m+\u001b[39m [\u001b[39m\"\u001b[39m\u001b[39mstate_gerry_score\u001b[39m\u001b[39m\"\u001b[39m]\n\u001b[1;32m 5\u001b[0m )\n", + "Cell \u001b[0;32mIn[44], line 2\u001b[0m, in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 1\u001b[0m data \u001b[39m=\u001b[39m pd\u001b[39m.\u001b[39mDataFrame(\n\u001b[0;32m----> 2\u001b[0m [\u001b[39msorted\u001b[39m(partition[\u001b[39m\"\u001b[39m\u001b[39mdistrict_gerry_scores\u001b[39m\u001b[39m\"\u001b[39m]) \u001b[39m+\u001b[39m [partition[\u001b[39m\"\u001b[39m\u001b[39mgerry_score\u001b[39m\u001b[39m\"\u001b[39m]]\n\u001b[1;32m 3\u001b[0m \u001b[39mfor\u001b[39;00m partition \u001b[39min\u001b[39;00m chain\u001b[39m.\u001b[39mwith_progress_bar()],\n\u001b[1;32m 4\u001b[0m columns\u001b[39m=\u001b[39m[\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mdistrict\u001b[39m\u001b[39m{\u001b[39;00mi\u001b[39m}\u001b[39;00m\u001b[39m\"\u001b[39m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(\u001b[39m1\u001b[39m, NUM_DISTRICTS \u001b[39m+\u001b[39m \u001b[39m1\u001b[39m)] \u001b[39m+\u001b[39m [\u001b[39m\"\u001b[39m\u001b[39mstate_gerry_score\u001b[39m\u001b[39m\"\u001b[39m]\n\u001b[1;32m 5\u001b[0m )\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/tqdm/notebook.py:259\u001b[0m, in \u001b[0;36mtqdm_notebook.__iter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 257\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m 258\u001b[0m it \u001b[39m=\u001b[39m \u001b[39msuper\u001b[39m(tqdm_notebook, \u001b[39mself\u001b[39m)\u001b[39m.\u001b[39m\u001b[39m__iter__\u001b[39m()\n\u001b[0;32m--> 259\u001b[0m \u001b[39mfor\u001b[39;00m obj \u001b[39min\u001b[39;00m it:\n\u001b[1;32m 260\u001b[0m \u001b[39m# return super(tqdm...) will not catch exception\u001b[39;00m\n\u001b[1;32m 261\u001b[0m \u001b[39myield\u001b[39;00m obj\n\u001b[1;32m 262\u001b[0m \u001b[39m# NB: except ... [ as ...] breaks IPython async KeyboardInterrupt\u001b[39;00m\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/tqdm/std.py:1195\u001b[0m, in \u001b[0;36mtqdm.__iter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 1192\u001b[0m time \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_time\n\u001b[1;32m 1194\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m-> 1195\u001b[0m \u001b[39mfor\u001b[39;00m obj \u001b[39min\u001b[39;00m iterable:\n\u001b[1;32m 1196\u001b[0m \u001b[39myield\u001b[39;00m obj\n\u001b[1;32m 1197\u001b[0m \u001b[39m# Update and possibly print the progressbar.\u001b[39;00m\n\u001b[1;32m 1198\u001b[0m \u001b[39m# Note: does not call self.update(1) for speed optimisation.\u001b[39;00m\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/gerrychain/chain.py:67\u001b[0m, in \u001b[0;36mMarkovChain.__next__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mstate\n\u001b[1;32m 66\u001b[0m \u001b[39mwhile\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcounter \u001b[39m<\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtotal_steps:\n\u001b[0;32m---> 67\u001b[0m proposed_next_state \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mproposal(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mstate)\n\u001b[1;32m 68\u001b[0m \u001b[39m# Erase the parent of the parent, to avoid memory leak\u001b[39;00m\n\u001b[1;32m 69\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mstate\u001b[39m.\u001b[39mparent \u001b[39m=\u001b[39m \u001b[39mNone\u001b[39;00m\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/gerrychain/proposals/tree_proposals.py:48\u001b[0m, in \u001b[0;36mrecom\u001b[0;34m(partition, pop_col, pop_target, epsilon, node_repeats, method)\u001b[0m\n\u001b[1;32m 42\u001b[0m parts_to_merge \u001b[39m=\u001b[39m (partition\u001b[39m.\u001b[39massignment\u001b[39m.\u001b[39mmapping[edge[\u001b[39m0\u001b[39m]], partition\u001b[39m.\u001b[39massignment\u001b[39m.\u001b[39mmapping[edge[\u001b[39m1\u001b[39m]])\n\u001b[1;32m 44\u001b[0m subgraph \u001b[39m=\u001b[39m partition\u001b[39m.\u001b[39mgraph\u001b[39m.\u001b[39msubgraph(\n\u001b[1;32m 45\u001b[0m partition\u001b[39m.\u001b[39mparts[parts_to_merge[\u001b[39m0\u001b[39m]] \u001b[39m|\u001b[39m partition\u001b[39m.\u001b[39mparts[parts_to_merge[\u001b[39m1\u001b[39m]]\n\u001b[1;32m 46\u001b[0m )\n\u001b[0;32m---> 48\u001b[0m flips \u001b[39m=\u001b[39m recursive_tree_part(\n\u001b[1;32m 49\u001b[0m subgraph\u001b[39m.\u001b[39;49mgraph,\n\u001b[1;32m 50\u001b[0m parts_to_merge,\n\u001b[1;32m 51\u001b[0m pop_col\u001b[39m=\u001b[39;49mpop_col,\n\u001b[1;32m 52\u001b[0m pop_target\u001b[39m=\u001b[39;49mpop_target,\n\u001b[1;32m 53\u001b[0m epsilon\u001b[39m=\u001b[39;49mepsilon,\n\u001b[1;32m 54\u001b[0m node_repeats\u001b[39m=\u001b[39;49mnode_repeats,\n\u001b[1;32m 55\u001b[0m method\u001b[39m=\u001b[39;49mmethod,\n\u001b[1;32m 56\u001b[0m )\n\u001b[1;32m 58\u001b[0m \u001b[39mreturn\u001b[39;00m partition\u001b[39m.\u001b[39mflip(flips)\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/gerrychain/tree.py:348\u001b[0m, in \u001b[0;36mrecursive_tree_part\u001b[0;34m(graph, parts, pop_target, pop_col, epsilon, node_repeats, method)\u001b[0m\n\u001b[1;32m 346\u001b[0m min_pop \u001b[39m=\u001b[39m \u001b[39mmax\u001b[39m(pop_target \u001b[39m*\u001b[39m (\u001b[39m1\u001b[39m \u001b[39m-\u001b[39m epsilon), pop_target \u001b[39m*\u001b[39m (\u001b[39m1\u001b[39m \u001b[39m-\u001b[39m epsilon) \u001b[39m-\u001b[39m debt)\n\u001b[1;32m 347\u001b[0m max_pop \u001b[39m=\u001b[39m \u001b[39mmin\u001b[39m(pop_target \u001b[39m*\u001b[39m (\u001b[39m1\u001b[39m \u001b[39m+\u001b[39m epsilon), pop_target \u001b[39m*\u001b[39m (\u001b[39m1\u001b[39m \u001b[39m+\u001b[39m epsilon) \u001b[39m-\u001b[39m debt)\n\u001b[0;32m--> 348\u001b[0m nodes \u001b[39m=\u001b[39m method(\n\u001b[1;32m 349\u001b[0m graph\u001b[39m.\u001b[39;49msubgraph(remaining_nodes),\n\u001b[1;32m 350\u001b[0m pop_col\u001b[39m=\u001b[39;49mpop_col,\n\u001b[1;32m 351\u001b[0m pop_target\u001b[39m=\u001b[39;49m(min_pop \u001b[39m+\u001b[39;49m max_pop) \u001b[39m/\u001b[39;49m \u001b[39m2\u001b[39;49m,\n\u001b[1;32m 352\u001b[0m epsilon\u001b[39m=\u001b[39;49m(max_pop \u001b[39m-\u001b[39;49m min_pop) \u001b[39m/\u001b[39;49m (\u001b[39m2\u001b[39;49m \u001b[39m*\u001b[39;49m pop_target),\n\u001b[1;32m 353\u001b[0m node_repeats\u001b[39m=\u001b[39;49mnode_repeats,\n\u001b[1;32m 354\u001b[0m )\n\u001b[1;32m 356\u001b[0m \u001b[39mif\u001b[39;00m nodes \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m 357\u001b[0m \u001b[39mraise\u001b[39;00m BalanceError()\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/gerrychain/tree.py:204\u001b[0m, in \u001b[0;36mbipartition_tree\u001b[0;34m(graph, pop_col, pop_target, epsilon, node_repeats, spanning_tree, spanning_tree_fn, balance_edge_fn, choice, max_attempts)\u001b[0m\n\u001b[1;32m 202\u001b[0m possible_cuts \u001b[39m=\u001b[39m []\n\u001b[1;32m 203\u001b[0m \u001b[39mif\u001b[39;00m spanning_tree \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m--> 204\u001b[0m spanning_tree \u001b[39m=\u001b[39m spanning_tree_fn(graph)\n\u001b[1;32m 206\u001b[0m restarts \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m\n\u001b[1;32m 207\u001b[0m attempts \u001b[39m=\u001b[39m \u001b[39m0\u001b[39m\n", + "File \u001b[0;32m~/Projects/GCRSEF2023/the-rebalancing-act/rba/util.py:68\u001b[0m, in \u001b[0;36mget_county_weighted_random_spanning_tree\u001b[0;34m(graph)\u001b[0m\n\u001b[1;32m 66\u001b[0m \u001b[39mfor\u001b[39;00m u, v \u001b[39min\u001b[39;00m graph\u001b[39m.\u001b[39medges:\n\u001b[1;32m 67\u001b[0m weight \u001b[39m=\u001b[39m random\u001b[39m.\u001b[39mrandom()\n\u001b[0;32m---> 68\u001b[0m \u001b[39mif\u001b[39;00m graph[u][\u001b[39m\"\u001b[39;49m\u001b[39mCOUNTYFP10\u001b[39;49m\u001b[39m\"\u001b[39;49m] \u001b[39m==\u001b[39m graph[v][\u001b[39m\"\u001b[39m\u001b[39mCOUNTYFP10\u001b[39m\u001b[39m\"\u001b[39m]:\n\u001b[1;32m 69\u001b[0m weight \u001b[39m*\u001b[39m\u001b[39m=\u001b[39m constants\u001b[39m.\u001b[39mSAME_COUNTY_PENALTY\n\u001b[1;32m 70\u001b[0m graph[u][v][\u001b[39m\"\u001b[39m\u001b[39mrandom_weight\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m weight\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/networkx/classes/coreviews.py:53\u001b[0m, in \u001b[0;36mAtlasView.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 52\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__getitem__\u001b[39m(\u001b[39mself\u001b[39m, key):\n\u001b[0;32m---> 53\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_atlas[key]\n", + "File \u001b[0;32m~/miniconda3/envs/rba/lib/python3.11/site-packages/networkx/classes/coreviews.py:286\u001b[0m, in \u001b[0;36mFilterAtlas.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 284\u001b[0m \u001b[39mif\u001b[39;00m key \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_atlas \u001b[39mand\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mNODE_OK(key):\n\u001b[1;32m 285\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_atlas[key]\n\u001b[0;32m--> 286\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mKeyError\u001b[39;00m(\u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mKey \u001b[39m\u001b[39m{\u001b[39;00mkey\u001b[39m}\u001b[39;00m\u001b[39m not found\u001b[39m\u001b[39m\"\u001b[39m)\n", + "\u001b[0;31mKeyError\u001b[0m: 'Key COUNTYFP10 not found'" + ] + } + ], + "source": [ + "data = pd.DataFrame(\n", + " [sorted(partition[\"district_gerry_scores\"]) + [partition[\"gerry_score\"]]\n", + " for partition in chain.with_progress_bar()],\n", + " columns=[f\"district{i}\" for i in range(1, NUM_DISTRICTS + 1)] + [\"state_gerry_score\"]\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(data[\"gerry_score\"])\n", + "plt.show()" + ] + } + ], + "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.11.0" + }, + "vscode": { + "interpreter": { + "hash": "3218fa316ded8cc4ffe9096b2baf1ef24e574309e5ed6f9629d67de0243fa942" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/vra_nh.json b/examples/vra_nh.json new file mode 100644 index 0000000..a982b73 --- /dev/null +++ b/examples/vra_nh.json @@ -0,0 +1,9 @@ +{ + "black": 0, + "hispanic": 0, + "asian": 0, + "native": 0, + "islander": 0, + "combined": 0, + "opportunity_threshold": 0.51 +} \ No newline at end of file diff --git a/rba/__main__.py b/rba/__main__.py index edc0525..0169c0a 100644 --- a/rba/__main__.py +++ b/rba/__main__.py @@ -7,7 +7,10 @@ Generates `num_thresholds` community maps based on the precinct graph `graph`, and writes to a file storing a list of individual communities, containing data on constituent precincts, birth and death times. - - districtgen + - districtgen [--graph_file] [--edge_lifetimes_file] [--vra_config] [--output_dir] + Runs simulated annealing algorithm, saves the ten best maps, as well as a dataframe keeping + track of various statistics for each state of the chain. `vra_config` is a JSON file + containing information about minority-opportunity district constraints. - ensemblegen - quantify - draw @@ -47,5 +50,11 @@ draw_parser.add_argument("--num_frames", type=int, default=50) draw_parser.set_defaults(func=rba.visualization.visualize) + optimize_parser = subparsers.add_parser("optimize") + optimize_parser.add_argument("--graph_file", type=str, default=os.path.join(package_dir, "data/2010/new_hampshire_geodata_merged.json")) + optimize_parser.add_argument("--edge_lifetime_file", type=str) + optimize_parser.add_argument("--vra_config", type=str) + optimize_parser.add_argument("--output_dir", type=str) + args = parser.parse_args() args.func(**{key: val for key, val in vars(args).items() if key != "func"}) \ No newline at end of file diff --git a/rba/constants.py b/rba/constants.py new file mode 100644 index 0000000..ceeaee0 --- /dev/null +++ b/rba/constants.py @@ -0,0 +1,11 @@ +"""Arbitrarily decided parameters for all algorithms. +""" + + +# Edges between nodes in the same county should be weighted less than those that cross, because +# the maximum spanning tree should be made more likely to choose an edge crossing county lines. +SAME_COUNTY_PENALTY = 0.5 + +POP_EQUALITY_THRESHOLD = 0.005 + +MINORITY_NAMES = ["black", "hispanic", "asian", "native", "islander"] \ No newline at end of file diff --git a/rba/district_optimization.py b/rba/district_optimization.py index e69de29..d45a25e 100644 --- a/rba/district_optimization.py +++ b/rba/district_optimization.py @@ -0,0 +1,64 @@ +""" +Given supercommunity edge lifetimes, uses simulated annealing to generate a map that minimizes +the average border edge lifetime while conforming to redistricting requirements. +""" + +from functools import partial +import random + +from gerrychain import (GeographicPartition, Partition, Graph, MarkovChain, + proposals, updaters, constraints, accept, Election) +from gerrychain.proposals import recom +from gerrychain.tree import bipartition_tree +from networkx import tree + +from .util import get_num_vra_districts, get_county_weighted_random_spanning_tree + + +class SimulatedAnnealingChain(MarkovChain): + """Augments gerrychain.MarkovChain to take both the current state and proposal in the `accept` + function. + """ + def __next__(self): + if self.counter == 0: + self.counter += 1 + return self.state + + while self.counter < self.total_steps: + proposed_next_state = self.proposal(self.state) + # Erase the parent of the parent, to avoid memory leak + if self.state is not None: + self.state.parent = None + + if self.is_valid(proposed_next_state): + if self.accept(self.state, proposed_next_state): + self.state = proposed_next_state + self.counter += 1 + return self.state + raise StopIteration + + +def accept_proposal(temperature, current_energy, proposed_energy): + """Simple simulated-annealing acceptance function. + """ + if current_energy > proposed_energy or random.random() < temperature: + return True + return False + + +def generate_districts_simulated_annealing(graph, edge_lifetimes, num_vra_districts, vra_threshold, + pop_equality_threshold): + """Returns the 10 best maps and a dataframe of statistics for the entire chain. + """ + + weighted_recom_proposal = partial( + recom, + method=partial( + bipartition_tree, + spanning_tree_fn=get_county_weighted_random_spanning_tree) + ) + + +def optimize(): + """ + """ \ No newline at end of file diff --git a/rba/util.py b/rba/util.py index 417c263..daf3dde 100644 --- a/rba/util.py +++ b/rba/util.py @@ -1,9 +1,12 @@ -from scipy.special import rel_entr -import networkx as nx +"""Miscellaneous utilities. +""" + +import random +import networkx as nx -def jenson_shannon_divergence(distribution1, distribution2): - average = [(distribution1[i] + distribution2[i])/2 for i in range(distribution1)] +from . import constants +from .district_quantification import quantify_gerrymandering def copy_adjacency(graph): @@ -14,4 +17,59 @@ def copy_adjacency(graph): copy_graph.add_node(node) for u, v in graph.edges: copy_graph.add_edge(u, v) - return copy_graph \ No newline at end of file + return copy_graph + + +def get_num_vra_districts(partition, label, threshold): + """Returns the number of minority-opportunity distrcts for a given minority and threshold. + + Parameters + ---------- + partition : gerrychain.Parition + Proposed district plan. + label : str + Node data key that returns the population of that minority. + threshold : float + Value between 0 and 1 indicating the percent population required for a district to be + considered minority opportunity. + """ + num_vra_districts = 0 + for part in partition.parts: + total_pop = 0 + minority_pop = 0 + for node in partition.parts[part]: + total_pop += partition.graph.nodes[node]["total_pop"] + if label == "total_combined": + for minority in constants.MINORITY_NAMES: + minority_pop += partition.graph.nodes[node][f"total_{minority}"] + else: + minority_pop += partition.graph.nodes[node][label] + if minority_pop / total_pop >= threshold: + num_vra_districts += 1 + return num_vra_districts + + +def get_gerrymandering_score(partition, edge_lifetimes): + """Returns the gerrymandering score of a partition. + """ + return quantify_gerrymandering(partition.graph, partition.subgraphs, edge_lifetimes)[1] + + +def get_district_gerrymandering_scores(partition, edge_lifetimes): + """Returns the gerrymandering scores of the districts in a partition""" + return quantify_gerrymandering(partition.graph, partition.subgraphs, edge_lifetimes)[0] + + +def get_county_weighted_random_spanning_tree(graph): + """Applies random edge weights to a graph, then multiplies those weights depending on whether or + not the edge crosses a county border. Then returns the maximum spanning tree for the graph.""" + for u, v in graph.edges: + weight = random.random() + if graph[u]["COUNTYFP10"] == graph[v]["COUNTYFP10"]: + weight *= constants.SAME_COUNTY_PENALTY + graph[u][v]["random_weight"] = weight + + spanning_tree = nx.tree.maximum_spanning_tree( + graph, algorithm="kruskal", weight="random_weight" + ) + return spanning_tree \ No newline at end of file diff --git a/rba/visualization.py b/rba/visualization.py index 5351e64..d16b3eb 100644 --- a/rba/visualization.py +++ b/rba/visualization.py @@ -10,6 +10,7 @@ import math from PIL import Image, ImageDraw, ImageFont +import geopandas import networkx as nx import shapely.geometry import shapely.ops @@ -104,6 +105,19 @@ def modify_coords(coords, bounds): return new_coords +def visualize_partition_geopandas(partition): + """Visualizes a gerrychain.Partition object using geopandas. + """ + data = {"assignment": [], "geometry": []} + for node in partition.graph: + data["assignment"].append(partition.assignment[node]) + data["geometry"].append(shapely.geometry.shape(partition.graph.nodes[node]['geometry'])) + + gdf = geopandas.GeoDataFrame(data) + del data + gdf.plot(column="assignment") + + def visualize_map(graph, output_fpath, node_coords, edge_coords, node_colors=None, edge_colors=None, edge_widths=None, node_list=None, additional_polygons=None, text=None, show=False): """Creates an image of a map and saves it to a file.