{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Max-P Regionalization for Multiple Components\n", "\n", "Authors: [Sergio Rey](https://github.com/sjsrey), [James Gaboardi](https://github.com/jGaboardi)\n", "\n", "\n", "The `max-p` problem involves the clustering of a set of geographic areas into the maximum number of homogeneous regions such that the value of a spatially extensive regional attribute is above a predefined threshold value. The spatially extensive attribute can be specified to ensure that each region contains sufficient population size, or a minimum number of enumeration units. The number of regions $p$ is endogenous to the problem and is useful for regionalization problems where the analyst does not require a fixed number of regions a-priori.\n", "\n", "Originally formulated as a mixed-integer problem in [Duque, Anselin, Rey (2012)](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9787.2011.00743.x), `max-p` is an [NP-hard problem](https://en.wikipedia.org/wiki/NP-hardness) and exact solutions are only feasible for small problem sizes. As such, a number of heuristic solution approaches have been suggested. PySAL implements the heuristic approach described in\n", "[Wei, Rey, and Knaap (2020)](https://www.tandfonline.com/doi/full/10.1080/13658816.2020.1759806).\n", "\n", "One issue with the current version of maxp is when it is applied to a collection of areas that have multiple connected components *and* some of the components do not allow for feasible region building.\n", "\n", "In this notebook we outline a number of alternatives to use when encountering infeasible components." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:18.700281Z", "start_time": "2022-10-29T17:23:18.651857Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:28.215532Z", "iopub.status.busy": "2025-04-06T22:39:28.215323Z", "iopub.status.idle": "2025-04-06T22:39:28.250500Z", "shell.execute_reply": "2025-04-06T22:39:28.250087Z", "shell.execute_reply.started": "2025-04-06T22:39:28.215507Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: 2025-04-06T18:39:28.239339-04:00\n", "\n", "Python implementation: CPython\n", "Python version : 3.12.9\n", "IPython version : 9.0.2\n", "\n", "Compiler : Clang 18.1.8 \n", "OS : Darwin\n", "Release : 24.4.0\n", "Machine : arm64\n", "Processor : arm\n", "CPU cores : 8\n", "Architecture: 64bit\n", "\n" ] } ], "source": [ "%config InlineBackend.figure_format = \"retina\"\n", "%load_ext watermark\n", "%watermark" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:20.967825Z", "start_time": "2022-10-29T17:23:18.704063Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:28.251319Z", "iopub.status.busy": "2025-04-06T22:39:28.251176Z", "iopub.status.idle": "2025-04-06T22:39:29.596273Z", "shell.execute_reply": "2025-04-06T22:39:29.596062Z", "shell.execute_reply.started": "2025-04-06T22:39:28.251301Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Watermark: 2.5.0\n", "\n", "matplotlib: 3.10.1\n", "spopt : 0.6.2.dev3+g13ca45e\n", "libpysal : 4.12.1\n", "shapely : 2.1.0\n", "scipy : 1.15.2\n", "numpy : 2.2.4\n", "geopandas : 1.0.1\n", "\n" ] } ], "source": [ "import warnings\n", "\n", "import geopandas\n", "import libpysal\n", "import matplotlib.pyplot as plt\n", "import numpy\n", "import scipy\n", "import shapely\n", "from scipy.spatial import KDTree\n", "from shapely.geometry import Polygon, box\n", "\n", "import spopt\n", "from spopt.region import MaxPHeuristic as MaxP\n", "\n", "plt.rcParams[\"figure.figsize\"] = [12, 8]\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "RANDOM_SEED = 123456\n", "\n", "%watermark -w\n", "%watermark -iv" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:20.979007Z", "start_time": "2022-10-29T17:23:20.972699Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.597507Z", "iopub.status.busy": "2025-04-06T22:39:29.597346Z", "iopub.status.idle": "2025-04-06T22:39:29.601067Z", "shell.execute_reply": "2025-04-06T22:39:29.600899Z", "shell.execute_reply.started": "2025-04-06T22:39:29.597500Z" } }, "outputs": [], "source": [ "n_cols = 5\n", "n_rows = 10\n", "b = 0\n", "h = w = 10\n", "component_0 = [box(l * w, b, l * w + w, b + h) for l in range(n_cols)] # noqa: E741\n", "b = b + h * 2\n", "component_1 = [\n", " box(l * w, b + h * r, l * w + w, b + h + h * r)\n", " for r in range(n_rows)\n", " for l in range(n_cols) # noqa: E741\n", "]\n", "geometries = component_0 + component_1" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:20.989502Z", "start_time": "2022-10-29T17:23:20.980637Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.601380Z", "iopub.status.busy": "2025-04-06T22:39:29.601322Z", "iopub.status.idle": "2025-04-06T22:39:29.603609Z", "shell.execute_reply": "2025-04-06T22:39:29.603462Z", "shell.execute_reply.started": "2025-04-06T22:39:29.601374Z" } }, "outputs": [ { "data": { "text/plain": [ "55" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(geometries)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.188034Z", "start_time": "2022-10-29T17:23:20.991878Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.603907Z", "iopub.status.busy": "2025-04-06T22:39:29.603855Z", "iopub.status.idle": "2025-04-06T22:39:29.668407Z", "shell.execute_reply": "2025-04-06T22:39:29.668203Z", "shell.execute_reply.started": "2025-04-06T22:39:29.603901Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAUnCAYAAAAy/TFcAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAewgAAHsIBbtB1PgAAU45JREFUeJzt3XucV3W96P/3wHBzhlS8/cJhe0ewXSceAeHGQttpWzQmSd26O4Xm1q4kltJVtOMlQRM7bk9uHnjrPAovbS8pSrnVwAscIjFtByWIxRBub1RyF2b9/qD5NsgwM8DM9ztveT4fDx6PxXzXWvOZ96C8WOv7/U5VURRFAACQSrdKLwAAgB0n4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACRUXekF7C7Wr18fzz33XERE7LffflFdbfQAsLvYtGlTvPLKKxER8e53vzt69+69y+dUEmXy3HPPxfDhwyu9DACgwubPnx/Dhg3b5fO4nQoAkJArcWWy3377lbbnz58f73znOyu4GgCgnFauXFm6I9e8CXaFiCuT5s+Be+c73xl1dXUVXA0AUCkd9bx4t1MBABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICEOjXiXn755XjggQdi0qRJceKJJ8a+++4bVVVVUVVVFWeddVa7zrF+/fq47777Yvz48fH+978/+vXrFz169Ih+/frF0UcfHZdeemmsXLmy3Wtau3ZtXH311TF8+PDo169f1NbWxuDBg+PCCy+MP/zhDzv5lQIAlFd1Z578gAMO2KXjn3322TjmmGPijTfe2OaxVatWxbx582LevHlx7bXXxvTp0+P0009v9XxLly6Nk046KX77299u9fHFixfH4sWLY/r06fGjH/0oRo8evUvrBgDobGW7nTpgwIA44YQTduiYv/zlL6WAGzlyZHznO9+Jhx9+OJ5++un46U9/Gp/5zGeie/fu8cYbb8S//Mu/xEMPPbTdc61evTpOPvnkUsCde+658cgjj8RTTz0VV1xxRdTW1saf//znOO200+LZZ5/d+S8UAKAMOvVK3KRJk2LYsGExbNiwOOCAA+LFF1+MQw45pN3Hd+vWLU4//fS45JJL4qijjtrm8RNOOCFOPPHEOOWUU2Lz5s0xfvz4eP7556Oqqmqbfa+55ppYvHhxRERMmTIlLrrootJjRx99dBx33HHxwQ9+MNauXRsTJkyIRx99dCe+YgCA8qgqiqIo1ydrHnHjxo2LW2+9tUPOe+qpp8Z//Md/RETE008/HUOGDNnq8TfffDP233//+NOf/hSDBw+OX//619Gt27YXIT/72c/Gv//7v0dExIIFC+J973tfh6wvIqKhoSEGDBgQERHLly+Purq6Djs3ANC1dUYHvC1enXrccceVtpcuXbrN4z//+c/jT3/6U0RsiceWAi4itnqxxd13392hawQA6Ehvi4jbsGFDabulQHv88cdL26NGjdrueYYOHRo1NTUREfHEE0904AoBADrW2yLiZs+eXdoeNGjQNo8vWrSo1cebVFdXx2GHHbbNMQAAXU2nvrChHH71q1/FzJkzIyLiXe96V4svgFi+fHlERNTU1MRee+3V6vkGDBgQzz77bLzyyiuxYcOG6NWrV7vW0dDQ0OrjO/JedgAAbUkdcRs2bIh//dd/jc2bN0dExJVXXtnifk1vU1JbW9vmOZtup0ZseVuS9kZc05MVK+muu+6KSZMmtfi+enSOdevWxerVq6O2tjb69OlT6eXsFsy8Msy9/My8Mvr27RuXXXZZnHrqqZVeSptSR9wXv/jFWLBgQURsecHCmDFjWtxv/fr1ERHRs2fPNs/ZPNrWrVvXAassn0mTJpXeRoXyev311yu9hN2OmVeGuZefmZffxRdfLOI603e+852YPn16RES8733vixtuuGG7+/bu3TsiIjZu3NjmeZu/SGJH/uXTdMt2e1auXBnDhw9v9/l2RukKXFW36F6zd6d+LrbYvPq1LRtmXjZmXhnmXn5mXn6b16yKKBrT3NFKGXH//u//Ht/4xjciIuLII4+Mhx56aKvboG/Vt2/fiNhye7Qta9asKW235/Zrk670vm/da/aOui/cVull7BZ+P2VMRNFo5mVk5pVh7uVn5uXXcMO4v8VzAulenTpjxoz4/Oc/HxERBx10UPznf/5n7Lfffq0e0xRYa9asKb1f3PY0XVHbb7/92v18OACAcksVcT/5yU/iU5/6VDQ2NsY73/nOeOSRR9p1Baz5K1Zbe87Ypk2bSm8WPHjw4F1fMABAJ0kTcY888kicfvrpsWnTpthnn33i4YcfLr2nW1uOOeaY0nbz95R7qwULFpRup44cOXLXFgwA0IlSRNxTTz0V9fX1sWHDhnjHO94RP/3pT+Nd73pXu48/9thjY88994yIiNtuuy229+Nim/8s11NOOWWX1gwA0Jm6fMQ988wzcdJJJ8WaNWuipqYmHnzwwR3+wfQ9e/aML33pSxGx5ScxXHPNNdvsM3fu3LjpppsiYsuP5ho2bNiuLx4AoJN06qtTn3jiiViyZEnp96+++mppe8mSJVtd+YrY+gfQR2z5YfYf+chHSi9GuPzyy2PPPfeMX//619v9nPvvv3/sv//+23z8oosuijvuuCN+97vfxcSJE2PJkiVxxhlnRJ8+feKxxx6LK6+8MjZt2hR9+vSJ6667boe/VgCAcurUiJs+fXrcdlvLL4t+8skn48knn9zqY2+NuMcffzxefvnl0u8vuOCCNj/nJZdcEpdeeuk2H+/bt2/MnDkzRo8eHc8//3xMmzYtpk2bttU+73jHO+KHP/xhvPe9723z8wAAVFKXv53akQ4//PBYuHBhTJ48OYYOHRp77bVX7LHHHnHkkUfGBRdcEM8++2ycfPLJlV4mAECbOvVK3K233rrNLdMdcdZZZ21zdW5X1dTUxMSJE2PixIkdel4AgHLara7EAQC8XYg4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABKqKoqiqPQidgcNDQ0xYMCAiIhYvnx51NXVdfjn2GeffeL111+PqOoW3Wv27vDzs63Nq1/bsmHmZWPmlWHu5Wfm5dc08379+sVrr73WoefujA6o3uUz0GWsXr16y0bR+Lf/+CkPMy8/M68Mcy8/My+76uoceZRjlbRLbW3tlitxERFV7pSXRdH4t20zLw8zrwxzLz8zL7+/zrxHjx4VXkj7iLi3kT59+kRERPfafaLuC7dVeDW7h99PGRNRNJp5GZl5ZZh7+Zl5+TXcMC7VVU9pDwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQUKdG3MsvvxwPPPBATJo0KU488cTYd999o6qqKqqqquKss87a4fPNmjUrxo4dG3V1ddGrV6+oq6uLsWPHxqxZs9p9jrVr18bVV18dw4cPj379+kVtbW0MHjw4LrzwwvjDH/6ww2sCAKiE6s48+QEHHNAh5ymKIj772c/GtGnTtvr4ihUr4p577ol77rknzjvvvLjxxhujqqpqu+dZunRpnHTSSfHb3/52q48vXrw4Fi9eHNOnT48f/ehHMXr06A5ZNwBAZynb7dQBAwbECSecsFPHfutb3yoF3JAhQ2LGjBkxf/78mDFjRgwZMiQiIqZNmxYXX3zxds+xevXqOPnkk0sBd+6558YjjzwSTz31VFxxxRVRW1sbf/7zn+O0006LZ599dqfWCQBQLp16JW7SpEkxbNiwGDZsWBxwwAHx4osvxiGHHLJD51iyZElMmTIlIiKGDh0ac+bMiT59+kRExLBhw2LMmDExatSoWLBgQUyePDnOPvvsOOyww7Y5zzXXXBOLFy+OiIgpU6bERRddVHrs6KOPjuOOOy4++MEPxtq1a2PChAnx6KOP7uyXDQDQ6Tr1Sty3v/3tOPnkk3fpturUqVNj06ZNERFx/fXXlwKuyR577BHXX399RERs2rQprrvuum3O8eabb8b3vve9iIgYPHhwfOUrX9lmn6OPPjrOOeeciIh47LHH4pe//OVOrxkAoLN16VenFkUR9913X0REDBo0KEaMGNHifiNGjIgjjzwyIiLuvffeKIpiq8d//vOfx5/+9KeIiBg3blx069byl938xRZ33333Lq4eAKDzdOmIW7ZsWaxYsSIiIkaNGtXqvk2PNzQ0xIsvvrjVY48//vg2+7Vk6NChUVNTExERTzzxxM4sGQCgLDr1OXG7atGiRaXtQYMGtbpv88cXLVq01XPv2nue6urqOOyww+LZZ5/d6pj2aGhoaPXxlStX7tD5AABa06Ujbvny5aXturq6VvcdMGBAi8c1/31NTU3stddebZ7n2WefjVdeeSU2bNgQvXr1atdam3/+Slm3bl1ERGxesyoabhhX4dXsJorGiDDzsjLzyjD38jPzstu8+rWI+Nvfp11dl464N954o7RdW1vb6r5Nt0EjtrydSEvnaescLZ2nvRHXFZS+7qKx9AeRMjHz8jPzyjD38jPzsquu7tJ5VNKlV7l+/frSds+ePVvdt3lsvbWgm87T1jnaOk9r3nr1761WrlwZw4cPb/f5dkZtbW28/vrrW35T1aWf7vj28dd/KUeEmZeLmVeGuZefmZffX2feo0ePCi+kfbp0xPXu3bu0vXHjxlb33bBhQ2n7rW9D0nSets7R1nla09bt3nJoWm/32n2i7gu3VXg1u4ffTxkTUTSaeRmZeWWYe/mZefk13DAu1VXPLp32ffv2LW2/9RbpW61Zs6a0/dbbpk3naescbZ0HAKCr6NIR1/zqVluv/mx+O/OtLzJoOs+aNWtK7xfX1nn222+/VM+HAwB2L1064o466qjSdtOPzNqe5o8PHjx4p86zadOmWLp0aYvnAADoSrp0xB1yyCHRv3//iIiYPXt2q/vOmTMnIiIOPPDAOPjgg7d67Jhjjiltt3aeBQsWlG6njhw5cmeWDABQFl064qqqqqK+vj4itlxBmzdvXov7zZs3r3SFrb6+PqqqqrZ6/Nhjj40999wzIiJuu+22bX4sV5Nbb721tH3KKafs6vIBADpNl464iIgJEyaU3q9l/Pjx27ztx7p162L8+PERseV9XSZMmLDNOXr27Blf+tKXImLLT2+45pprttln7ty5cdNNN0XElh/NNWzYsI78MgAAOlSnvsXIE088EUuWLCn9/tVXXy1tL1myZKsrXxFb/wD6JgMHDowLL7wwrrrqqliwYEGMHDkyvvrVr8Zhhx0WS5cujcmTJ8fChQsjIuKiiy6KI444osW1XHTRRXHHHXfE7373u5g4cWIsWbIkzjjjjOjTp0889thjceWVV8amTZuiT58+cd111+3y1w4A0Jk6NeKmT58et93W8nvbPPnkk/Hkk09u9bGWIi4i4oorroiXX345br755li4cGGcccYZ2+xzzjnnxOWXX77dtfTt2zdmzpwZo0ePjueffz6mTZsW06ZN22qfd7zjHfHDH/4w3vve97b+hQEAVFiXv50aEdGtW7e46aabYubMmVFfXx/9+/ePnj17Rv/+/aO+vj4efPDBmD59enTr1vqXc/jhh8fChQtj8uTJMXTo0Nhrr71ijz32iCOPPDIuuOCCePbZZ+Pkk08u01cFALDzOvVK3K233rrNLdNdMXr06Bg9evQunaOmpiYmTpwYEydO7KBVAQCUX4orcQAAbE3EAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJBQVVEURaUXsTtoaGiIAQMGRETE8uXLo66ursM/xz777BOvv/56RFW36F6zd4efn21tXv3alg0zLxszrwxzLz8zL7+mmffr1y9ee+21Dj13Z3RA9S6fgS5j9erVWzaKxr/9x095mHn5mXllmHv5mXnZVVfnyKMcq6Rdamtrt1yJi4iocqe8LIrGv22beXmYeWWYe/mZefn9deY9evSo8ELaR8S9jfTp0yciIrrX7hN1X7itwqvZPfx+ypiIotHMy8jMK8Pcy8/My6/hhnGprnpKewCAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmlibiNGzfGTTfdFP/0T/8U73znO6NXr15RW1sbRx55ZHz605+OefPmtes8s2bNirFjx0ZdXV306tUr6urqYuzYsTFr1qxO/goAADpOdaUX0B7Lly+Pk046KZ577rmtPr5x48b43e9+F7/73e/illtuiQsuuCC++93vRlVV1TbnKIoiPvvZz8a0adO2+viKFSvinnvuiXvuuSfOO++8uPHGG1s8HgCgK+nyV+I2bdq0VcC95z3viVtvvTXmzp0bP/vZz2LSpElRU1MTERFTp06Na665psXzfOtb3yoF3JAhQ2LGjBkxf/78mDFjRgwZMiQiIqZNmxYXX3xxGb4qAIBd0+WvxN13332lgDv66KPj8ccfj+7du5ceP/7442PMmDFx9NFHx5tvvhnf+c534oILLojq6r99aUuWLIkpU6ZERMTQoUNjzpw50adPn4iIGDZsWIwZMyZGjRoVCxYsiMmTJ8fZZ58dhx12WBm/SgCAHdPlr8Q9+eSTpe2vf/3rWwVck/e9731x8sknR0TEqlWrYvHixVs9PnXq1Ni0aVNERFx//fWlgGuyxx57xPXXXx8RW678XXfddR35JQAAdLguH3EbN24sbR966KHb3a/5lbMNGzaUtouiiPvuuy8iIgYNGhQjRoxo8fgRI0bEkUceGRER9957bxRFsUvrBgDoTF0+4gYOHFjafuGFF7a739KlSyMioqqqKo444ojSx5ctWxYrVqyIiIhRo0a1+rmaHm9oaIgXX3xxZ5cMANDpunzEnXnmmfGOd7wjIiImT54cmzdv3mafhQsXxsyZMyMi4owzzijtHxGxaNGi0vagQYNa/VzNH29+HABAV9PlX9iw3377xa233hqf+MQn4sknn4xhw4bFhAkTYuDAgbF69ep48skn47vf/W5s3Lgx3vve98a111671fHLly8vbdfV1bX6uQYMGNDice3R0NDQ6uMrV67cofMBALSmy0dcRMQpp5wSCxYsiGuvvTZuvvnmGDdu3FaPH3DAAfHtb387zjvvvNLbjTR54403Stu1tbWtfp7mx65evXqH1tg8ACtl3bp1ERGxec2qaLhhXBt70yGKxogw87Iy88ow9/Iz87LbvPq1iPjb36ddXYqIe/PNN+NHP/pR3H///S2+4OC///u/Y8aMGTFw4MA46aSTtnps/fr1pe2ePXu2+nl69epV2s7yDWyuFJ5FY+kPImVi5uVn5pVh7uVn5mXX/G3KurIuv8o1a9bE6NGjY86cOdG9e/eYOHFinH322XHooYfG+vXr4//9v/8X/+t//a944okn4qMf/WhMnTo1zj///NLxvXv3Lm03f6VrS5q/qvWtb0PSlrZuv65cuTKGDx++Q+fcUbW1tfH6669v+U1Vl3+649vDX/+lHBFmXi5mXhnmXn5mXn5/nXmPHj0qvJD26fIRd8kll8ScOXMiIuKmm27a6lZqz5494/jjj4/jjjsuTjjhhHjsscfiy1/+chx33HHxnve8JyIi+vbtW9q/rVuka9asKW23dev1rdp6vl05NIVn99p9ou4Lt1V4NbuH308ZE1E0mnkZmXllmHv5mXn5NdwwLtVVzy6d9kVRxC233BIRW95q5K3PhWtSXV0dl112WURENDY2lo6J2Dqu2nrxQfOraV3hOW4AANvTpSPuv//7v0u3B5t+vun2vO997yttN/+JDUcddVSLH29J88cHDx68Q2sFACinLh1xzZ9Y2PRjs7bnzTffbPG4Qw45JPr37x8REbNnz271HE23bQ888MA4+OCDd3S5AABl06Ujrl+/fqU37p07d26rIdc80A455JDSdlVVVdTX10fElitt8+bNa/H4efPmla7E1dfXR1VV1S6vHwCgs3TpiOvWrVvpLUP++Mc/xhVXXNHifqtWrYqvfvWrpd+ffPLJWz0+YcKE0tW58ePHb/P2IevWrYvx48dHxJareBMmTOioLwEAoFN06YiLiJg0aVLsscceERFx6aWXxpgxY+I//uM/YuHChTF37tyYOnVqvPe9743f/OY3ERHxj//4j3HCCSdsdY6BAwfGhRdeGBERCxYsiJEjR8Ydd9wRCxYsiDvuuCNGjhwZCxYsiIiIiy66aKufvQoA0BV1+bcYGTRoUNx3331x5plnxquvvhr3339/3H///S3u+6EPfSjuuuuuFh+74oor4uWXX46bb745Fi5cGGecccY2+5xzzjlx+eWXd+j6AQA6Q5e/EhcR8eEPfzgWL14ckydPjmOPPTb222+/6NGjR/Tp0ycOOeSQOP300+Pee++N//zP/4y99967xXN069Ytbrrpppg5c2bU19dH//79o2fPntG/f/+or6+PBx98MKZPnx7duqUYCQCwm+vyV+Ka7LPPPjFx4sSYOHHiLp1n9OjRMXr06A5aFQBAZbjsBACQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AIKGqoiiKSi9id9DQ0BADBgyIiIjly5dHXV1dh3+OffbZJ15//fWIqm7RvWbvDj8/29q8+rUtG2ZeNmZeGeZefmZefk0z79evX7z22msdeu7O6IDqXT4DXcbq1au3bBSNf/uPn/Iw8/Iz88ow9/Iz87Krrs6RRzlWSbvU1tZuuRIXEVHlTnlZFI1/2zbz8jDzyjD38jPz8vvrzHv06FHhhbSPiHsb6dOnT0REdK/dJ+q+cFuFV7N7+P2UMRFFo5mXkZlXhrmXn5mXX8MN41Jd9ZT2AAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJpYu4V199NaZMmRIjR46M/+//+/+iV69e0b9//3j/+98fF110UcydO7fNc8yaNSvGjh0bdXV10atXr6irq4uxY8fGrFmzyvAVAADsuupKL2BH3HXXXfG5z30uXnvtta0+vnLlyli5cmXMnz8/nn/++bj33ntbPL4oivjsZz8b06ZN2+rjK1asiHvuuSfuueeeOO+88+LGG2+MqqqqzvoyAAB2WZqI+8EPfhBnn312NDY2xv777x+f+9zn4phjjol+/frFSy+9FEuXLo37778/evTosd1zfOtb3yoF3JAhQ2LixIlx2GGHxdKlS2PKlCmxcOHCmDZtWuy3335x+eWXl+tLAwDYYSkibtGiRXHeeedFY2NjfOADH4j7778/9txzz232Gz9+fGzcuLHFcyxZsiSmTJkSERFDhw6NOXPmRJ8+fSIiYtiwYTFmzJgYNWpULFiwICZPnhxnn312HHbYYZ33RQEA7IIUz4kbP358bNiwIfbdd9+4++67Wwy4Jj179mzx41OnTo1NmzZFRMT1119fCrgme+yxR1x//fUREbFp06a47rrrOmbxAACdoMtH3OLFi+ORRx6JiIgvfvGLse++++7wOYqiiPvuuy8iIgYNGhQjRoxocb8RI0bEkUceGRER9957bxRFsZOrBgDoXF0+4u66667S9mmnnVbaXrVqVTz//PPbvMihJcuWLYsVK1ZERMSoUaNa3bfp8YaGhnjxxRd3YsUAAJ2vy0fcvHnzIiJizz33jMGDB8cPf/jD+B//439Ev379YuDAgbHvvvvGoYceGt/+9rdj9erVLZ5j0aJFpe1Bgwa1+vmaP978OACArqTLv7DhN7/5TUREHHzwwTF+/Pi44YYbttln2bJlcemll8aPf/zj+OlPfxr9+/ff6vHly5eXtuvq6lr9fAMGDGjxuLY0NDS0+vjKlSvbfS4AgLZ0+Yh7/fXXI2LLc+N+9atfxV577RVXXXVVjB07Nt7xjnfEc889F5MmTYqHHnoofv3rX8dpp50Wjz/+eHTr9reLjG+88UZpu7a2ttXPV1NTU9re3pW9ljSPv0pZt25dRERsXrMqGm4YV+HV7CaKxogw87Iy88ow9/Iz87LbvHrLU7Sa/j7t6rp8xK1ZsyYiIjZs2BDdu3ePhx56aKsXJgwdOjQeeOCBOPnkk+Ohhx6Kp556Ku6+++449dRTS/usX7++tL29V6826dWrV2k7yzexSSk6i8bSH0TKxMzLz8wrw9zLz8zLrrq6y+dRRCSIuN69e5dC7rTTTmvxlaXdunWLq6++Oh566KGIiJgxY8ZWEde7d+/S9vbeR67Jhg0bSttvfRuS1rR163XlypUxfPjwdp9vZ9TW1pauXEZVl3+649vDX/+lHBFmXi5mXhnmXn5mXn5/nXlrPzigK+nyEde3b99SxJ144onb3e9d73pXHHjggbFixYr4xS9+sc05mrR1i7Tpc0W0feu1ubaea1cOTdHZvXafqPvCbRVeze7h91PGRBSNZl5GZl4Z5l5+Zl5+DTeMS3XVs8unffPnmrX3RQkvv/zyVh9vflxbL0BofkWtKzzPDQCgJV0+4t71rneVtjdv3tzqvk2Pv/Ve9lFHHVXaXrx4cavnaP744MGD271OAIBy6vIR98EPfrC0vXTp0lb3feGFFyIi4sADD9zq44ccckjpbUdmz57d6jnmzJlTOsfBBx+8o8sFACiLLh9xY8aMKT3B8O67797ufrNnzy799IYPfOADWz1WVVUV9fX1EbHlSlvTGwi/1bx580pX4urr66OqqmqX1w8A0Bm6fMTts88+8a//+q8REfHwww/H7bffvs0+b7zxRkyYMKH0+8985jPb7DNhwoTSbdbx48dv8/Yh69ati/Hjx0fEltuxzc8HANDVdPmIi4j49re/HX/3d38XERGf/OQnY/z48fHYY4/FL3/5y7j11ltj+PDh8cwzz0RExOc+97kYNmzYNucYOHBgXHjhhRERsWDBghg5cmTccccdsWDBgrjjjjti5MiRsWDBgoiIuOiii+KII44ozxcHALATuvxbjERE7LfffjFr1qwYM2ZMLFmyJP7t3/4t/u3f/m2b/T796U/H9773ve2e54orroiXX345br755li4cGGcccYZ2+xzzjnnxOWXX96h6wcA6GgprsRFbHml6DPPPBNXX311vP/9749+/fpFz549o66uLv75n/85Hn300bjppptafYO+bt26xU033RQzZ86M+vr66N+/f/Ts2TP69+8f9fX18eCDD8b06dO3+pFdAABdUYorcU1qamriwgsvLN0W3VmjR4+O0aNHd9CqAADKzyUnAICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACVUVRVFUehG7g4aGhhgwYEBERCxfvjzq6uo6/HPss88+8frrr0dUdYvuNXt3+PnZ1ubVr23ZMPOyMfPKMPfyM/Pya5p5v3794rXXXuvQc3dGB1Tv8hnoMlavXr1lo2j823/8lIeZl5+ZV4a5l5+Zl111dY48yrFK2qW2tnbLlbiIiCp3ysuiaPzbtpmXh5lXhrmXn5mX319n3qNHjwovpH1E3NtInz59IiKie+0+UfeF2yq8mt3D76eMiSgazbyMzLwyzL38zLz8Gm4Yl+qqp7QHAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEgodcRNnDgxqqqqSr9+/vOft3nMrFmzYuzYsVFXVxe9evWKurq6GDt2bMyaNavzFwwA0EHSRtyvfvWrmDp1arv3L4oiPvOZz8SJJ54Y99xzT6xYsSI2btwYK1asiHvuuSdOPPHE+MxnPhNFUXTiqgEAOkbKiGtsbIxzzz03Nm3aFPvvv3+7jvnWt74V06ZNi4iIIUOGxIwZM2L+/PkxY8aMGDJkSERETJs2LS6++OJOWzcAQEdJGXH/+3//7/jFL34RgwYNinPOOafN/ZcsWRJTpkyJiIihQ4fGk08+GWeccUYMGzYszjjjjHjiiSdi6NChERExefLkWLp0aaeuHwBgV6WLuOXLl5euln3/+9+Pnj17tnnM1KlTY9OmTRERcf3110efPn22enyPPfaI66+/PiIiNm3aFNddd13HLhoAoIOli7jPf/7zsXr16hg3blwce+yxbe5fFEXcd999ERExaNCgGDFiRIv7jRgxIo488siIiLj33ns9Nw4A6NJSRdydd94ZDzzwQPTr1y+uvvrqdh2zbNmyWLFiRUREjBo1qtV9mx5vaGiIF198cZfWCgDQmaorvYD2+tOf/hTnn39+RGx53tp+++3XruMWLVpU2h40aFCr+zZ/fNGiRXHIIYe0e30NDQ2tPr5y5cp2nwsAoC1pIm7ixInx0ksvxT/8wz+068UMTZYvX17arqura3XfAQMGtHhcezQ/tlLWrVsXERGb16yKhhvGVXg1u4miMSLMvKzMvDLMvfzMvOw2r34tIv7292lXlyLinnjiiZg+fXpUV1fHjTfeGFVVVe0+9o033iht19bWtrpvTU1NaXv16tU7vtAKK625aCz9QaRMzLz8zLwyzL38zLzsqqtT5FHXj7iNGzfGeeedF0VRxAUXXBDvfve7d+j49evXl7bbeiVrr169Sts7WuFtXblbuXJlDB8+fIfOuaNqa2vj9ddf3/KbqlRPd8zrr/9SjggzLxczrwxzLz8zL7+/zrxHjx4VXkj7dPmIu/LKK2PRokXxd3/3d3HJJZfs8PG9e/cubW/cuLHVfTds2FDafuvbkLSlrVu15dC05u61+0TdF26r8Gp2D7+fMiaiaDTzMjLzyjD38jPz8mu4YVyqq55dOu0XL14c3/nOdyJiy/u7Nb/d2V59+/Ytbbd1i3TNmjWl7bZuvQIAVFKXvhI3derU2LhxYxx66KGxdu3auP3227fZ59e//nVp+9FHH42XXnopIiI++tGPRk1NzVZXyNp6BWnzW6Jd4YUKAADb06Ujrun25gsvvBBnnnlmm/tfdtllpe1ly5ZFTU1NHHXUUaWPLV68uNXjmz8+ePDgHV0uAEDZdOnbqR3hkEMOif79+0dExOzZs1vdd86cORERceCBB8bBBx/c2UsDANhpXTribr311iiKotVfzV/s8Nhjj5U+3hRhVVVVUV9fHxFbrrTNmzevxc81b9680pW4+vr6HXobEwCAcuvSEddRJkyYUHrPl/Hjx2/z9iHr1q2L8ePHR8SW94aZMGFCuZcIALBDdouIGzhwYFx44YUREbFgwYIYOXJk3HHHHbFgwYK44447YuTIkbFgwYKIiLjoooviiCOOqORyAQDa1KVf2NCRrrjiinj55Zfj5ptvjoULF8YZZ5yxzT7nnHNOXH755RVYHQDAjtktrsRFRHTr1i1uuummmDlzZtTX10f//v2jZ8+e0b9//6ivr48HH3wwpk+fHt267TYjAQASS38l7tJLL41LL7203fuPHj06Ro8e3XkLAgAoA5edAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJFRVFEVR6UXsDhoaGmLAgAEREbF8+fKoq6vr8M+xzz77xOuvvx5R1S261+zd4ednW5tXv7Zlw8zLxswrw9zLz8zLr2nm/fr1i9dee61Dz90ZHVC9y2egy1i9evWWjaLxb//xUx5mXn5mXhnmXn5mXnbV1TnyKMcqaZfa2totV+IiIqrcKS+LovFv22ZeHmZeGeZefmZefn+deY8ePSq8kPYRcW8jffr0iYiI7rX7RN0XbqvwanYPv58yJqJoNPMyMvPKMPfyM/Pya7hhXKqrntIeACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiki7umnn44rr7wyTjzxxBgwYED06tUramtrY+DAgXHWWWfF448/vkPnmzVrVowdOzbq6uqiV69eUVdXF2PHjo1Zs2Z10lcAANCxqiu9gLaMGjUq5syZs83HN27cGM8//3w8//zzcdttt8UnP/nJmD59evTs2XO75yqKIj772c/GtGnTtvr4ihUr4p577ol77rknzjvvvLjxxhujqqqqw78WAICO0uWvxK1YsSIiIvr37x/nn39+/PjHP4758+fH3Llz49prr40DDzwwIiL+7//9v3HWWWe1eq5vfetbpYAbMmRIzJgxI+bPnx8zZsyIIUOGRETEtGnT4uKLL+68LwgAoAN0+StxgwYNiiuvvDI+/vGPR/fu3bd6bMSIEfHJT34yRo4cGb/73e9ixowZ8bnPfS4+8IEPbHOeJUuWxJQpUyIiYujQoTFnzpzo06dPREQMGzYsxowZE6NGjYoFCxbE5MmT4+yzz47DDjus879AAICd0OWvxD3wwANx+umnbxNwTfbdd9/47ne/W/r9j3/84xb3mzp1amzatCkiIq6//vpSwDXZY4894vrrr4+IiE2bNsV1113XAasHAOgcXT7i2uPYY48tbS9dunSbx4uiiPvuuy8itlzZGzFiRIvnGTFiRBx55JEREXHvvfdGURQdv1gAgA7wtoi4jRs3lra7ddv2S1q2bFnpuXWjRo1q9VxNjzc0NMSLL77YcYsEAOhAb4uImz17dml70KBB2zy+aNGiVh9vrvnjzY8DAOhKuvwLG9rS2NgYV111Ven3p59++jb7LF++vLRdV1fX6vkGDBjQ4nFtaWhoaPXxlStXtvtcAABtSR9xU6dOjfnz50dExCmnnBJDhw7dZp833nijtF1bW9vq+Wpqakrbq1evbvc6msdfpW1esyoabhhX6WXsHorGiDDzsjLzyjD38jPzstu8+rVKL2GHpI642bNnx9e+9rWIiNh///3j+9//fov7rV+/vrTd2psBR0T06tWrtL1u3boOWGX59O3bd8tG0ZjuD2J6Zl5+Zl4Z5l5+Zl52X/7ylyu9hHZJG3H/9V//Faecckps2rQpevXqFXfeeWcccMABLe7bu3fv0nbzF0G0ZMOGDaXtt74NSWvauvW6cuXKGD58eLvPtzMuu+yyuPjii7e68kjnWrduXVRXV0ePHj0qvZTdhplXhrmXn5lXxpe//GUR15mWLVsWJ5xwQqxatSq6d+8eM2bMaPVVp6UrVNH2LdI1a9aUttu69dpcW8+1K4dTTz01Tj311EovAwAog3SvTv3jH/8YH/7wh+OPf/xjVFVVxc033xynnHJKq8c0D6y2XoDQ/IpaV3qeGwBAc6ki7tVXX43jjz8+XnjhhYjY8pMXPvWpT7V53FFHHVXaXrx4cav7Nn988ODBO7lSAIDOlSbi/vznP8dHPvKR+M1vfhMREVdddVV84QtfaNexhxxySPTv3z8itn5PuZbMmTMnIiIOPPDAOPjgg3d+wQAAnShFxK1duzZOOumkePrppyMi4pvf/GZ89atfbffxVVVVUV9fHxFbrrTNmzevxf3mzZtXuhJXX18fVVVVu7hyAIDO0eUjbuPGjXHKKafEk08+GRER559/flx++eU7fJ4JEyZEdfWW13GMHz9+m7cPWbduXYwfPz4iIqqrq2PChAm7tnAAgE7U5V+deuaZZ8bPfvaziIj40Ic+FOecc078+te/3u7+PXv2jIEDB27z8YEDB8aFF14YV111VSxYsCBGjhwZX/3qV+Owww6LpUuXxuTJk2PhwoUREXHRRRfFEUcc0TlfEABAB6gqiqKo9CJas6O3NA866KDt/uD6xsbGOPfcc+Pmm2/e7vHnnHNOTJs2Lbp169iLlA0NDaVXuy5fvrxLvCUJAFAendEBXf52akfq1q1b3HTTTTFz5syor6+P/v37R8+ePaN///5RX18fDz74YEyfPr3DAw4AoKN1+dupnXGhcPTo0TF69OgOPy8AQLm45AQAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACCh6kovgI5z1113xaRJk+KNN96o9FJ2G+vWrYvVq1dHbW1t9OnTp9LL2S2YeWWYe/mZeWX07ds3Lrvssjj11FMrvZQ2VRVFUVR6EbuDhoaGGDBgQERELF++POrq6jr8cwwePDgWL17c4ecFgN3JoEGDYtGiRR16zs7oAFfi3kZKV+CqukX3mr0ru5jdxObVr23ZMPOyMfPKMPfyM/Py27xmVUTRmOaOloh7G+pes3fUfeG2Si9jt/D7KWMiikYzLyMzrwxzLz8zL7+GG8b9LZ4T8MIGAICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICEdsuI+8Mf/hAXXnhhDB48OGpqaqJfv34xfPjwuOaaa2Lt2rWVXh4AQJuqK72Acps5c2Z84hOfiD//+c+lj61duzZ+8YtfxC9+8YuYPn16PPjgg3HooYdWcJUAAK3bra7E/epXv4rTTz89/vznP0dtbW1cccUV8dRTT8UjjzwS5557bkRE/Pa3v42TTjopVq9eXeHVAgBs3251JW7ChAmxdu3aqK6ujp/97Gdx9NFHlx770Ic+FEcccURMnDgxFi9eHNdee21MmjSpgqsFANi+3eZK3C9+8Yv4+c9/HhER55xzzlYB1+QrX/lKDB48OCIirrvuunjzzTfLuUQAgHbbbSLu3nvvLW2fffbZLe7TrVu3+NSnPhUREatWrSpFHwBAV7PbRNzjjz8eERE1NTXxvve9b7v7jRo1qrT9xBNPdPq6AAB2xm4TcYsWLYqIiMMPPzyqq7f/VMBBgwZtcwwAQFezW7ywYf369fHqq69GRERdXV2r++69995RU1MTa9asieXLl7f7czQ0NLT6+MqVK9t9LgCAtuwWEffGG2+Utmtra9vcvyniduRtRgYMGLBTa+sMm9esioYbxlV6GbuHojEizLyszLwyzL38zLzsNq9+rdJL2CG7RcStX7++tN2zZ8829+/Vq1dERKxbt67T1tQZ+vbtu2WjaEz3BzE9My8/M68Mcy8/My+7L3/5y5VeQrvsFhHXu3fv0vbGjRvb3H/Dhg0REdGnT592f462br2uXLkyhg8f3u7z7YzLLrssLr744q2uPNK51q1bF9XV1dGjR49KL2W3YeaVYe7lZ+aV8eUvf1nEdSWlK1QR7bpFumbNmoho363XJm09164cTj311Dj11FMrvQwAoAx2i1en9u7dO/bdd9+IaPsFCKtWrSpFXFd6nhsAQHO7RcRFROknMSxZsiQ2bdq03f0WL168zTEAAF3NbhNxxxxzTERsuVX6y1/+crv7zZ49u7Q9cuTITl8XAMDO2G0i7mMf+1hp+5Zbbmlxn8bGxvjBD34QERF77bVXHHfcceVYGgDADtttIm748OHxgQ98ICIibrrpppg7d+42+3z3u98t/ZSG888/3yuCAIAua7d4dWqT733vezFy5MhYt25dnHDCCfGNb3wjjjvuuFi3bl3cfvvtMW3atIiIGDhwYHzlK1+p8GoBALZvt4q4IUOGxB133BH/83/+z/jLX/4S3/jGN7bZZ+DAgTFz5syt3pYEAKCr2W1upzb56Ec/Gs8++2xccMEFMXDgwNhjjz1ir732iqFDh8bkyZNj4cKFcfjhh1d6mQAAraoqiqKo9CJ2Bw0NDaX3nVu+fHmXeHNgAKA8OqMDdrsrcQAAbwciDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgoepKL2B3sWnTptL2ypUrK7gSAKDcmv/d37wJdoWIK5NXXnmltD18+PAKrgQAqKRXXnklDj744F0+j9upAAAJVRVFUVR6EbuD9evXx3PPPRcREfvtt19UV3fsRdCVK1eWrvDNnz8/3vnOd3bo+dmWmZefmVeGuZefmZdfZ89806ZNpbty7373u6N37967fE63U8ukd+/eMWzYsLJ8rne+851RV1dXls/FFmZefmZeGeZefmZefp018464hdqc26kAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJOTNfgEAEnIlDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIRH3NvCHP/whLrzwwhg8eHDU1NREv379Yvjw4XHNNdfE2rVrK728FF5++eV44IEHYtKkSXHiiSfGvvvuG1VVVVFVVRVnnXXWDp9v1qxZMXbs2Kirq4tevXpFXV1djB07NmbNmtXxi0/s6aefjiuvvDJOPPHEGDBgQPTq1Stqa2tj4MCBcdZZZ8Xjjz++Q+cz99b95S9/idtvvz2+8pWvxKhRo+Lwww+PPffcM3r27Bn7779/HHvssTFlypR47bXX2nU+8951EydOLP2/pqqqKn7+85+3eYy5t0/zubb269hjj23zXF125gWpPfDAA8Wee+5ZRESLv4488shi6dKllV5ml7e9+UVEMW7cuHafp7GxsTjvvPNaPd95551XNDY2dt4Xk8QHP/jBVufU9OuTn/xksWHDhlbPZe7t8/DDD7dr5vvuu28xa9as7Z7HvDvGM888U1RXV281t8cee2y7+5v7jmnPn/WIKEaNGrXdc3T1mYu4xJ555plijz32KCKiqK2tLa644oriqaeeKh555JHi3HPPLf0BGzRoUPHGG29UerldWvP/IAcMGFCccMIJOxVx3/jGN0rHDRkypJgxY0Yxf/78YsaMGcWQIUNKj33zm9/svC8micMOO6yIiKJ///7F+eefX/z4xz8u5s+fX8ydO7e49tpriwMPPLA0rzPPPLPVc5l7+zz88MPFgAEDik996lPF9773veLuu+8u5s6dWzz55JPFHXfcUZx22mlF9+7di4goevbsWfzqV79q8Tzmves2b95cDBs2rIiIYv/9929XxJn7jmmax+c+97niueee2+6vF154Ybvn6OozF3GJHXvssUVEFNXV1cVTTz21zeNTpkwp/QH79re/XYEV5jFp0qTi/vvvL1566aWiKIpi2bJlOxxxzz//fOlf1UOHDi3Wrl271eNr1qwphg4dWvqeLVmypKO/jFROOumk4o477ig2bdrU4uOvvPJKMXDgwNL3Yc6cOS3uZ+7tt71ZN3fPPfeUZj527NhtHjfvjjF16tTSP7K//vWvtxlx5r7jmmZ6ySWX7NTxGWYu4pKaP39+6Q/oZz7zmRb32bx5czF48OAiIoq999672LhxY5lXmdfORNznP//50jFz585tcZ+5c+eW9vniF7/YgSt+e7r//vtL8/rSl77U4j7m3vEGDRpUuq36Vua96/7whz8UtbW1pWi75JJL2ow4c99xuxpxGWbuhQ1J3XvvvaXts88+u8V9unXrFp/61KciImLVqlXtesIsO6coirjvvvsiImLQoEExYsSIFvcbMWJEHHnkkRGx5XtYFEXZ1phR8yccL126dJvHzb1z1NTURETE+vXrt/q4eXeMz3/+87F69eoYN25cu55Ub+7ll2XmIi6pplft1dTUxPve977t7jdq1KjS9hNPPNHp69pdLVu2LFasWBERW8+8JU2PNzQ0xIsvvtjZS0tt48aNpe1u3bb935W5d7xFixbFM888ExFb/vJqzrx33Z133hkPPPBA9OvXL66++up2HWPu5Zdl5iIuqUWLFkVExOGHHx7V1dXb3a/5/4SbjqHjNZ/tW//ieyvfk/abPXt2abuluZp7x1i7dm08//zzce2118Zxxx0XmzdvjoiI888/f6v9zHvX/OlPfyrNdPLkybHffvu16zhz3zV33XVXHHnkkdGnT5/o27dvHHHEETFu3Lh47LHHtntMlplv/29/uqz169fHq6++GhERdXV1re679957R01NTaxZsyaWL19ejuXtlprPtq3vyYABA1o8jq01NjbGVVddVfr96aefvs0+5r7zbr311u0+FSMi4sILL4xPfOITW33MvHfNxIkT46WXXop/+Id/iHPOOafdx5n7rvnNb36z1e+XLFkSS5YsiR/84AfxsY99LG699dbYc889t9ony8xFXEJvvPFGabu2trbN/ZsibvXq1Z25rN3ajnxPmp5vFBG+J62YOnVqzJ8/PyIiTjnllBg6dOg2+5h7x3vve98bN954Y7z//e/f5jHz3nlPPPFETJ8+Paqrq+PGG2+Mqqqqdh9r7jtnjz32iDFjxsQ//uM/xqBBg6K2tjZeeeWVmD17dtx4443x2muvxb333hv19fXx8MMPR48ePUrHZpm5iEuo+ZONe/bs2eb+vXr1ioiIdevWddqadnc78j1p+n5E+J5sz+zZs+NrX/taRETsv//+8f3vf7/F/cx9533sYx8rhfG6deti6dKlceedd8Y999wTn/jEJ+K6666Lk08+eatjzHvnbNy4Mc4777woiiIuuOCCePe7371Dx5v7zlmxYkXstdde23z8+OOPj/Hjx8eJJ54YCxcujNmzZ8f3v//9+NKXvlTaJ8vMPScuod69e5e2mz/xe3s2bNgQERF9+vTptDXt7nbke9L0/YjwPWnJf/3Xf8Upp5wSmzZtil69esWdd94ZBxxwQIv7mvvO22uvveLv//7v4+///u9j2LBhccYZZ8Tdd98dP/jBD+KFF16I+vr6uPXWW7c6xrx3zpVXXhmLFi2Kv/u7v4tLLrlkh483953TUsA1OeCAA+LHP/5xKdCuv/76rR7PMnMRl1Dfvn1L2+25dLtmzZqIaN+tV3bOjnxPmr4fEb4nb7Vs2bI44YQTYtWqVdG9e/eYMWNGq68MM/eO98lPfjJOO+20aGxsjC9+8YuxatWq0mPmveMWL14c3/nOdyJiSyg0v/XWXubeOQ499NA4/vjjI2LL8+T++Mc/lh7LMnO3UxPq3bt37LvvvvHqq69GQ0NDq/uuWrWq9Aes+ZMv6VjNn/ja1vek+RNffU/+5o9//GN8+MMfjj/+8Y9RVVUVN998c5xyyimtHmPunaO+vj7uvPPOWLNmTTz00EPxL//yLxFh3jtj6tSpsXHjxjj00ENj7dq1cfvtt2+zz69//evS9qOPPhovvfRSRER89KMfjZqaGnPvREcddVTMnDkzIrbcfu3fv39E5PmzLuKSGjx4cDz++OOxZMmS2LRp03bfZmTx4sVbHUPnOOqoo0rbzWfeEt+Tbb366qtx/PHHxwsvvBARW65YNL1RdWvMvXM0f+uL3//+96Vt895xTbfaXnjhhTjzzDPb3P+yyy4rbS9btixqamrMvRNt7815s8zc7dSkjjnmmIjYchn3l7/85Xb3a/4+WyNHjuz0de2uDjnkkNK/4JrPvCVz5syJiIgDDzwwDj744M5eWpf35z//OT7ykY+U3gbgqquuii984QvtOtbcO0fTm5xGbH17yLwrw9w7T/O3H2macUSemYu4pD72sY+Vtm+55ZYW92lsbIwf/OAHEbHlCZ7HHXdcOZa2W6qqqor6+vqI2PKvsnnz5rW437x580r/aquvr9+htxl4O1q7dm2cdNJJ8fTTT0dExDe/+c346le/2u7jzb1z3HXXXaXt5q+kNO8dd+utt0ax5eeUb/dX8xc7PPbYY6WPNwWBuXeOF154IR5++OGI2PL8uAMPPLD0WJqZl/MHtdKxPvCBDxQRUVRXVxdPPfXUNo9PmTJll38A8O5q2bJlpdmNGzeuXcf89re/Laqrq4uIKIYOHVqsXbt2q8fXrl1bDB06tPQ9+93vftcJK89jw4YNxQknnFCa8/nnn79T5zH39rvllluKdevWtbrPtddeW/qeHHzwwcWbb7651ePm3fEuueSS0swfe+yxFvcx9x3zk5/8ZJs/u8299NJLxZAhQ0pz/+53v7vNPhlmLuISe/rpp4s+ffoUEVHU1tYWV155ZTF37tzi0UcfLc4777zSH86BAwcWf/nLXyq93C7t8ccfL2655ZbSr6uvvro0v5EjR2712C233LLd83zta18rHTdkyJDi9ttvL37xi18Ut99++1b/w/j6179evi+uixo7dmxpHh/60IeKZ599tnjuuee2++u3v/3tds9l7u1z0EEHFf369SvOPffc4rbbbiueeOKJ4plnnikef/zx4v/8n/9TjBw5sjSrnj17Fg8//HCL5zHvjtWeiCsKc98RBx10UNG/f/9i/PjxxY9+9KPiqaeeKhYuXFg8/PDDxTe/+c1in332Kc3rmGOOKdavX9/iebr6zEVccj/5yU+Kd7zjHaU/SG/9NXDgwOL555+v9DK7vHHjxm13hi392p7NmzcXn/70p1s99pxzzik2b95cxq+ua9qReUdEcdBBB233XObePgcddFC7Zl1XV1f87Gc/2+55zLtjtTfizL392vtn/eMf/3ixatWq7Z6nq89cxL0NvPjii8UFF1xQDBw4sNhjjz2Kvfbaqxg6dGgxefLkYs2aNZVeXgodFXFNZs6cWdTX1xf9+/cvevbsWfTv37+or68vHnzwwTJ8NTl0ZMQ1MffWLVmypLjxxhuLf/7nfy7e8573FAcccEBRXV1d1NbWFocddljx8Y9/vLjlllva/f8N8+4Y7Y24Jubetp///OfFt7/97eKf/umfioEDBxb9+vUrqquri7322qt497vfXXzmM59p8WlI29NVZ15VFNt5fS0AAF2WV6cCACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACT0/wNkgixrW/CnQwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 659, "width": 312 } }, "output_type": "display_data" } ], "source": [ "gdf = geopandas.GeoDataFrame(\n", " geometry=geometries,\n", " data=numpy.ones((n_cols * n_rows + n_cols, 1), int),\n", " columns=[\"var\"],\n", ")\n", "gdf.plot(edgecolor=\"black\");" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.202068Z", "start_time": "2022-10-29T17:23:21.191055Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.668838Z", "iopub.status.busy": "2025-04-06T22:39:29.668757Z", "iopub.status.idle": "2025-04-06T22:39:29.672701Z", "shell.execute_reply": "2025-04-06T22:39:29.672538Z", "shell.execute_reply.started": "2025-04-06T22:39:29.668829Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
vargeometry
01POLYGON ((10 0, 10 10, 0 10, 0 0, 10 0))
11POLYGON ((20 0, 20 10, 10 10, 10 0, 20 0))
21POLYGON ((30 0, 30 10, 20 10, 20 0, 30 0))
31POLYGON ((40 0, 40 10, 30 10, 30 0, 40 0))
41POLYGON ((50 0, 50 10, 40 10, 40 0, 50 0))
\n", "
" ], "text/plain": [ " var geometry\n", "0 1 POLYGON ((10 0, 10 10, 0 10, 0 0, 10 0))\n", "1 1 POLYGON ((20 0, 20 10, 10 10, 10 0, 20 0))\n", "2 1 POLYGON ((30 0, 30 10, 20 10, 20 0, 30 0))\n", "3 1 POLYGON ((40 0, 40 10, 30 10, 30 0, 40 0))\n", "4 1 POLYGON ((50 0, 50 10, 40 10, 40 0, 50 0))" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf.head()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.214821Z", "start_time": "2022-10-29T17:23:21.204646Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.673119Z", "iopub.status.busy": "2025-04-06T22:39:29.673021Z", "iopub.status.idle": "2025-04-06T22:39:29.676584Z", "shell.execute_reply": "2025-04-06T22:39:29.676419Z", "shell.execute_reply.started": "2025-04-06T22:39:29.673113Z" } }, "outputs": [], "source": [ "w = libpysal.weights.Queen.from_dataframe(gdf)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.223053Z", "start_time": "2022-10-29T17:23:21.217867Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.676973Z", "iopub.status.busy": "2025-04-06T22:39:29.676883Z", "iopub.status.idle": "2025-04-06T22:39:29.678636Z", "shell.execute_reply": "2025-04-06T22:39:29.678477Z", "shell.execute_reply.started": "2025-04-06T22:39:29.676967Z" } }, "outputs": [ { "data": { "text/plain": [ "array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int32)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w.component_labels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we have two components, one with 5 areas in the south and a larger northern component with 50 areas." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.229719Z", "start_time": "2022-10-29T17:23:21.225802Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.678977Z", "iopub.status.busy": "2025-04-06T22:39:29.678905Z", "iopub.status.idle": "2025-04-06T22:39:29.680323Z", "shell.execute_reply": "2025-04-06T22:39:29.680175Z", "shell.execute_reply.started": "2025-04-06T22:39:29.678970Z" } }, "outputs": [], "source": [ "numpy.unique(w.component_labels)\n", "rng = numpy.random.default_rng(2021) # set seed for random numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first attempt to solve where the threshold=5. The small southern component is feasible in this case so we should get a solution:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.535670Z", "start_time": "2022-10-29T17:23:21.232947Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.680759Z", "iopub.status.busy": "2025-04-06T22:39:29.680622Z", "iopub.status.idle": "2025-04-06T22:39:29.772828Z", "shell.execute_reply": "2025-04-06T22:39:29.772525Z", "shell.execute_reply.started": "2025-04-06T22:39:29.680753Z" } }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = MaxP(gdf, w, \"var\", \"var\", 5, 2, policy=\"keep\")\n", "model.solve()\n", "gdf[\"region\"] = model.labels_\n", "gdf.explore(column=\"region\", categorical=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we increase the threshold beyond the level of feasibility for this small component:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.549175Z", "start_time": "2022-10-29T17:23:21.538673Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.773539Z", "iopub.status.busy": "2025-04-06T22:39:29.773347Z", "iopub.status.idle": "2025-04-06T22:39:29.777256Z", "shell.execute_reply": "2025-04-06T22:39:29.777050Z", "shell.execute_reply.started": "2025-04-06T22:39:29.773530Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No feasible solution!\n" ] } ], "source": [ "try:\n", " model = MaxP(gdf, w, \"var\", \"var\", 6, 2, policy=\"keep\")\n", " model.solve()\n", "except IndexError as e:\n", " if str(e).startswith(\"list index out of range\"):\n", " print(\"No feasible solution!\")\n", " else:\n", " raise e" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution: `multiple attach`\n", "We can handle the small component by attaching each area in the infeasible component with its nearest neighbor area belonging to a feasible component. This is done with setting the policy to `multiple`." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.644765Z", "start_time": "2022-10-29T17:23:21.555283Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.778906Z", "iopub.status.busy": "2025-04-06T22:39:29.778835Z", "iopub.status.idle": "2025-04-06T22:39:29.805935Z", "shell.execute_reply": "2025-04-06T22:39:29.805649Z", "shell.execute_reply.started": "2025-04-06T22:39:29.778899Z" }, "tags": [] }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = MaxP(gdf, w, \"var\", \"var\", 6, 2, policy=\"multiple\")\n", "model.solve()\n", "gdf[\"region\"] = model.labels_\n", "gdf.explore(column=\"region\", categorical=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `multiple` may allow areas from the same original infeasible component to be assigned to different regions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution: `single attach`\n", "We can handle the small component by attaching an infeasible component to a feasible component using a single join. That join is based on the minimum nearest neighbor distance separating an area in the infeasible component with an area from a feasible component. This will always result in the areas from the original infeasible component being assigned to the same region." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.714892Z", "start_time": "2022-10-29T17:23:21.647742Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.806470Z", "iopub.status.busy": "2025-04-06T22:39:29.806380Z", "iopub.status.idle": "2025-04-06T22:39:29.831489Z", "shell.execute_reply": "2025-04-06T22:39:29.831239Z", "shell.execute_reply.started": "2025-04-06T22:39:29.806462Z" }, "tags": [] }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = MaxP(gdf, w, \"var\", \"var\", 6, 2, policy=\"single\")\n", "model.solve()\n", "gdf[\"region\"] = model.labels_\n", "gdf.explore(column=\"region\", categorical=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution: `drop`\n", "\n", "An alternative approach is to drop the areas belonging to infeasible components, and then solve for the remaining feasible components." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:21.881777Z", "start_time": "2022-10-29T17:23:21.717282Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.832204Z", "iopub.status.busy": "2025-04-06T22:39:29.832116Z", "iopub.status.idle": "2025-04-06T22:39:29.884579Z", "shell.execute_reply": "2025-04-06T22:39:29.884312Z", "shell.execute_reply.started": "2025-04-06T22:39:29.832196Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAUnCAYAAAAy/TFcAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAewgAAHsIBbtB1PgAAU45JREFUeJzt3XucV3W96P/3wHBzhlS8/cJhe0ewXSceAeHGQttpWzQmSd26O4Xm1q4kltJVtOMlQRM7bk9uHnjrPAovbS8pSrnVwAscIjFtByWIxRBub1RyF2b9/qD5NsgwM8DM9ztveT4fDx6PxXzXWvOZ96C8WOv7/U5VURRFAACQSrdKLwAAgB0n4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACRUXekF7C7Wr18fzz33XERE7LffflFdbfQAsLvYtGlTvPLKKxER8e53vzt69+69y+dUEmXy3HPPxfDhwyu9DACgwubPnx/Dhg3b5fO4nQoAkJArcWWy3377lbbnz58f73znOyu4GgCgnFauXFm6I9e8CXaFiCuT5s+Be+c73xl1dXUVXA0AUCkd9bx4t1MBABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICEOjXiXn755XjggQdi0qRJceKJJ8a+++4bVVVVUVVVFWeddVa7zrF+/fq47777Yvz48fH+978/+vXrFz169Ih+/frF0UcfHZdeemmsXLmy3Wtau3ZtXH311TF8+PDo169f1NbWxuDBg+PCCy+MP/zhDzv5lQIAlFd1Z578gAMO2KXjn3322TjmmGPijTfe2OaxVatWxbx582LevHlx7bXXxvTp0+P0009v9XxLly6Nk046KX77299u9fHFixfH4sWLY/r06fGjH/0oRo8evUvrBgDobGW7nTpgwIA44YQTduiYv/zlL6WAGzlyZHznO9+Jhx9+OJ5++un46U9/Gp/5zGeie/fu8cYbb8S//Mu/xEMPPbTdc61evTpOPvnkUsCde+658cgjj8RTTz0VV1xxRdTW1saf//znOO200+LZZ5/d+S8UAKAMOvVK3KRJk2LYsGExbNiwOOCAA+LFF1+MQw45pN3Hd+vWLU4//fS45JJL4qijjtrm8RNOOCFOPPHEOOWUU2Lz5s0xfvz4eP7556Oqqmqbfa+55ppYvHhxRERMmTIlLrrootJjRx99dBx33HHxwQ9+MNauXRsTJkyIRx99dCe+YgCA8qgqiqIo1ydrHnHjxo2LW2+9tUPOe+qpp8Z//Md/RETE008/HUOGDNnq8TfffDP233//+NOf/hSDBw+OX//619Gt27YXIT/72c/Gv//7v0dExIIFC+J973tfh6wvIqKhoSEGDBgQERHLly+Purq6Djs3ANC1dUYHvC1enXrccceVtpcuXbrN4z//+c/jT3/6U0RsiceWAi4itnqxxd13392hawQA6Ehvi4jbsGFDabulQHv88cdL26NGjdrueYYOHRo1NTUREfHEE0904AoBADrW2yLiZs+eXdoeNGjQNo8vWrSo1cebVFdXx2GHHbbNMQAAXU2nvrChHH71q1/FzJkzIyLiXe96V4svgFi+fHlERNTU1MRee+3V6vkGDBgQzz77bLzyyiuxYcOG6NWrV7vW0dDQ0OrjO/JedgAAbUkdcRs2bIh//dd/jc2bN0dExJVXXtnifk1vU1JbW9vmOZtup0ZseVuS9kZc05MVK+muu+6KSZMmtfi+enSOdevWxerVq6O2tjb69OlT6eXsFsy8Msy9/My8Mvr27RuXXXZZnHrqqZVeSptSR9wXv/jFWLBgQURsecHCmDFjWtxv/fr1ERHRs2fPNs/ZPNrWrVvXAassn0mTJpXeRoXyev311yu9hN2OmVeGuZefmZffxRdfLOI603e+852YPn16RES8733vixtuuGG7+/bu3TsiIjZu3NjmeZu/SGJH/uXTdMt2e1auXBnDhw9v9/l2RukKXFW36F6zd6d+LrbYvPq1LRtmXjZmXhnmXn5mXn6b16yKKBrT3NFKGXH//u//Ht/4xjciIuLII4+Mhx56aKvboG/Vt2/fiNhye7Qta9asKW235/Zrk670vm/da/aOui/cVull7BZ+P2VMRNFo5mVk5pVh7uVn5uXXcMO4v8VzAulenTpjxoz4/Oc/HxERBx10UPznf/5n7Lfffq0e0xRYa9asKb1f3PY0XVHbb7/92v18OACAcksVcT/5yU/iU5/6VDQ2NsY73/nOeOSRR9p1Baz5K1Zbe87Ypk2bSm8WPHjw4F1fMABAJ0kTcY888kicfvrpsWnTpthnn33i4YcfLr2nW1uOOeaY0nbz95R7qwULFpRup44cOXLXFgwA0IlSRNxTTz0V9fX1sWHDhnjHO94RP/3pT+Nd73pXu48/9thjY88994yIiNtuuy229+Nim/8s11NOOWWX1gwA0Jm6fMQ988wzcdJJJ8WaNWuipqYmHnzwwR3+wfQ9e/aML33pSxGx5ScxXHPNNdvsM3fu3LjpppsiYsuP5ho2bNiuLx4AoJN06qtTn3jiiViyZEnp96+++mppe8mSJVtd+YrY+gfQR2z5YfYf+chHSi9GuPzyy2PPPfeMX//619v9nPvvv3/sv//+23z8oosuijvuuCN+97vfxcSJE2PJkiVxxhlnRJ8+feKxxx6LK6+8MjZt2hR9+vSJ6667boe/VgCAcurUiJs+fXrcdlvLL4t+8skn48knn9zqY2+NuMcffzxefvnl0u8vuOCCNj/nJZdcEpdeeuk2H+/bt2/MnDkzRo8eHc8//3xMmzYtpk2bttU+73jHO+KHP/xhvPe9723z8wAAVFKXv53akQ4//PBYuHBhTJ48OYYOHRp77bVX7LHHHnHkkUfGBRdcEM8++2ycfPLJlV4mAECbOvVK3K233rrNLdMdcdZZZ21zdW5X1dTUxMSJE2PixIkdel4AgHLara7EAQC8XYg4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABKqKoqiqPQidgcNDQ0xYMCAiIhYvnx51NXVdfjn2GeffeL111+PqOoW3Wv27vDzs63Nq1/bsmHmZWPmlWHu5Wfm5dc08379+sVrr73WoefujA6o3uUz0GWsXr16y0bR+Lf/+CkPMy8/M68Mcy8/My+76uoceZRjlbRLbW3tlitxERFV7pSXRdH4t20zLw8zrwxzLz8zL7+/zrxHjx4VXkj7iLi3kT59+kRERPfafaLuC7dVeDW7h99PGRNRNJp5GZl5ZZh7+Zl5+TXcMC7VVU9pDwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQUKdG3MsvvxwPPPBATJo0KU488cTYd999o6qqKqqqquKss87a4fPNmjUrxo4dG3V1ddGrV6+oq6uLsWPHxqxZs9p9jrVr18bVV18dw4cPj379+kVtbW0MHjw4LrzwwvjDH/6ww2sCAKiE6s48+QEHHNAh5ymKIj772c/GtGnTtvr4ihUr4p577ol77rknzjvvvLjxxhujqqpqu+dZunRpnHTSSfHb3/52q48vXrw4Fi9eHNOnT48f/ehHMXr06A5ZNwBAZynb7dQBAwbECSecsFPHfutb3yoF3JAhQ2LGjBkxf/78mDFjRgwZMiQiIqZNmxYXX3zxds+xevXqOPnkk0sBd+6558YjjzwSTz31VFxxxRVRW1sbf/7zn+O0006LZ599dqfWCQBQLp16JW7SpEkxbNiwGDZsWBxwwAHx4osvxiGHHLJD51iyZElMmTIlIiKGDh0ac+bMiT59+kRExLBhw2LMmDExatSoWLBgQUyePDnOPvvsOOyww7Y5zzXXXBOLFy+OiIgpU6bERRddVHrs6KOPjuOOOy4++MEPxtq1a2PChAnx6KOP7uyXDQDQ6Tr1Sty3v/3tOPnkk3fpturUqVNj06ZNERFx/fXXlwKuyR577BHXX399RERs2rQprrvuum3O8eabb8b3vve9iIgYPHhwfOUrX9lmn6OPPjrOOeeciIh47LHH4pe//OVOrxkAoLN16VenFkUR9913X0REDBo0KEaMGNHifiNGjIgjjzwyIiLuvffeKIpiq8d//vOfx5/+9KeIiBg3blx069byl938xRZ33333Lq4eAKDzdOmIW7ZsWaxYsSIiIkaNGtXqvk2PNzQ0xIsvvrjVY48//vg2+7Vk6NChUVNTExERTzzxxM4sGQCgLDr1OXG7atGiRaXtQYMGtbpv88cXLVq01XPv2nue6urqOOyww+LZZ5/d6pj2aGhoaPXxlStX7tD5AABa06Ujbvny5aXturq6VvcdMGBAi8c1/31NTU3stddebZ7n2WefjVdeeSU2bNgQvXr1atdam3/+Slm3bl1ERGxesyoabhhX4dXsJorGiDDzsjLzyjD38jPzstu8+rWI+Nvfp11dl464N954o7RdW1vb6r5Nt0EjtrydSEvnaescLZ2nvRHXFZS+7qKx9AeRMjHz8jPzyjD38jPzsquu7tJ5VNKlV7l+/frSds+ePVvdt3lsvbWgm87T1jnaOk9r3nr1761WrlwZw4cPb/f5dkZtbW28/vrrW35T1aWf7vj28dd/KUeEmZeLmVeGuZefmZffX2feo0ePCi+kfbp0xPXu3bu0vXHjxlb33bBhQ2n7rW9D0nSets7R1nla09bt3nJoWm/32n2i7gu3VXg1u4ffTxkTUTSaeRmZeWWYe/mZefk13DAu1VXPLp32ffv2LW2/9RbpW61Zs6a0/dbbpk3naescbZ0HAKCr6NIR1/zqVluv/mx+O/OtLzJoOs+aNWtK7xfX1nn222+/VM+HAwB2L1064o466qjSdtOPzNqe5o8PHjx4p86zadOmWLp0aYvnAADoSrp0xB1yyCHRv3//iIiYPXt2q/vOmTMnIiIOPPDAOPjgg7d67Jhjjiltt3aeBQsWlG6njhw5cmeWDABQFl064qqqqqK+vj4itlxBmzdvXov7zZs3r3SFrb6+PqqqqrZ6/Nhjj40999wzIiJuu+22bX4sV5Nbb721tH3KKafs6vIBADpNl464iIgJEyaU3q9l/Pjx27ztx7p162L8+PERseV9XSZMmLDNOXr27Blf+tKXImLLT2+45pprttln7ty5cdNNN0XElh/NNWzYsI78MgAAOlSnvsXIE088EUuWLCn9/tVXXy1tL1myZKsrXxFb/wD6JgMHDowLL7wwrrrqqliwYEGMHDkyvvrVr8Zhhx0WS5cujcmTJ8fChQsjIuKiiy6KI444osW1XHTRRXHHHXfE7373u5g4cWIsWbIkzjjjjOjTp0889thjceWVV8amTZuiT58+cd111+3y1w4A0Jk6NeKmT58et93W8nvbPPnkk/Hkk09u9bGWIi4i4oorroiXX345br755li4cGGcccYZ2+xzzjnnxOWXX77dtfTt2zdmzpwZo0ePjueffz6mTZsW06ZN22qfd7zjHfHDH/4w3vve97b+hQEAVFiXv50aEdGtW7e46aabYubMmVFfXx/9+/ePnj17Rv/+/aO+vj4efPDBmD59enTr1vqXc/jhh8fChQtj8uTJMXTo0Nhrr71ijz32iCOPPDIuuOCCePbZZ+Pkk08u01cFALDzOvVK3K233rrNLdNdMXr06Bg9evQunaOmpiYmTpwYEydO7KBVAQCUX4orcQAAbE3EAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJBQVVEURaUXsTtoaGiIAQMGRETE8uXLo66ursM/xz777BOvv/56RFW36F6zd4efn21tXv3alg0zLxszrwxzLz8zL7+mmffr1y9ee+21Dj13Z3RA9S6fgS5j9erVWzaKxr/9x095mHn5mXllmHv5mXnZVVfnyKMcq6Rdamtrt1yJi4iocqe8LIrGv22beXmYeWWYe/mZefn9deY9evSo8ELaR8S9jfTp0yciIrrX7hN1X7itwqvZPfx+ypiIotHMy8jMK8Pcy8/My6/hhnGprnpKewCAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmlibiNGzfGTTfdFP/0T/8U73znO6NXr15RW1sbRx55ZHz605+OefPmtes8s2bNirFjx0ZdXV306tUr6urqYuzYsTFr1qxO/goAADpOdaUX0B7Lly+Pk046KZ577rmtPr5x48b43e9+F7/73e/illtuiQsuuCC++93vRlVV1TbnKIoiPvvZz8a0adO2+viKFSvinnvuiXvuuSfOO++8uPHGG1s8HgCgK+nyV+I2bdq0VcC95z3viVtvvTXmzp0bP/vZz2LSpElRU1MTERFTp06Na665psXzfOtb3yoF3JAhQ2LGjBkxf/78mDFjRgwZMiQiIqZNmxYXX3xxGb4qAIBd0+WvxN13332lgDv66KPj8ccfj+7du5ceP/7442PMmDFx9NFHx5tvvhnf+c534oILLojq6r99aUuWLIkpU6ZERMTQoUNjzpw50adPn4iIGDZsWIwZMyZGjRoVCxYsiMmTJ8fZZ58dhx12WBm/SgCAHdPlr8Q9+eSTpe2vf/3rWwVck/e9731x8sknR0TEqlWrYvHixVs9PnXq1Ni0aVNERFx//fWlgGuyxx57xPXXXx8RW678XXfddR35JQAAdLguH3EbN24sbR966KHb3a/5lbMNGzaUtouiiPvuuy8iIgYNGhQjRoxo8fgRI0bEkUceGRER9957bxRFsUvrBgDoTF0+4gYOHFjafuGFF7a739KlSyMioqqqKo444ojSx5ctWxYrVqyIiIhRo0a1+rmaHm9oaIgXX3xxZ5cMANDpunzEnXnmmfGOd7wjIiImT54cmzdv3mafhQsXxsyZMyMi4owzzijtHxGxaNGi0vagQYNa/VzNH29+HABAV9PlX9iw3377xa233hqf+MQn4sknn4xhw4bFhAkTYuDAgbF69ep48skn47vf/W5s3Lgx3vve98a111671fHLly8vbdfV1bX6uQYMGNDice3R0NDQ6uMrV67cofMBALSmy0dcRMQpp5wSCxYsiGuvvTZuvvnmGDdu3FaPH3DAAfHtb387zjvvvNLbjTR54403Stu1tbWtfp7mx65evXqH1tg8ACtl3bp1ERGxec2qaLhhXBt70yGKxogw87Iy88ow9/Iz87LbvPq1iPjb36ddXYqIe/PNN+NHP/pR3H///S2+4OC///u/Y8aMGTFw4MA46aSTtnps/fr1pe2ePXu2+nl69epV2s7yDWyuFJ5FY+kPImVi5uVn5pVh7uVn5mXX/G3KurIuv8o1a9bE6NGjY86cOdG9e/eYOHFinH322XHooYfG+vXr4//9v/8X/+t//a944okn4qMf/WhMnTo1zj///NLxvXv3Lm03f6VrS5q/qvWtb0PSlrZuv65cuTKGDx++Q+fcUbW1tfH6669v+U1Vl3+649vDX/+lHBFmXi5mXhnmXn5mXn5/nXmPHj0qvJD26fIRd8kll8ScOXMiIuKmm27a6lZqz5494/jjj4/jjjsuTjjhhHjsscfiy1/+chx33HHxnve8JyIi+vbtW9q/rVuka9asKW23dev1rdp6vl05NIVn99p9ou4Lt1V4NbuH308ZE1E0mnkZmXllmHv5mXn5NdwwLtVVzy6d9kVRxC233BIRW95q5K3PhWtSXV0dl112WURENDY2lo6J2Dqu2nrxQfOraV3hOW4AANvTpSPuv//7v0u3B5t+vun2vO997yttN/+JDUcddVSLH29J88cHDx68Q2sFACinLh1xzZ9Y2PRjs7bnzTffbPG4Qw45JPr37x8REbNnz271HE23bQ888MA4+OCDd3S5AABl06Ujrl+/fqU37p07d26rIdc80A455JDSdlVVVdTX10fElitt8+bNa/H4efPmla7E1dfXR1VV1S6vHwCgs3TpiOvWrVvpLUP++Mc/xhVXXNHifqtWrYqvfvWrpd+ffPLJWz0+YcKE0tW58ePHb/P2IevWrYvx48dHxJareBMmTOioLwEAoFN06YiLiJg0aVLsscceERFx6aWXxpgxY+I//uM/YuHChTF37tyYOnVqvPe9743f/OY3ERHxj//4j3HCCSdsdY6BAwfGhRdeGBERCxYsiJEjR8Ydd9wRCxYsiDvuuCNGjhwZCxYsiIiIiy66aKufvQoA0BV1+bcYGTRoUNx3331x5plnxquvvhr3339/3H///S3u+6EPfSjuuuuuFh+74oor4uWXX46bb745Fi5cGGecccY2+5xzzjlx+eWXd+j6AQA6Q5e/EhcR8eEPfzgWL14ckydPjmOPPTb222+/6NGjR/Tp0ycOOeSQOP300+Pee++N//zP/4y99967xXN069Ytbrrpppg5c2bU19dH//79o2fPntG/f/+or6+PBx98MKZPnx7duqUYCQCwm+vyV+Ka7LPPPjFx4sSYOHHiLp1n9OjRMXr06A5aFQBAZbjsBACQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AIKGqoiiKSi9id9DQ0BADBgyIiIjly5dHXV1dh3+OffbZJ15//fWIqm7RvWbvDj8/29q8+rUtG2ZeNmZeGeZefmZefk0z79evX7z22msdeu7O6IDqXT4DXcbq1au3bBSNf/uPn/Iw8/Iz88ow9/Iz87Krrs6RRzlWSbvU1tZuuRIXEVHlTnlZFI1/2zbz8jDzyjD38jPz8vvrzHv06FHhhbSPiHsb6dOnT0REdK/dJ+q+cFuFV7N7+P2UMRFFo5mXkZlXhrmXn5mXX8MN41Jd9ZT2AAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJpYu4V199NaZMmRIjR46M/+//+/+iV69e0b9//3j/+98fF110UcydO7fNc8yaNSvGjh0bdXV10atXr6irq4uxY8fGrFmzyvAVAADsuupKL2BH3HXXXfG5z30uXnvtta0+vnLlyli5cmXMnz8/nn/++bj33ntbPL4oivjsZz8b06ZN2+rjK1asiHvuuSfuueeeOO+88+LGG2+MqqqqzvoyAAB2WZqI+8EPfhBnn312NDY2xv777x+f+9zn4phjjol+/frFSy+9FEuXLo37778/evTosd1zfOtb3yoF3JAhQ2LixIlx2GGHxdKlS2PKlCmxcOHCmDZtWuy3335x+eWXl+tLAwDYYSkibtGiRXHeeedFY2NjfOADH4j7778/9txzz232Gz9+fGzcuLHFcyxZsiSmTJkSERFDhw6NOXPmRJ8+fSIiYtiwYTFmzJgYNWpULFiwICZPnhxnn312HHbYYZ33RQEA7IIUz4kbP358bNiwIfbdd9+4++67Wwy4Jj179mzx41OnTo1NmzZFRMT1119fCrgme+yxR1x//fUREbFp06a47rrrOmbxAACdoMtH3OLFi+ORRx6JiIgvfvGLse++++7wOYqiiPvuuy8iIgYNGhQjRoxocb8RI0bEkUceGRER9957bxRFsZOrBgDoXF0+4u66667S9mmnnVbaXrVqVTz//PPbvMihJcuWLYsVK1ZERMSoUaNa3bfp8YaGhnjxxRd3YsUAAJ2vy0fcvHnzIiJizz33jMGDB8cPf/jD+B//439Ev379YuDAgbHvvvvGoYceGt/+9rdj9erVLZ5j0aJFpe1Bgwa1+vmaP978OACArqTLv7DhN7/5TUREHHzwwTF+/Pi44YYbttln2bJlcemll8aPf/zj+OlPfxr9+/ff6vHly5eXtuvq6lr9fAMGDGjxuLY0NDS0+vjKlSvbfS4AgLZ0+Yh7/fXXI2LLc+N+9atfxV577RVXXXVVjB07Nt7xjnfEc889F5MmTYqHHnoofv3rX8dpp50Wjz/+eHTr9reLjG+88UZpu7a2ttXPV1NTU9re3pW9ljSPv0pZt25dRERsXrMqGm4YV+HV7CaKxogw87Iy88ow9/Iz87LbvHrLU7Sa/j7t6rp8xK1ZsyYiIjZs2BDdu3ePhx56aKsXJgwdOjQeeOCBOPnkk+Ohhx6Kp556Ku6+++449dRTS/usX7++tL29V6826dWrV2k7yzexSSk6i8bSH0TKxMzLz8wrw9zLz8zLrrq6y+dRRCSIuN69e5dC7rTTTmvxlaXdunWLq6++Oh566KGIiJgxY8ZWEde7d+/S9vbeR67Jhg0bSttvfRuS1rR163XlypUxfPjwdp9vZ9TW1pauXEZVl3+649vDX/+lHBFmXi5mXhnmXn5mXn5/nXlrPzigK+nyEde3b99SxJ144onb3e9d73pXHHjggbFixYr4xS9+sc05mrR1i7Tpc0W0feu1ubaea1cOTdHZvXafqPvCbRVeze7h91PGRBSNZl5GZl4Z5l5+Zl5+DTeMS3XVs8unffPnmrX3RQkvv/zyVh9vflxbL0BofkWtKzzPDQCgJV0+4t71rneVtjdv3tzqvk2Pv/Ve9lFHHVXaXrx4cavnaP744MGD271OAIBy6vIR98EPfrC0vXTp0lb3feGFFyIi4sADD9zq44ccckjpbUdmz57d6jnmzJlTOsfBBx+8o8sFACiLLh9xY8aMKT3B8O67797ufrNnzy799IYPfOADWz1WVVUV9fX1EbHlSlvTGwi/1bx580pX4urr66OqqmqX1w8A0Bm6fMTts88+8a//+q8REfHwww/H7bffvs0+b7zxRkyYMKH0+8985jPb7DNhwoTSbdbx48dv8/Yh69ati/Hjx0fEltuxzc8HANDVdPmIi4j49re/HX/3d38XERGf/OQnY/z48fHYY4/FL3/5y7j11ltj+PDh8cwzz0RExOc+97kYNmzYNucYOHBgXHjhhRERsWDBghg5cmTccccdsWDBgrjjjjti5MiRsWDBgoiIuOiii+KII44ozxcHALATuvxbjERE7LfffjFr1qwYM2ZMLFmyJP7t3/4t/u3f/m2b/T796U/H9773ve2e54orroiXX345br755li4cGGcccYZ2+xzzjnnxOWXX96h6wcA6GgprsRFbHml6DPPPBNXX311vP/9749+/fpFz549o66uLv75n/85Hn300bjppptafYO+bt26xU033RQzZ86M+vr66N+/f/Ts2TP69+8f9fX18eCDD8b06dO3+pFdAABdUYorcU1qamriwgsvLN0W3VmjR4+O0aNHd9CqAADKzyUnAICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACVUVRVFUehG7g4aGhhgwYEBERCxfvjzq6uo6/HPss88+8frrr0dUdYvuNXt3+PnZ1ubVr23ZMPOyMfPKMPfyM/Pya5p5v3794rXXXuvQc3dGB1Tv8hnoMlavXr1lo2j823/8lIeZl5+ZV4a5l5+Zl111dY48yrFK2qW2tnbLlbiIiCp3ysuiaPzbtpmXh5lXhrmXn5mX319n3qNHjwovpH1E3NtInz59IiKie+0+UfeF2yq8mt3D76eMiSgazbyMzLwyzL38zLz8Gm4Yl+qqp7QHAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEgodcRNnDgxqqqqSr9+/vOft3nMrFmzYuzYsVFXVxe9evWKurq6GDt2bMyaNavzFwwA0EHSRtyvfvWrmDp1arv3L4oiPvOZz8SJJ54Y99xzT6xYsSI2btwYK1asiHvuuSdOPPHE+MxnPhNFUXTiqgEAOkbKiGtsbIxzzz03Nm3aFPvvv3+7jvnWt74V06ZNi4iIIUOGxIwZM2L+/PkxY8aMGDJkSERETJs2LS6++OJOWzcAQEdJGXH/+3//7/jFL34RgwYNinPOOafN/ZcsWRJTpkyJiIihQ4fGk08+GWeccUYMGzYszjjjjHjiiSdi6NChERExefLkWLp0aaeuHwBgV6WLuOXLl5euln3/+9+Pnj17tnnM1KlTY9OmTRERcf3110efPn22enyPPfaI66+/PiIiNm3aFNddd13HLhoAoIOli7jPf/7zsXr16hg3blwce+yxbe5fFEXcd999ERExaNCgGDFiRIv7jRgxIo488siIiLj33ns9Nw4A6NJSRdydd94ZDzzwQPTr1y+uvvrqdh2zbNmyWLFiRUREjBo1qtV9mx5vaGiIF198cZfWCgDQmaorvYD2+tOf/hTnn39+RGx53tp+++3XruMWLVpU2h40aFCr+zZ/fNGiRXHIIYe0e30NDQ2tPr5y5cp2nwsAoC1pIm7ixInx0ksvxT/8wz+068UMTZYvX17arqura3XfAQMGtHhcezQ/tlLWrVsXERGb16yKhhvGVXg1u4miMSLMvKzMvDLMvfzMvOw2r34tIv7292lXlyLinnjiiZg+fXpUV1fHjTfeGFVVVe0+9o033iht19bWtrpvTU1NaXv16tU7vtAKK625aCz9QaRMzLz8zLwyzL38zLzsqqtT5FHXj7iNGzfGeeedF0VRxAUXXBDvfve7d+j49evXl7bbeiVrr169Sts7WuFtXblbuXJlDB8+fIfOuaNqa2vj9ddf3/KbqlRPd8zrr/9SjggzLxczrwxzLz8zL7+/zrxHjx4VXkj7dPmIu/LKK2PRokXxd3/3d3HJJZfs8PG9e/cubW/cuLHVfTds2FDafuvbkLSlrVu15dC05u61+0TdF26r8Gp2D7+fMiaiaDTzMjLzyjD38jPz8mu4YVyqq55dOu0XL14c3/nOdyJiy/u7Nb/d2V59+/Ytbbd1i3TNmjWl7bZuvQIAVFKXvhI3derU2LhxYxx66KGxdu3auP3227fZ59e//nVp+9FHH42XXnopIiI++tGPRk1NzVZXyNp6BWnzW6Jd4YUKAADb06Ujrun25gsvvBBnnnlmm/tfdtllpe1ly5ZFTU1NHHXUUaWPLV68uNXjmz8+ePDgHV0uAEDZdOnbqR3hkEMOif79+0dExOzZs1vdd86cORERceCBB8bBBx/c2UsDANhpXTribr311iiKotVfzV/s8Nhjj5U+3hRhVVVVUV9fHxFbrrTNmzevxc81b9680pW4+vr6HXobEwCAcuvSEddRJkyYUHrPl/Hjx2/z9iHr1q2L8ePHR8SW94aZMGFCuZcIALBDdouIGzhwYFx44YUREbFgwYIYOXJk3HHHHbFgwYK44447YuTIkbFgwYKIiLjoooviiCOOqORyAQDa1KVf2NCRrrjiinj55Zfj5ptvjoULF8YZZ5yxzT7nnHNOXH755RVYHQDAjtktrsRFRHTr1i1uuummmDlzZtTX10f//v2jZ8+e0b9//6ivr48HH3wwpk+fHt267TYjAQASS38l7tJLL41LL7203fuPHj06Ro8e3XkLAgAoA5edAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJFRVFEVR6UXsDhoaGmLAgAEREbF8+fKoq6vr8M+xzz77xOuvvx5R1S261+zd4ednW5tXv7Zlw8zLxswrw9zLz8zLr2nm/fr1i9dee61Dz90ZHVC9y2egy1i9evWWjaLxb//xUx5mXn5mXhnmXn5mXnbV1TnyKMcqaZfa2totV+IiIqrcKS+LovFv22ZeHmZeGeZefmZefn+deY8ePSq8kPYRcW8jffr0iYiI7rX7RN0XbqvwanYPv58yJqJoNPMyMvPKMPfyM/Pya7hhXKqrntIeACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiki7umnn44rr7wyTjzxxBgwYED06tUramtrY+DAgXHWWWfF448/vkPnmzVrVowdOzbq6uqiV69eUVdXF2PHjo1Zs2Z10lcAANCxqiu9gLaMGjUq5syZs83HN27cGM8//3w8//zzcdttt8UnP/nJmD59evTs2XO75yqKIj772c/GtGnTtvr4ihUr4p577ol77rknzjvvvLjxxhujqqqqw78WAICO0uWvxK1YsSIiIvr37x/nn39+/PjHP4758+fH3Llz49prr40DDzwwIiL+7//9v3HWWWe1eq5vfetbpYAbMmRIzJgxI+bPnx8zZsyIIUOGRETEtGnT4uKLL+68LwgAoAN0+StxgwYNiiuvvDI+/vGPR/fu3bd6bMSIEfHJT34yRo4cGb/73e9ixowZ8bnPfS4+8IEPbHOeJUuWxJQpUyIiYujQoTFnzpzo06dPREQMGzYsxowZE6NGjYoFCxbE5MmT4+yzz47DDjus879AAICd0OWvxD3wwANx+umnbxNwTfbdd9/47ne/W/r9j3/84xb3mzp1amzatCkiIq6//vpSwDXZY4894vrrr4+IiE2bNsV1113XAasHAOgcXT7i2uPYY48tbS9dunSbx4uiiPvuuy8itlzZGzFiRIvnGTFiRBx55JEREXHvvfdGURQdv1gAgA7wtoi4jRs3lra7ddv2S1q2bFnpuXWjRo1q9VxNjzc0NMSLL77YcYsEAOhAb4uImz17dml70KBB2zy+aNGiVh9vrvnjzY8DAOhKuvwLG9rS2NgYV111Ven3p59++jb7LF++vLRdV1fX6vkGDBjQ4nFtaWhoaPXxlStXtvtcAABtSR9xU6dOjfnz50dExCmnnBJDhw7dZp833nijtF1bW9vq+Wpqakrbq1evbvc6msdfpW1esyoabhhX6WXsHorGiDDzsjLzyjD38jPzstu8+rVKL2GHpI642bNnx9e+9rWIiNh///3j+9//fov7rV+/vrTd2psBR0T06tWrtL1u3boOWGX59O3bd8tG0ZjuD2J6Zl5+Zl4Z5l5+Zl52X/7ylyu9hHZJG3H/9V//Faecckps2rQpevXqFXfeeWcccMABLe7bu3fv0nbzF0G0ZMOGDaXtt74NSWvauvW6cuXKGD58eLvPtzMuu+yyuPjii7e68kjnWrduXVRXV0ePHj0qvZTdhplXhrmXn5lXxpe//GUR15mWLVsWJ5xwQqxatSq6d+8eM2bMaPVVp6UrVNH2LdI1a9aUttu69dpcW8+1K4dTTz01Tj311EovAwAog3SvTv3jH/8YH/7wh+OPf/xjVFVVxc033xynnHJKq8c0D6y2XoDQ/IpaV3qeGwBAc6ki7tVXX43jjz8+XnjhhYjY8pMXPvWpT7V53FFHHVXaXrx4cav7Nn988ODBO7lSAIDOlSbi/vznP8dHPvKR+M1vfhMREVdddVV84QtfaNexhxxySPTv3z8itn5PuZbMmTMnIiIOPPDAOPjgg3d+wQAAnShFxK1duzZOOumkePrppyMi4pvf/GZ89atfbffxVVVVUV9fHxFbrrTNmzevxf3mzZtXuhJXX18fVVVVu7hyAIDO0eUjbuPGjXHKKafEk08+GRER559/flx++eU7fJ4JEyZEdfWW13GMHz9+m7cPWbduXYwfPz4iIqqrq2PChAm7tnAAgE7U5V+deuaZZ8bPfvaziIj40Ic+FOecc078+te/3u7+PXv2jIEDB27z8YEDB8aFF14YV111VSxYsCBGjhwZX/3qV+Owww6LpUuXxuTJk2PhwoUREXHRRRfFEUcc0TlfEABAB6gqiqKo9CJas6O3NA866KDt/uD6xsbGOPfcc+Pmm2/e7vHnnHNOTJs2Lbp169iLlA0NDaVXuy5fvrxLvCUJAFAendEBXf52akfq1q1b3HTTTTFz5syor6+P/v37R8+ePaN///5RX18fDz74YEyfPr3DAw4AoKN1+dupnXGhcPTo0TF69OgOPy8AQLm45AQAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACCh6kovgI5z1113xaRJk+KNN96o9FJ2G+vWrYvVq1dHbW1t9OnTp9LL2S2YeWWYe/mZeWX07ds3Lrvssjj11FMrvZQ2VRVFUVR6EbuDhoaGGDBgQERELF++POrq6jr8cwwePDgWL17c4ecFgN3JoEGDYtGiRR16zs7oAFfi3kZKV+CqukX3mr0ru5jdxObVr23ZMPOyMfPKMPfyM/Py27xmVUTRmOaOloh7G+pes3fUfeG2Si9jt/D7KWMiikYzLyMzrwxzLz8zL7+GG8b9LZ4T8MIGAICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICEdsuI+8Mf/hAXXnhhDB48OGpqaqJfv34xfPjwuOaaa2Lt2rWVXh4AQJuqK72Acps5c2Z84hOfiD//+c+lj61duzZ+8YtfxC9+8YuYPn16PPjgg3HooYdWcJUAAK3bra7E/epXv4rTTz89/vznP0dtbW1cccUV8dRTT8UjjzwS5557bkRE/Pa3v42TTjopVq9eXeHVAgBs3251JW7ChAmxdu3aqK6ujp/97Gdx9NFHlx770Ic+FEcccURMnDgxFi9eHNdee21MmjSpgqsFANi+3eZK3C9+8Yv4+c9/HhER55xzzlYB1+QrX/lKDB48OCIirrvuunjzzTfLuUQAgHbbbSLu3nvvLW2fffbZLe7TrVu3+NSnPhUREatWrSpFHwBAV7PbRNzjjz8eERE1NTXxvve9b7v7jRo1qrT9xBNPdPq6AAB2xm4TcYsWLYqIiMMPPzyqq7f/VMBBgwZtcwwAQFezW7ywYf369fHqq69GRERdXV2r++69995RU1MTa9asieXLl7f7czQ0NLT6+MqVK9t9LgCAtuwWEffGG2+Utmtra9vcvyniduRtRgYMGLBTa+sMm9esioYbxlV6GbuHojEizLyszLwyzL38zLzsNq9+rdJL2CG7RcStX7++tN2zZ8829+/Vq1dERKxbt67T1tQZ+vbtu2WjaEz3BzE9My8/M68Mcy8/My+7L3/5y5VeQrvsFhHXu3fv0vbGjRvb3H/Dhg0REdGnT592f462br2uXLkyhg8f3u7z7YzLLrssLr744q2uPNK51q1bF9XV1dGjR49KL2W3YeaVYe7lZ+aV8eUvf1nEdSWlK1QR7bpFumbNmoho363XJm09164cTj311Dj11FMrvQwAoAx2i1en9u7dO/bdd9+IaPsFCKtWrSpFXFd6nhsAQHO7RcRFROknMSxZsiQ2bdq03f0WL168zTEAAF3NbhNxxxxzTERsuVX6y1/+crv7zZ49u7Q9cuTITl8XAMDO2G0i7mMf+1hp+5Zbbmlxn8bGxvjBD34QERF77bVXHHfcceVYGgDADtttIm748OHxgQ98ICIibrrpppg7d+42+3z3u98t/ZSG888/3yuCAIAua7d4dWqT733vezFy5MhYt25dnHDCCfGNb3wjjjvuuFi3bl3cfvvtMW3atIiIGDhwYHzlK1+p8GoBALZvt4q4IUOGxB133BH/83/+z/jLX/4S3/jGN7bZZ+DAgTFz5syt3pYEAKCr2W1upzb56Ec/Gs8++2xccMEFMXDgwNhjjz1ir732iqFDh8bkyZNj4cKFcfjhh1d6mQAAraoqiqKo9CJ2Bw0NDaX3nVu+fHmXeHNgAKA8OqMDdrsrcQAAbwciDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgoepKL2B3sWnTptL2ypUrK7gSAKDcmv/d37wJdoWIK5NXXnmltD18+PAKrgQAqKRXXnklDj744F0+j9upAAAJVRVFUVR6EbuD9evXx3PPPRcREfvtt19UV3fsRdCVK1eWrvDNnz8/3vnOd3bo+dmWmZefmVeGuZefmZdfZ89806ZNpbty7373u6N37967fE63U8ukd+/eMWzYsLJ8rne+851RV1dXls/FFmZefmZeGeZefmZefp018464hdqc26kAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJOTNfgEAEnIlDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIRH3NvCHP/whLrzwwhg8eHDU1NREv379Yvjw4XHNNdfE2rVrK728FF5++eV44IEHYtKkSXHiiSfGvvvuG1VVVVFVVRVnnXXWDp9v1qxZMXbs2Kirq4tevXpFXV1djB07NmbNmtXxi0/s6aefjiuvvDJOPPHEGDBgQPTq1Stqa2tj4MCBcdZZZ8Xjjz++Q+cz99b95S9/idtvvz2+8pWvxKhRo+Lwww+PPffcM3r27Bn7779/HHvssTFlypR47bXX2nU+8951EydOLP2/pqqqKn7+85+3eYy5t0/zubb269hjj23zXF125gWpPfDAA8Wee+5ZRESLv4488shi6dKllV5ml7e9+UVEMW7cuHafp7GxsTjvvPNaPd95551XNDY2dt4Xk8QHP/jBVufU9OuTn/xksWHDhlbPZe7t8/DDD7dr5vvuu28xa9as7Z7HvDvGM888U1RXV281t8cee2y7+5v7jmnPn/WIKEaNGrXdc3T1mYu4xJ555plijz32KCKiqK2tLa644oriqaeeKh555JHi3HPPLf0BGzRoUPHGG29UerldWvP/IAcMGFCccMIJOxVx3/jGN0rHDRkypJgxY0Yxf/78YsaMGcWQIUNKj33zm9/svC8micMOO6yIiKJ///7F+eefX/z4xz8u5s+fX8ydO7e49tpriwMPPLA0rzPPPLPVc5l7+zz88MPFgAEDik996lPF9773veLuu+8u5s6dWzz55JPFHXfcUZx22mlF9+7di4goevbsWfzqV79q8Tzmves2b95cDBs2rIiIYv/9929XxJn7jmmax+c+97niueee2+6vF154Ybvn6OozF3GJHXvssUVEFNXV1cVTTz21zeNTpkwp/QH79re/XYEV5jFp0qTi/vvvL1566aWiKIpi2bJlOxxxzz//fOlf1UOHDi3Wrl271eNr1qwphg4dWvqeLVmypKO/jFROOumk4o477ig2bdrU4uOvvPJKMXDgwNL3Yc6cOS3uZ+7tt71ZN3fPPfeUZj527NhtHjfvjjF16tTSP7K//vWvtxlx5r7jmmZ6ySWX7NTxGWYu4pKaP39+6Q/oZz7zmRb32bx5czF48OAiIoq999672LhxY5lXmdfORNznP//50jFz585tcZ+5c+eW9vniF7/YgSt+e7r//vtL8/rSl77U4j7m3vEGDRpUuq36Vua96/7whz8UtbW1pWi75JJL2ow4c99xuxpxGWbuhQ1J3XvvvaXts88+u8V9unXrFp/61KciImLVqlXtesIsO6coirjvvvsiImLQoEExYsSIFvcbMWJEHHnkkRGx5XtYFEXZ1phR8yccL126dJvHzb1z1NTURETE+vXrt/q4eXeMz3/+87F69eoYN25cu55Ub+7ll2XmIi6pplft1dTUxPve977t7jdq1KjS9hNPPNHp69pdLVu2LFasWBERW8+8JU2PNzQ0xIsvvtjZS0tt48aNpe1u3bb935W5d7xFixbFM888ExFb/vJqzrx33Z133hkPPPBA9OvXL66++up2HWPu5Zdl5iIuqUWLFkVExOGHHx7V1dXb3a/5/4SbjqHjNZ/tW//ieyvfk/abPXt2abuluZp7x1i7dm08//zzce2118Zxxx0XmzdvjoiI888/f6v9zHvX/OlPfyrNdPLkybHffvu16zhz3zV33XVXHHnkkdGnT5/o27dvHHHEETFu3Lh47LHHtntMlplv/29/uqz169fHq6++GhERdXV1re679957R01NTaxZsyaWL19ejuXtlprPtq3vyYABA1o8jq01NjbGVVddVfr96aefvs0+5r7zbr311u0+FSMi4sILL4xPfOITW33MvHfNxIkT46WXXop/+Id/iHPOOafdx5n7rvnNb36z1e+XLFkSS5YsiR/84AfxsY99LG699dbYc889t9ony8xFXEJvvPFGabu2trbN/ZsibvXq1Z25rN3ajnxPmp5vFBG+J62YOnVqzJ8/PyIiTjnllBg6dOg2+5h7x3vve98bN954Y7z//e/f5jHz3nlPPPFETJ8+Paqrq+PGG2+Mqqqqdh9r7jtnjz32iDFjxsQ//uM/xqBBg6K2tjZeeeWVmD17dtx4443x2muvxb333hv19fXx8MMPR48ePUrHZpm5iEuo+ZONe/bs2eb+vXr1ioiIdevWddqadnc78j1p+n5E+J5sz+zZs+NrX/taRETsv//+8f3vf7/F/cx9533sYx8rhfG6deti6dKlceedd8Y999wTn/jEJ+K6666Lk08+eatjzHvnbNy4Mc4777woiiIuuOCCePe7371Dx5v7zlmxYkXstdde23z8+OOPj/Hjx8eJJ54YCxcujNmzZ8f3v//9+NKXvlTaJ8vMPScuod69e5e2mz/xe3s2bNgQERF9+vTptDXt7nbke9L0/YjwPWnJf/3Xf8Upp5wSmzZtil69esWdd94ZBxxwQIv7mvvO22uvveLv//7v4+///u9j2LBhccYZZ8Tdd98dP/jBD+KFF16I+vr6uPXWW7c6xrx3zpVXXhmLFi2Kv/u7v4tLLrlkh483953TUsA1OeCAA+LHP/5xKdCuv/76rR7PMnMRl1Dfvn1L2+25dLtmzZqIaN+tV3bOjnxPmr4fEb4nb7Vs2bI44YQTYtWqVdG9e/eYMWNGq68MM/eO98lPfjJOO+20aGxsjC9+8YuxatWq0mPmveMWL14c3/nOdyJiSyg0v/XWXubeOQ499NA4/vjjI2LL8+T++Mc/lh7LMnO3UxPq3bt37LvvvvHqq69GQ0NDq/uuWrWq9Aes+ZMv6VjNn/ja1vek+RNffU/+5o9//GN8+MMfjj/+8Y9RVVUVN998c5xyyimtHmPunaO+vj7uvPPOWLNmTTz00EPxL//yLxFh3jtj6tSpsXHjxjj00ENj7dq1cfvtt2+zz69//evS9qOPPhovvfRSRER89KMfjZqaGnPvREcddVTMnDkzIrbcfu3fv39E5PmzLuKSGjx4cDz++OOxZMmS2LRp03bfZmTx4sVbHUPnOOqoo0rbzWfeEt+Tbb366qtx/PHHxwsvvBARW65YNL1RdWvMvXM0f+uL3//+96Vt895xTbfaXnjhhTjzzDPb3P+yyy4rbS9btixqamrMvRNt7815s8zc7dSkjjnmmIjYchn3l7/85Xb3a/4+WyNHjuz0de2uDjnkkNK/4JrPvCVz5syJiIgDDzwwDj744M5eWpf35z//OT7ykY+U3gbgqquuii984QvtOtbcO0fTm5xGbH17yLwrw9w7T/O3H2macUSemYu4pD72sY+Vtm+55ZYW92lsbIwf/OAHEbHlCZ7HHXdcOZa2W6qqqor6+vqI2PKvsnnz5rW437x580r/aquvr9+htxl4O1q7dm2cdNJJ8fTTT0dExDe/+c346le/2u7jzb1z3HXXXaXt5q+kNO8dd+utt0ax5eeUb/dX8xc7PPbYY6WPNwWBuXeOF154IR5++OGI2PL8uAMPPLD0WJqZl/MHtdKxPvCBDxQRUVRXVxdPPfXUNo9PmTJll38A8O5q2bJlpdmNGzeuXcf89re/Laqrq4uIKIYOHVqsXbt2q8fXrl1bDB06tPQ9+93vftcJK89jw4YNxQknnFCa8/nnn79T5zH39rvllluKdevWtbrPtddeW/qeHHzwwcWbb7651ePm3fEuueSS0swfe+yxFvcx9x3zk5/8ZJs/u8299NJLxZAhQ0pz/+53v7vNPhlmLuISe/rpp4s+ffoUEVHU1tYWV155ZTF37tzi0UcfLc4777zSH86BAwcWf/nLXyq93C7t8ccfL2655ZbSr6uvvro0v5EjR2712C233LLd83zta18rHTdkyJDi9ttvL37xi18Ut99++1b/w/j6179evi+uixo7dmxpHh/60IeKZ599tnjuuee2++u3v/3tds9l7u1z0EEHFf369SvOPffc4rbbbiueeOKJ4plnnikef/zx4v/8n/9TjBw5sjSrnj17Fg8//HCL5zHvjtWeiCsKc98RBx10UNG/f/9i/PjxxY9+9KPiqaeeKhYuXFg8/PDDxTe/+c1in332Kc3rmGOOKdavX9/iebr6zEVccj/5yU+Kd7zjHaU/SG/9NXDgwOL555+v9DK7vHHjxm13hi392p7NmzcXn/70p1s99pxzzik2b95cxq+ua9qReUdEcdBBB233XObePgcddFC7Zl1XV1f87Gc/2+55zLtjtTfizL392vtn/eMf/3ixatWq7Z6nq89cxL0NvPjii8UFF1xQDBw4sNhjjz2Kvfbaqxg6dGgxefLkYs2aNZVeXgodFXFNZs6cWdTX1xf9+/cvevbsWfTv37+or68vHnzwwTJ8NTl0ZMQ1MffWLVmypLjxxhuLf/7nfy7e8573FAcccEBRXV1d1NbWFocddljx8Y9/vLjlllva/f8N8+4Y7Y24Jubetp///OfFt7/97eKf/umfioEDBxb9+vUrqquri7322qt497vfXXzmM59p8WlI29NVZ15VFNt5fS0AAF2WV6cCACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACT0/wNkgixrW/CnQwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 659, "width": 312 } }, "output_type": "display_data" } ], "source": [ "from spopt.region.base import infeasible_components\n", "\n", "gdf = geopandas.GeoDataFrame(\n", " geometry=geometries,\n", " data=numpy.ones((n_cols * n_rows + n_cols, 1), int),\n", " columns=[\"var\"],\n", ")\n", "gdf.plot(edgecolor=\"black\");" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:22.121043Z", "start_time": "2022-10-29T17:23:21.884332Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.885071Z", "iopub.status.busy": "2025-04-06T22:39:29.884988Z", "iopub.status.idle": "2025-04-06T22:39:29.933194Z", "shell.execute_reply": "2025-04-06T22:39:29.932961Z", "shell.execute_reply.started": "2025-04-06T22:39:29.885063Z" }, "tags": [] }, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = MaxP(gdf, w, \"var\", \"var\", 6, 2, policy=\"drop\")\n", "model.solve()\n", "ifcs = infeasible_components(gdf, w, \"var\", 6)\n", "keep_ids = numpy.where(~numpy.isin(w.component_labels, ifcs))[0]\n", "gdf[\"region\"] = -1\n", "gdf.loc[keep_ids, \"region\"] = model.labels_\n", "gdf.explore(column=\"region\", categorical=True) # areas in region -1 are dropped" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## No Feasible Components\n", "We also handle the case when none of the components are feasible by rasing an exception." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "ExecuteTime": { "end_time": "2022-10-29T17:23:22.135196Z", "start_time": "2022-10-29T17:23:22.124966Z" }, "execution": { "iopub.execute_input": "2025-04-06T22:39:29.933563Z", "iopub.status.busy": "2025-04-06T22:39:29.933498Z", "iopub.status.idle": "2025-04-06T22:39:29.936748Z", "shell.execute_reply": "2025-04-06T22:39:29.936552Z", "shell.execute_reply.started": "2025-04-06T22:39:29.933556Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No feasible components found in input.\n" ] } ], "source": [ "model = MaxP(gdf, w, \"var\", \"var\", 55, 2, policy=\"drop\")\n", "try:\n", " model.solve()\n", "except ValueError as e:\n", " if str(e).startswith(\"No feasible\"):\n", " print(e)\n", " else:\n", " raise e" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----------------" ] } ], "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.12.9" } }, "nbformat": 4, "nbformat_minor": 4 }