{ "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": "2023-12-10T18:30:11.076957Z", "iopub.status.busy": "2023-12-10T18:30:11.076532Z", "iopub.status.idle": "2023-12-10T18:30:11.115228Z", "shell.execute_reply": "2023-12-10T18:30:11.114103Z", "shell.execute_reply.started": "2023-12-10T18:30:11.076891Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: 2023-12-10T13:30:11.099582-05:00\n", "\n", "Python implementation: CPython\n", "Python version : 3.12.0\n", "IPython version : 8.18.0\n", "\n", "Compiler : Clang 15.0.7 \n", "OS : Darwin\n", "Release : 23.1.0\n", "Machine : x86_64\n", "Processor : i386\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": "2023-12-10T18:30:11.118847Z", "iopub.status.busy": "2023-12-10T18:30:11.118579Z", "iopub.status.idle": "2023-12-10T18:30:13.035696Z", "shell.execute_reply": "2023-12-10T18:30:13.034767Z", "shell.execute_reply.started": "2023-12-10T18:30:11.118819Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Watermark: 2.4.3\n", "\n", "libpysal : 4.9.2\n", "geopandas : 0.14.1\n", "numpy : 1.26.2\n", "shapely : 2.0.2\n", "spopt : 0.5.1.dev53+g5cadae7\n", "scipy : 1.11.4\n", "matplotlib: 3.8.2\n", "\n" ] } ], "source": [ "import geopandas\n", "import libpysal\n", "import matplotlib.pyplot as plt\n", "import numpy\n", "import scipy\n", "from scipy.spatial import KDTree\n", "import shapely\n", "from shapely.geometry import Polygon, box\n", "import spopt\n", "from spopt.region import MaxPHeuristic as MaxP\n", "import warnings\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": "2023-12-10T18:30:13.037354Z", "iopub.status.busy": "2023-12-10T18:30:13.036678Z", "iopub.status.idle": "2023-12-10T18:30:13.044805Z", "shell.execute_reply": "2023-12-10T18:30:13.043855Z", "shell.execute_reply.started": "2023-12-10T18:30:13.037328Z" } }, "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)]\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)\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": "2023-12-10T18:30:13.046881Z", "iopub.status.busy": "2023-12-10T18:30:13.046336Z", "iopub.status.idle": "2023-12-10T18:30:13.054633Z", "shell.execute_reply": "2023-12-10T18:30:13.054011Z", "shell.execute_reply.started": "2023-12-10T18:30:13.046857Z" } }, "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": "2023-12-10T18:30:13.055766Z", "iopub.status.busy": "2023-12-10T18:30:13.055531Z", "iopub.status.idle": "2023-12-10T18:30:13.242008Z", "shell.execute_reply": "2023-12-10T18:30:13.241186Z", "shell.execute_reply.started": "2023-12-10T18:30:13.055746Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAUnCAYAAAAy/TFcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAB7CAAAewgFu0HU+AABTjklEQVR4nO3de5xXdb3o//fAcHOGVLz9wmF7R7BdJx4B4cZC22lbNCZJ3bo7hebWriSW0lW04yVBEztuT24eeOs8Ci9tLylKudXACxwiMW0HJYjFEG5vVHIXZv3+oPk2yDAzwMz3O295Ph8PHo/FfNda85n3oLxY6/v9TlVRFEUAAJBKt0ovAACAHSfiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJFRd6QXsLtavXx/PPfdcRETst99+UV1t9ACwu9i0aVO88sorERHx7ne/O3r37r3L51QSZfLcc8/F8OHDK70MAKDC5s+fH8OGDdvl87idCgCQkCtxZbLffvuVtufPnx/vfOc7K7gaAKCcVq5cWboj17wJdoWIK5Pmz4F75zvfGXV1dRVcDQBQKR31vHi3UwEAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIQ6NeJefvnleOCBB2LSpElx4oknxr777htVVVVRVVUVZ511VrvOsX79+rjvvvti/Pjx8f73vz/69esXPXr0iH79+sXRRx8dl156aaxcubLda1q7dm1cffXVMXz48OjXr1/U1tbG4MGD48ILL4w//OEPO/mVAgCUV3VnnvyAAw7YpeOfffbZOOaYY+KNN97Y5rFVq1bFvHnzYt68eXHttdfG9OnT4/TTT2/1fEuXLo2TTjopfvvb32718cWLF8fixYtj+vTp8aMf/ShGjx69S+sGAOhsZbudOmDAgDjhhBN26Ji//OUvpYAbOXJkfOc734mHH344nn766fjpT38an/nMZ6J79+7xxhtvxL/8y7/EQw89tN1zrV69Ok4++eRSwJ177rnxyCOPxFNPPRVXXHFF1NbWxp///Oc47bTT4tlnn935LxQAoAw69UrcpEmTYtiwYTFs2LA44IAD4sUXX4xDDjmk3cd369YtTj/99LjkkkviqKOO2ubxE044IU488cQ45ZRTYvPmzTF+/Ph4/vnno6qqapt9r7nmmli8eHFEREyZMiUuuuii0mNHH310HHfccfHBD34w1q5dGxMmTIhHH310J75iAIDyqCqKoijXJ2secePGjYtbb721Q8576qmnxn/8x39ERMTTTz8dQ4YM2erxN998M/bff//405/+FIMHD45f//rX0a3bthchP/vZz8a///u/R0TEggUL4n3ve1+HrC8ioqGhIQYMGBAREcuXL4+6uroOOzcA0LV1Rge8LV6detxxx5W2ly5dus3jP//5z+NPf/pTRGyJx5YCLiK2erHF3Xff3aFrBADoSG+LiNuwYUNpu6VAe/zxx0vbo0aN2u55hg4dGjU1NRER8cQTT3TgCgEAOtbbIuJmz55d2h40aNA2jy9atKjVx5tUV1fHYYcdts0xAABdTae+sKEcfvWrX8XMmTMjIuJd73pXiy+AWL58eURE1NTUxF577dXq+QYMGBDPPvtsvPLKK7Fhw4bo1atXu9bR0NDQ6uM78l52AABtSR1xGzZsiH/913+NzZs3R0TElVde2eJ+TW9TUltb2+Y5m26nRmx5W5L2RlzTkxUr6a677opJkya1+L56dI5169bF6tWro7a2Nvr06VPp5ewWzLwyzL38zLwy+vbtG5dddlmceuqplV5Km1JH3Be/+MVYsGBBRGx5wcKYMWNa3G/9+vUREdGzZ882z9k82tatW9cBqyyfSZMmld5GhfJ6/fXXK72E3Y6ZV4a5l5+Zl9/FF18s4jrTd77znZg+fXpERLzvfe+LG264Ybv79u7dOyIiNm7c2OZ5m79IYkf+5dN0y3Z7Vq5cGcOHD2/3+XZG6QpcVbfoXrN3p34utti8+rUtG2ZeNmZeGeZefmZefpvXrIooGtPc0UoZcf/+7/8e3/jGNyIi4sgjj4yHHnpoq9ugb9W3b9+I2HJ7tC1r1qwpbbfn9muTrvS+b91r9o66L9xW6WXsFn4/ZUxE0WjmZWTmlWHu5Wfm5ddww7i/xXMC6V6dOmPGjPj85z8fEREHHXRQ/Od//mfst99+rR7TFFhr1qwpvV/c9jRdUdtvv/3a/Xw4AIBySxVxP/nJT+JTn/pUNDY2xjvf+c545JFH2nUFrPkrVlt7ztimTZtKbxY8ePDgXV8wAEAnSRNxjzzySJx++umxadOm2GeffeLhhx8uvadbW4455pjSdvP3lHurBQsWlG6njhw5ctcWDADQiVJE3FNPPRX19fWxYcOGeMc73hE//elP413vele7jz/22GNjzz33jIiI2267Lbb342Kb/yzXU045ZZfWDADQmbp8xD3zzDNx0kknxZo1a6KmpiYefPDBHf7B9D179owvfelLEbHlJzFcc8012+wzd+7cuOmmmyJiy4/mGjZs2K4vHgCgk3Tqq1OfeOKJWLJkSen3r776aml7yZIlW135itj6B9BHbPlh9h/5yEdKL0a4/PLLY88994xf//rX2/2c+++/f+y///7bfPyiiy6KO+64I373u9/FxIkTY8mSJXHGGWdEnz594rHHHosrr7wyNm3aFH369Inrrrtuh79WAIBy6tSImz59etx2W8svi37yySfjySef3Opjb424xx9/PF5++eXS7y+44II2P+cll1wSl1566TYf79u3b8ycOTNGjx4dzz//fEybNi2mTZu21T7veMc74oc//GG8973vbfPzAABUUpe/ndqRDj/88Fi4cGFMnjw5hg4dGnvttVfsscceceSRR8YFF1wQzz77bJx88smVXiYAQJs69Urcrbfeus0t0x1x1llnbXN1blfV1NTExIkTY+LEiR16XgCActqtrsQBALxdiDgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEqoqiqKo9CJ2Bw0NDTFgwICIiFi+fHnU1dV1+OfYZ5994vXXX4+o6hbda/bu8POzrc2rX9uyYeZlY+aVYe7lZ+bl1zTzfv36xWuvvdah5+6MDqje5TPQZaxevXrLRtH4t//4KQ8zLz8zrwxzLz8zL7vq6hx5lGOVtEttbe2WK3EREVXulJdF0fi3bTMvDzOvDHMvPzMvv7/OvEePHhVeSPuIuLeRPn36RERE99p9ou4Lt1V4NbuH308ZE1E0mnkZmXllmHv5mXn5NdwwLtVVT2kPAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJBQp0bcyy+/HA888EBMmjQpTjzxxNh3332jqqoqqqqq4qyzztrh882aNSvGjh0bdXV10atXr6irq4uxY8fGrFmz2n2OtWvXxtVXXx3Dhw+Pfv36RW1tbQwePDguvPDC+MMf/rDDawIAqITqzjz5AQcc0CHnKYoiPvvZz8a0adO2+viKFSvinnvuiXvuuSfOO++8uPHGG6Oqqmq751m6dGmcdNJJ8dvf/narjy9evDgWL14c06dPjx/96EcxevToDlk3AEBnKdvt1AEDBsQJJ5ywU8d+61vfKgXckCFDYsaMGTF//vyYMWNGDBkyJCIipk2bFhdffPF2z7F69eo4+eSTSwF37rnnxiOPPBJPPfVUXHHFFVFbWxt//vOf47TTTotnn312p9YJAFAunXolbtKkSTFs2LAYNmxYHHDAAfHiiy/GIYccskPnWLJkSUyZMiUiIoYOHRpz5syJPn36RETEsGHDYsyYMTFq1KhYsGBBTJ48Oc4+++w47LDDtjnPNddcE4sXL46IiClTpsRFF11Ueuzoo4+O4447Lj74wQ/G2rVrY8KECfHoo4/u7JcNANDpOvVK3Le//e04+eSTd+m26tSpU2PTpk0REXH99deXAq7JHnvsEddff31ERGzatCmuu+66bc7x5ptvxve+972IiBg8eHB85Stf2Wafo48+Os4555yIiHjsscfil7/85U6vGQCgs3XpV6cWRRH33XdfREQMGjQoRowY0eJ+I0aMiCOPPDIiIu69994oimKrx3/+85/Hn/70p4iIGDduXHTr1vKX3fzFFnffffcurh4AoPN06YhbtmxZrFixIiIiRo0a1eq+TY83NDTEiy++uNVjjz/++Db7tWTo0KFRU1MTERFPPPHEziwZAKAsOvU5cbtq0aJFpe1Bgwa1um/zxxctWrTVc+/ae57q6uo47LDD4tlnn93qmPZoaGho9fGVK1fu0PkAAFrTpSNu+fLlpe26urpW9x0wYECLxzX/fU1NTey1115tnufZZ5+NV155JTZs2BC9evVq11qbf/5KWbduXUREbF6zKhpuGFfh1ewmisaIMPOyMvPKMPfyM/Oy27z6tYj429+nXV2Xjrg33nijtF1bW9vqvk23QSO2vJ1IS+dp6xwtnae9EdcVlL7uorH0B5EyMfPyM/PKMPfyM/Oyq67u0nlU0qVXuX79+tJ2z549W923eWy9taCbztPWOdo6T2veevXvrVauXBnDhw9v9/l2Rm1tbbz++utbflPVpZ/u+Pbx138pR4SZl4uZV4a5l5+Zl99fZ96jR48KL6R9unTE9e7du7S9cePGVvfdsGFDafutb0PSdJ62ztHWeVrT1u3ecmhab/fafaLuC7dVeDW7h99PGRNRNJp5GZl5ZZh7+Zl5+TXcMC7VVc8unfZ9+/Ytbb/1FulbrVmzprT91tumTedp6xxtnQcAoKvo0hHX/OpWW6/+bH47860vMmg6z5o1a0rvF9fWefbbb79Uz4cDAHYvXTrijjrqqNJ204/M2p7mjw8ePHinzrNp06ZYunRpi+cAAOhKunTEHXLIIdG/f/+IiJg9e3ar+86ZMyciIg488MA4+OCDt3rsmGOOKW23dp4FCxaUbqeOHDlyZ5YMAFAWXTriqqqqor6+PiK2XEGbN29ei/vNmzevdIWtvr4+qqqqtnr82GOPjT333DMiIm677bZtfixXk1tvvbW0fcopp+zq8gEAOk2XjriIiAkTJpTer2X8+PHbvO3HunXrYvz48RGx5X1dJkyYsM05evbsGV/60pciYstPb7jmmmu22Wfu3Llx0003RcSWH801bNiwjvwyAAA6VKe+xcgTTzwRS5YsKf3+1VdfLW0vWbJkqytfEVv/APomAwcOjAsvvDCuuuqqWLBgQYwcOTK++tWvxmGHHRZLly6NyZMnx8KFCyMi4qKLLoojjjiixbVcdNFFcccdd8Tvfve7mDhxYixZsiTOOOOM6NOnTzz22GNx5ZVXxqZNm6JPnz5x3XXX7fLXDgDQmTo14qZPnx633dbye9s8+eST8eSTT271sZYiLiLiiiuuiJdffjluvvnmWLhwYZxxxhnb7HPOOefE5Zdfvt219O3bN2bOnBmjR4+O559/PqZNmxbTpk3bap93vOMd8cMf/jDe+973tv6FAQBUWJe/nRoR0a1bt7jpppti5syZUV9fH/3794+ePXtG//79o76+Ph588MGYPn16dOvW+pdz+OGHx8KFC2Py5MkxdOjQ2GuvvWKPPfaII488Mi644IJ49tln4+STTy7TVwUAsPM69Urcrbfeus0t010xevToGD169C6do6amJiZOnBgTJ07soFUBAJRfiitxAABsTcQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkFBVURRFpRexO2hoaIgBAwZERMTy5cujrq6uwz/HPvvsE6+//npEVbfoXrN3h5+fbW1e/dqWDTMvGzOvDHMvPzMvv6aZ9+vXL1577bUOPXdndED1Lp+BLmP16tVbNorGv/3HT3mYefmZeWWYe/mZedlVV+fIoxyrpF1qa2u3XImLiKhyp7wsisa/bZt5eZh5ZZh7+Zl5+f115j169KjwQtpHxL2N9OnTJyIiutfuE3VfuK3Cq9k9/H7KmIii0czLyMwrw9zLz8zLr+GGcamuekp7AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACaWJuI0bN8ZNN90U//RP/xTvfOc7o1evXlFbWxtHHnlkfPrTn4558+a16zyzZs2KsWPHRl1dXfTq1Svq6upi7NixMWvWrE7+CgAAOk51pRfQHsuXL4+TTjopnnvuua0+vnHjxvjd734Xv/vd7+KWW26JCy64IL773e9GVVXVNucoiiI++9nPxrRp07b6+IoVK+Kee+6Je+65J84777y48cYbWzweAKAr6fJX4jZt2rRVwL3nPe+JW2+9NebOnRs/+9nPYtKkSVFTUxMREVOnTo1rrrmmxfN861vfKgXckCFDYsaMGTF//vyYMWNGDBkyJCIipk2bFhdffHEZvioAgF3T5a/E3XfffaWAO/roo+Pxxx+P7t27lx4//vjjY8yYMXH00UfHm2++Gd/5znfiggsuiOrqv31pS5YsiSlTpkRExNChQ2POnDnRp0+fiIgYNmxYjBkzJkaNGhULFiyIyZMnx9lnnx2HHXZYGb9KAIAd0+WvxD355JOl7a9//etbBVyT973vfXHyySdHRMSqVati8eLFWz0+derU2LRpU0REXH/99aWAa7LHHnvE9ddfHxFbrvxdd911HfklAAB0uC4fcRs3bixtH3roodvdr/mVsw0bNpS2i6KI++67LyIiBg0aFCNGjGjx+BEjRsSRRx4ZERH33ntvFEWxS+sGAOhMXT7iBg4cWNp+4YUXtrvf0qVLIyKiqqoqjjjiiNLHly1bFitWrIiIiFGjRrX6uZoeb2hoiBdffHFnlwwA0Om6fMSdeeaZ8Y53vCMiIiZPnhybN2/eZp+FCxfGzJkzIyLijDPOKO0fEbFo0aLS9qBBg1r9XM0fb34cAEBX0+Vf2LDffvvFrbfeGp/4xCfiySefjGHDhsWECRNi4MCBsXr16njyySfju9/9bmzcuDHe+973xrXXXrvV8cuXLy9t19XVtfq5BgwY0OJx7dHQ0NDq4ytXrtyh8wEAtKbLR1xExCmnnBILFiyIa6+9Nm6++eYYN27cVo8fcMAB8e1vfzvOO++80tuNNHnjjTdK27W1ta1+nubHrl69eofW2DwAK2XdunUREbF5zapouGFcG3vTIYrGiDDzsjLzyjD38jPzstu8+rWI+Nvfp11dioh7880340c/+lHcf//9Lb7g4L//+79jxowZMXDgwDjppJO2emz9+vWl7Z49e7b6eXr16lXazvINbK4UnkVj6Q8iZWLm5WfmlWHu5WfmZdf8bcq6si6/yjVr1sTo0aNjzpw50b1795g4cWKcffbZceihh8b69evj//2//xf/63/9r3jiiSfiox/9aEydOjXOP//80vG9e/cubTd/pWtLmr+q9a1vQ9KWtm6/rly5MoYPH75D59xRtbW18frrr2/5TVWXf7rj28Nf/6UcEWZeLmZeGeZefmZefn+deY8ePSq8kPbp8hF3ySWXxJw5cyIi4qabbtrqVmrPnj3j+OOPj+OOOy5OOOGEeOyxx+LLX/5yHHfccfGe97wnIiL69u1b2r+tW6Rr1qwpbbd16/Wt2nq+XTk0hWf32n2i7gu3VXg1u4ffTxkTUTSaeRmZeWWYe/mZefk13DAu1VXPLp32RVHELbfcEhFb3mrkrc+Fa1JdXR2XXXZZREQ0NjaWjonYOq7aevFB86tpXeE5bgAA29OlI+6///u/S7cHm36+6fa8733vK203/4kNRx11VIsfb0nzxwcPHrxDawUAKKcuHXHNn1jY9GOztufNN99s8bhDDjkk+vfvHxERs2fPbvUcTbdtDzzwwDj44IN3dLkAAGXTpSOuX79+pTfunTt3bqsh1zzQDjnkkNJ2VVVV1NfXR8SWK23z5s1r8fh58+aVrsTV19dHVVXVLq8fAKCzdOmI69atW+ktQ/74xz/GFVdc0eJ+q1atiq9+9aul35988slbPT5hwoTS1bnx48dv8/Yh69ati/Hjx0fElqt4EyZM6KgvAQCgU3TpiIuImDRpUuyxxx4REXHppZfGmDFj4j/+4z9i4cKFMXfu3Jg6dWq8973vjd/85jcREfGP//iPccIJJ2x1joEDB8aFF14YERELFiyIkSNHxh133BELFiyIO+64I0aOHBkLFiyIiIiLLrpoq5+9CgDQFXX5txgZNGhQ3HfffXHmmWfGq6++Gvfff3/cf//9Le77oQ99KO66664WH7viiivi5ZdfjptvvjkWLlwYZ5xxxjb7nHPOOXH55Zd36PoBADpDl78SFxHx4Q9/OBYvXhyTJ0+OY489Nvbbb7/o0aNH9OnTJw455JA4/fTT4957743//M//jL333rvFc3Tr1i1uuummmDlzZtTX10f//v2jZ8+e0b9//6ivr48HH3wwpk+fHt26pRgJALCb6/JX4prss88+MXHixJg4ceIunWf06NExevToDloVAEBluOwEAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgoaqiKIpKL2J30NDQEAMGDIiIiOXLl0ddXV2Hf4599tknXn/99YiqbtG9Zu8OPz/b2rz6tS0bZl42Zl4Z5l5+Zl5+TTPv169fvPbaax167s7ogOpdPgNdxurVq7dsFI1/+4+f8jDz8jPzyjD38jPzsquuzpFHOVZJu9TW1m65EhcRUeVOeVkUjX/bNvPyMPPKMPfyM/Py++vMe/ToUeGFtI+Iexvp06dPRER0r90n6r5wW4VXs3v4/ZQxEUWjmZeRmVeGuZefmZdfww3jUl31lPYAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmli7hXX301pkyZEiNHjoz/7//7/6JXr17Rv3//eP/73x8XXXRRzJ07t81zzJo1K8aOHRt1dXXRq1evqKuri7Fjx8asWbPK8BUAAOy66kovYEfcdddd8bnPfS5ee+21rT6+cuXKWLlyZcyfPz+ef/75uPfee1s8viiK+OxnPxvTpk3b6uMrVqyIe+65J+65554477zz4sYbb4yqqqrO+jIAAHZZmoj7wQ9+EGeffXY0NjbG/vvvH5/73OfimGOOiX79+sVLL70US5cujfvvvz969Oix3XN861vfKgXckCFDYuLEiXHYYYfF0qVLY8qUKbFw4cKYNm1a7LfffnH55ZeX60sDANhhKSJu0aJFcd5550VjY2N84AMfiPvvvz/23HPPbfYbP358bNy4scVzLFmyJKZMmRIREUOHDo05c+ZEnz59IiJi2LBhMWbMmBg1alQsWLAgJk+eHGeffXYcdthhnfdFAQDsghTPiRs/fnxs2LAh9t1337j77rtbDLgmPXv2bPHjU6dOjU2bNkVExPXXX18KuCZ77LFHXH/99RERsWnTprjuuus6ZvEAAJ2gy0fc4sWL45FHHomIiC9+8Yux77777vA5iqKI++67LyIiBg0aFCNGjGhxvxEjRsSRRx4ZERH33ntvFEWxk6sGAOhcXT7i7rrrrtL2aaedVtpetWpVPP/889u8yKEly5YtixUrVkRExKhRo1rdt+nxhoaGePHFF3dixQAAna/LR9y8efMiImLPPfeMwYMHxw9/+MP4H//jf0S/fv1i4MCBse+++8ahhx4a3/72t2P16tUtnmPRokWl7UGDBrX6+Zo/3vw4AICupMu/sOE3v/lNREQcfPDBMX78+Ljhhhu22WfZsmVx6aWXxo9//OP46U9/Gv3799/q8eXLl5e26+rqWv18AwYMaPG4tjQ0NLT6+MqVK9t9LgCAtnT5iHv99dcjYstz4371q1/FXnvtFVdddVWMHTs23vGOd8Rzzz0XkyZNioceeih+/etfx2mnnRaPP/54dOv2t4uMb7zxRmm7tra21c9XU1NT2t7elb2WNI+/Slm3bl1ERGxesyoabhhX4dXsJorGiDDzsjLzyjD38jPzstu8estTtJr+Pu3qunzErVmzJiIiNmzYEN27d4+HHnpoqxcmDB06NB544IE4+eST46GHHoqnnnoq7r777jj11FNL+6xfv760vb1Xrzbp1atXaTvLN7FJKTqLxtIfRMrEzMvPzCvD3MvPzMuuurrL51FEJIi43r17l0LutNNOa/GVpd26dYurr746HnrooYiImDFjxlYR17t379L29t5HrsmGDRtK2299G5LWtHXrdeXKlTF8+PB2n29n1NbWlq5cRlWXf7rj28Nf/6UcEWZeLmZeGeZefmZefn+deWs/OKAr6fIR17dv31LEnXjiidvd713velcceOCBsWLFivjFL36xzTmatHWLtOlzRbR967W5tp5rVw5N0dm9dp+o+8JtFV7N7uH3U8ZEFI1mXkZmXhnmXn5mXn4NN4xLddWzy6d98+eatfdFCS+//PJWH29+XFsvQGh+Ra0rPM8NAKAlXT7i3vWud5W2N2/e3Oq+TY+/9V72UUcdVdpevHhxq+do/vjgwYPbvU4AgHLq8hH3wQ9+sLS9dOnSVvd94YUXIiLiwAMP3OrjhxxySOltR2bPnt3qOebMmVM6x8EHH7yjywUAKIsuH3FjxowpPcHw7rvv3u5+s2fPLv30hg984ANbPVZVVRX19fURseVKW9MbCL/VvHnzSlfi6uvro6qqapfXDwDQGbp8xO2zzz7xr//6rxER8fDDD8ftt9++zT5vvPFGTJgwofT7z3zmM9vsM2HChNJt1vHjx2/z9iHr1q2L8ePHR8SW27HNzwcA0NV0+YiLiPj2t78df/d3fxcREZ/85Cdj/Pjx8dhjj8Uvf/nLuPXWW2P48OHxzDPPRETE5z73uRg2bNg25xg4cGBceOGFERGxYMGCGDlyZNxxxx2xYMGCuOOOO2LkyJGxYMGCiIi46KKL4ogjjijPFwcAsBO6/FuMRETst99+MWvWrBgzZkwsWbIk/u3f/i3+7d/+bZv9Pv3pT8f3vve97Z7niiuuiJdffjluvvnmWLhwYZxxxhnb7HPOOefE5Zdf3qHrBwDoaCmuxEVseaXoM888E1dffXW8//3vj379+kXPnj2jrq4u/vmf/zkeffTRuOmmm1p9g75u3brFTTfdFDNnzoz6+vro379/9OzZM/r37x/19fXx4IMPxvTp07f6kV0AAF1RiitxTWpqauLCCy8s3RbdWaNHj47Ro0d30KoAAMrPJScAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJVRVFUVR6EbuDhoaGGDBgQERELF++POrq6jr8c+yzzz7x+uuvR1R1i+41e3f4+dnW5tWvbdkw87Ix88ow9/Iz8/Jrmnm/fv3itdde69Bzd0YHVO/yGegyVq9evWWjaPzbf/yUh5mXn5lXhrmXn5mXXXV1jjzKsUrapba2dsuVuIiIKnfKy6Jo/Nu2mZeHmVeGuZefmZffX2feo0ePCi+kfUTc20ifPn0iIqJ77T5R94XbKrya3cPvp4yJKBrNvIzMvDLMvfzMvPwabhiX6qqntAcASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASCh1xE2cODGqqqpKv37+85+3ecysWbNi7NixUVdXF7169Yq6uroYO3ZszJo1q/MXDADQQdJG3K9+9auYOnVqu/cviiI+85nPxIknnhj33HNPrFixIjZu3BgrVqyIe+65J0488cT4zGc+E0VRdOKqAQA6RsqIa2xsjHPPPTc2bdoU+++/f7uO+da3vhXTpk2LiIghQ4bEjBkzYv78+TFjxowYMmRIRERMmzYtLr744k5bNwBAR0kZcf/7f//v+MUvfhGDBg2Kc845p839lyxZElOmTImIiKFDh8aTTz4ZZ5xxRgwbNizOOOOMeOKJJ2Lo0KERETF58uRYunRpp64fAGBXpYu45cuXl66Wff/734+ePXu2eczUqVNj06ZNERFx/fXXR58+fbZ6fI899ojrr78+IiI2bdoU1113XccuGgCgg6WLuM9//vOxevXqGDduXBx77LFt7l8URdx3330RETFo0KAYMWJEi/uNGDEijjzyyIiIuPfeez03DgDo0lJF3J133hkPPPBA9OvXL66++up2HbNs2bJYsWJFRESMGjWq1X2bHm9oaIgXX3xxl9YKANCZqiu9gPb605/+FOeff35EbHne2n777deu4xYtWlTaHjRoUKv7Nn980aJFccghh7R7fQ0NDa0+vnLlynafCwCgLWkibuLEifHSSy/FP/zDP7TrxQxNli9fXtquq6trdd8BAwa0eFx7ND+2UtatWxcREZvXrIqGG8ZVeDW7iaIxIsy8rMy8Msy9/My87Davfi0i/vb3aVeXIuKeeOKJmD59elRXV8eNN94YVVVV7T72jTfeKG3X1ta2um9NTU1pe/Xq1Tu+0AorrbloLP1BpEzMvPzMvDLMvfzMvOyqq1PkUdePuI0bN8Z5550XRVHEBRdcEO9+97t36Pj169eXttt6JWuvXr1K2zta4W1duVu5cmUMHz58h865o2pra+P111/f8puqVE93zOuv/1KOCDMvFzOvDHMvPzMvv7/OvEePHhVeSPt0+Yi78sorY9GiRfF3f/d3cckll+zw8b179y5tb9y4sdV9N2zYUNp+69uQtKWtW7Xl0LTm7rX7RN0XbqvwanYPv58yJqJoNPMyMvPKMPfyM/Pya7hhXKqrnl067RcvXhzf+c53ImLL+7s1v93ZXn379i1tt3WLdM2aNaXttm69AgBUUpe+Ejd16tTYuHFjHHroobF27dq4/fbbt9nn17/+dWn70UcfjZdeeikiIj760Y9GTU3NVlfI2noFafNbol3hhQoAANvTpSOu6fbmCy+8EGeeeWab+1922WWl7WXLlkVNTU0cddRRpY8tXry41eObPz548OAdXS4AQNl06dupHeGQQw6J/v37R0TE7NmzW913zpw5ERFx4IEHxsEHH9zZSwMA2GldOuJuvfXWKIqi1V/NX+zw2GOPlT7eFGFVVVVRX18fEVuutM2bN6/FzzVv3rzSlbj6+vodehsTAIBy69IR11EmTJhQes+X8ePHb/P2IevWrYvx48dHxJb3hpkwYUK5lwgAsEN2i4gbOHBgXHjhhRERsWDBghg5cmTccccdsWDBgrjjjjti5MiRsWDBgoiIuOiii+KII46o5HIBANrUpV/Y0JGuuOKKePnll+Pmm2+OhQsXxhlnnLHNPuecc05cfvnlFVgdAMCO2S2uxEVEdOvWLW666aaYOXNm1NfXR//+/aNnz57Rv3//qK+vjwcffDCmT58e3brtNiMBABJLfyXu0ksvjUsvvbTd+48ePTpGjx7deQsCACgDl50AABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkVFUURVHpRewOGhoaYsCAARERsXz58qirq+vwz7HPPvvE66+/HlHVLbrX7N3h52dbm1e/tmXDzMvGzCvD3MvPzMuvaeb9+vWL1157rUPP3RkdUL3LZ6DLWL169ZaNovFv//FTHmZefmZeGeZefmZedtXVOfIoxyppl9ra2i1X4iIiqtwpL4ui8W/bZl4eZl4Z5l5+Zl5+f515jx49KryQ9hFxbyN9+vSJiIjutftE3Rduq/Bqdg+/nzImomg08zIy88ow9/Iz8/JruGFcqque0h4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCKSLu6aefjiuvvDJOPPHEGDBgQPTq1Stqa2tj4MCBcdZZZ8Xjjz++Q+ebNWtWjB07Nurq6qJXr15RV1cXY8eOjVmzZnXSVwAA0LGqK72AtowaNSrmzJmzzcc3btwYzz//fDz//PNx2223xSc/+cmYPn169OzZc7vnKooiPvvZz8a0adO2+viKFSvinnvuiXvuuSfOO++8uPHGG6OqqqrDvxYAgI7S5a/ErVixIiIi+vfvH+eff378+Mc/jvnz58fcuXPj2muvjQMPPDAiIv7v//2/cdZZZ7V6rm9961ulgBsyZEjMmDEj5s+fHzNmzIghQ4ZERMS0adPi4osv7rwvCACgA3T5K3GDBg2KK6+8Mj7+8Y9H9+7dt3psxIgR8clPfjJGjhwZv/vd72LGjBnxuc99Lj7wgQ9sc54lS5bElClTIiJi6NChMWfOnOjTp09ERAwbNizGjBkTo0aNigULFsTkyZPj7LPPjsMOO6zzv0AAgJ3Q5a/EPfDAA3H66advE3BN9t133/jud79b+v2Pf/zjFvebOnVqbNq0KSIirr/++lLANdljjz3i+uuvj4iITZs2xXXXXdcBqwcA6BxdPuLa49hjjy1tL126dJvHi6KI++67LyK2XNkbMWJEi+cZMWJEHHnkkRERce+990ZRFB2/WACADvC2iLiNGzeWtrt12/ZLWrZsWem5daNGjWr1XE2PNzQ0xIsvvthxiwQA6EBvi4ibPXt2aXvQoEHbPL5o0aJWH2+u+ePNjwMA6Eq6/Asb2tLY2BhXXXVV6fenn376NvssX768tF1XV9fq+QYMGNDicW1paGho9fGVK1e2+1wAAG1JH3FTp06N+fPnR0TEKaecEkOHDt1mnzfeeKO0XVtb2+r5ampqSturV69u9zqax1+lbV6zKhpuGFfpZeweisaIMPOyMvPKMPfyM/Oy27z6tUovYYekjrjZs2fH1772tYiI2H///eP73/9+i/utX7++tN3amwFHRPTq1au0vW7dug5YZfn07dt3y0bRmO4PYnpmXn5mXhnmXn5mXnZf/vKXK72Edkkbcf/1X/8Vp5xySmzatCl69eoVd955ZxxwwAEt7tu7d+/SdvMXQbRkw4YNpe23vg1Ja9q69bpy5coYPnx4u8+3My677LK4+OKLt7rySOdat25dVFdXR48ePSq9lN2GmVeGuZefmVfGl7/8ZRHXmZYtWxYnnHBCrFq1Krp37x4zZsxo9VWnpStU0fYt0jVr1pS227r12lxbz7Urh1NPPTVOPfXUSi8DACiDdK9O/eMf/xgf/vCH449//GNUVVXFzTffHKecckqrxzQPrLZegND8ilpXep4bAEBzqSLu1VdfjeOPPz5eeOGFiNjykxc+9alPtXncUUcdVdpevHhxq/s2f3zw4ME7uVIAgM6VJuL+/Oc/x0c+8pH4zW9+ExERV111VXzhC19o17GHHHJI9O/fPyK2fk+5lsyZMyciIg488MA4+OCDd37BAACdKEXErV27Nk466aR4+umnIyLim9/8Znz1q19t9/FVVVVRX18fEVuutM2bN6/F/ebNm1e6EldfXx9VVVW7uHIAgM7R5SNu48aNccopp8STTz4ZERHnn39+XH755Tt8ngkTJkR19ZbXcYwfP36btw9Zt25djB8/PiIiqqurY8KECbu2cACATtTlX5165plnxs9+9rOIiPjQhz4U55xzTvz617/e7v49e/aMgQMHbvPxgQMHxoUXXhhXXXVVLFiwIEaOHBlf/epX47DDDoulS5fG5MmTY+HChRERcdFFF8URRxzROV8QAEAHqCqKoqj0Ilqzo7c0DzrooO3+4PrGxsY499xz4+abb97u8eecc05MmzYtunXr2IuUDQ0NpVe7Ll++vEu8JQkAUB6d0QFd/nZqR+rWrVvcdNNNMXPmzKivr4/+/ftHz549o3///lFfXx8PPvhgTJ8+vcMDDgCgo3X526mdcaFw9OjRMXr06A4/LwBAubjkBACQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AIKHqSi+AjnPXXXfFpEmT4o033qj0UnYb69ati9WrV0dtbW306dOn0svZLZh5ZZh7+Zl5ZfTt2zcuu+yyOPXUUyu9lDZVFUVRVHoRu4OGhoYYMGBAREQsX7486urqOvxzDB48OBYvXtzh5wWA3cmgQYNi0aJFHXrOzugAV+LeRkpX4Kq6RfeavSu7mN3E5tWvbdkw87Ix88ow9/Iz8/LbvGZVRNGY5o6WiHsb6l6zd9R94bZKL2O38PspYyKKRjMvIzOvDHMvPzMvv4Ybxv0tnhPwwgYAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIR2y4j7wx/+EBdeeGEMHjw4ampqol+/fjF8+PC45pprYu3atZVeHgBAm6orvYBymzlzZnziE5+IP//5z6WPrV27Nn7xi1/EL37xi5g+fXo8+OCDceihh1ZwlQAArdutrsT96le/itNPPz3+/Oc/R21tbVxxxRXx1FNPxSOPPBLnnntuRET89re/jZNOOilWr15d4dUCAGzfbnUlbsKECbF27dqorq6On/3sZ3H00UeXHvvQhz4URxxxREycODEWL14c1157bUyaNKmCqwUA2L7d5krcL37xi/j5z38eERHnnHPOVgHX5Ctf+UoMHjw4IiKuu+66ePPNN8u5RACAdtttIu7ee+8tbZ999tkt7tOtW7f41Kc+FRERq1atKkUfAEBXs9tE3OOPPx4RETU1NfG+971vu/uNGjWqtP3EE090+roAAHbGbhNxixYtioiIww8/PKqrt/9UwEGDBm1zDABAV7NbvLBh/fr18eqrr0ZERF1dXav77r333lFTUxNr1qyJ5cuXt/tzNDQ0tPr4ypUr230uAIC27BYR98Ybb5S2a2tr29y/KeJ25G1GBgwYsFNr6wyb16yKhhvGVXoZu4eiMSLMvKzMvDLMvfzMvOw2r36t0kvYIbtFxK1fv7603bNnzzb379WrV0RErFu3rtPW1Bn69u27ZaNoTPcHMT0zLz8zrwxzLz8zL7svf/nLlV5Cu+wWEde7d+/S9saNG9vcf8OGDRER0adPn3Z/jrZuva5cuTKGDx/e7vPtjMsuuywuvvjira480rnWrVsX1dXV0aNHj0ovZbdh5pVh7uVn5pXx5S9/WcR1JaUrVBHtukW6Zs2aiGjfrdcmbT3XrhxOPfXUOPXUUyu9DACgDHaLV6f27t079t1334ho+wUIq1atKkVcV3qeGwBAc7tFxEVE6ScxLFmyJDZt2rTd/RYvXrzNMQAAXc1uE3HHHHNMRGy5VfrLX/5yu/vNnj27tD1y5MhOXxcAwM7YbSLuYx/7WGn7lltuaXGfxsbG+MEPfhAREXvttVccd9xx5VgaAMAO220ibvjw4fGBD3wgIiJuuummmDt37jb7fPe73y39lIbzzz/fK4IAgC5rt3h1apPvfe97MXLkyFi3bl2ccMIJ8Y1vfCOOO+64WLduXdx+++0xbdq0iIgYOHBgfOUrX6nwagEAtm+3irghQ4bEHXfcEf/zf/7P+Mtf/hLf+MY3ttln4MCBMXPmzK3elgQAoKvZbW6nNvnoRz8azz77bFxwwQUxcODA2GOPPWKvvfaKoUOHxuTJk2PhwoVx+OGHV3qZAACtqiqKoqj0InYHDQ0NpfedW758eZd4c2AAoDw6owN2uytxAABvByIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACCh6kovYHexadOm0vbKlSsruBIAoNya/93fvAl2hYgrk1deeaW0PXz48AquBACopFdeeSUOPvjgXT6P26kAAAlVFUVRVHoRu4P169fHc889FxER++23X1RXd+xF0JUrV5au8M2fPz/e+c53duj52ZaZl5+ZV4a5l5+Zl19nz3zTpk2lu3Lvfve7o3fv3rt8TrdTy6R3794xbNiwsnyud77znVFXV1eWz8UWZl5+Zl4Z5l5+Zl5+nTXzjriF2pzbqQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAk5M1+AQASciUOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEfc28Ic//CEuvPDCGDx4cNTU1ES/fv1i+PDhcc0118TatWsrvbwUXn755XjggQdi0qRJceKJJ8a+++4bVVVVUVVVFWedddYOn2/WrFkxduzYqKuri169ekVdXV2MHTs2Zs2a1fGLT+zpp5+OK6+8Mk488cQYMGBA9OrVK2pra2PgwIFx1llnxeOPP75D5zP31v3lL3+J22+/Pb7yla/EqFGj4vDDD48999wzevbsGfvvv38ce+yxMWXKlHjttdfadT7z3nUTJ04s/b+mqqoqfv7zn7d5jLm3T/O5tvbr2GOPbfNcXXbmBak98MADxZ577llERIu/jjzyyGLp0qWVXmaXt735RUQxbty4dp+nsbGxOO+881o933nnnVc0NjZ23heTxAc/+MFW59T065Of/GSxYcOGVs9l7u3z8MMPt2vm++67bzFr1qztnse8O8YzzzxTVFdXbzW3xx57bLv7m/uOac+f9YgoRo0atd1zdPWZi7jEnnnmmWKPPfYoIqKora0trrjiiuKpp54qHnnkkeLcc88t/QEbNGhQ8cYbb1R6uV1a8/8gBwwYUJxwwgk7FXHf+MY3SscNGTKkmDFjRjF//vxixowZxZAhQ0qPffOb3+y8LyaJww47rIiIon///sX5559f/PjHPy7mz59fzJ07t7j22muLAw88sDSvM888s9VzmXv7PPzww8WAAQOKT33qU8X3vve94u677y7mzp1bPPnkk8Udd9xRnHbaaUX37t2LiCh69uxZ/OpXv2rxPOa96zZv3lwMGzasiIhi//33b1fEmfuOaZrH5z73ueK5557b7q8XXnhhu+fo6jMXcYkde+yxRUQU1dXVxVNPPbXN41OmTCn9Afv2t79dgRXmMWnSpOL+++8vXnrppaIoimLZsmU7HHHPP/986V/VQ4cOLdauXbvV42vWrCmGDh1a+p4tWbKko7+MVE466aTijjvuKDZt2tTi46+88koxcODA0vdhzpw5Le5n7u23vVk3d88995RmPnbs2G0eN++OMXXq1NI/sr/+9a+3GXHmvuOaZnrJJZfs1PEZZi7ikpo/f37pD+hnPvOZFvfZvHlzMXjw4CIiir333rvYuHFjmVeZ185E3Oc///nSMXPnzm1xn7lz55b2+eIXv9iBK357uv/++0vz+tKXvtTiPube8QYNGlS6rfpW5r3r/vCHPxS1tbWlaLvkkkvajDhz33G7GnEZZu6FDUnde++9pe2zzz67xX26desWn/rUpyIiYtWqVe16wiw7pyiKuO+++yIiYtCgQTFixIgW9xsxYkQceeSREbHle1gURdnWmFHzJxwvXbp0m8fNvXPU1NRERMT69eu3+rh5d4zPf/7zsXr16hg3bly7nlRv7uWXZeYiLqmmV+3V1NTE+973vu3uN2rUqNL2E0880enr2l0tW7YsVqxYERFbz7wlTY83NDTEiy++2NlLS23jxo2l7W7dtv3flbl3vEWLFsUzzzwTEVv+8mrOvHfdnXfeGQ888ED069cvrr766nYdY+7ll2XmIi6pRYsWRUTE4YcfHtXV1dvdr/n/hJuOoeM1n+1b/+J7K9+T9ps9e3Zpu6W5mnvHWLt2bTz//PNx7bXXxnHHHRebN2+OiIjzzz9/q/3Me9f86U9/Ks108uTJsd9++7XrOHPfNXfddVcceeSR0adPn+jbt28cccQRMW7cuHjssce2e0yWmW//b3+6rPXr18err74aERF1dXWt7rv33ntHTU1NrFmzJpYvX16O5e2Wms+2re/JgAEDWjyOrTU2NsZVV11V+v3pp5++zT7mvvNuvfXW7T4VIyLiwgsvjE984hNbfcy8d83EiRPjpZdein/4h3+Ic845p93Hmfuu+c1vfrPV75csWRJLliyJH/zgB/Gxj30sbr311thzzz232ifLzEVcQm+88UZpu7a2ts39myJu9erVnbms3dqOfE+anm8UEb4nrZg6dWrMnz8/IiJOOeWUGDp06Db7mHvHe+973xs33nhjvP/979/mMfPeeU888URMnz49qqur48Ybb4yqqqp2H2vuO2ePPfaIMWPGxD/+4z/GoEGDora2Nl555ZWYPXt23HjjjfHaa6/FvffeG/X19fHwww9Hjx49SsdmmbmIS6j5k4179uzZ5v69evWKiIh169Z12pp2dzvyPWn6fkT4nmzP7Nmz42tf+1pEROy///7x/e9/v8X9zH3nfexjHyuF8bp162Lp0qVx5513xj333BOf+MQn4rrrrouTTz55q2PMe+ds3LgxzjvvvCiKIi644IJ497vfvUPHm/vOWbFiRey1117bfPz444+P8ePHx4knnhgLFy6M2bNnx/e///340pe+VNony8w9Jy6h3r17l7abP/F7ezZs2BAREX369Om0Ne3uduR70vT9iPA9acl//dd/xSmnnBKbNm2KXr16xZ133hkHHHBAi/ua+87ba6+94u///u/j7//+72PYsGFxxhlnxN133x0/+MEP4oUXXoj6+vq49dZbtzrGvHfOlVdeGYsWLYq/+7u/i0suuWSHjzf3ndNSwDU54IAD4sc//nEp0K6//vqtHs8ycxGXUN++fUvb7bl0u2bNmoho361Xds6OfE+avh8RvidvtWzZsjjhhBNi1apV0b1795gxY0arrwwz9473yU9+Mk477bRobGyML37xi7Fq1arSY+a94xYvXhzf+c53ImJLKDS/9dZe5t45Dj300Dj++OMjYsvz5P74xz+WHssyc7dTE+rdu3fsu+++8eqrr0ZDQ0Or+65atar0B6z5ky/pWM2f+NrW96T5E199T/7mj3/8Y3z4wx+OP/7xj1FVVRU333xznHLKKa0eY+6do76+Pu68885Ys2ZNPPTQQ/Ev//IvEWHeO2Pq1KmxcePGOPTQQ2Pt2rVx++23b7PPr3/969L2o48+Gi+99FJERHz0ox+Nmpoac+9ERx11VMycOTMittx+7d+/f0Tk+bMu4pIaPHhwPP7447FkyZLYtGnTdt9mZPHixVsdQ+c46qijStvNZ94S35Ntvfrqq3H88cfHCy+8EBFbrlg0vVF1a8y9czR/64vf//73pW3z3nFNt9peeOGFOPPMM9vc/7LLLittL1u2LGpqasy9E23vzXmzzNzt1KSOOeaYiNhyGfeXv/zldvdr/j5bI0eO7PR17a4OOeSQ0r/gms+8JXPmzImIiAMPPDAOPvjgzl5al/fnP/85PvKRj5TeBuCqq66KL3zhC+061tw7R9ObnEZsfXvIvCvD3DtP87cfaZpxRJ6Zi7ikPvaxj5W2b7nllhb3aWxsjB/84AcRseUJnscdd1w5lrZbqqqqivr6+ojY8q+yefPmtbjfvHnzSv9qq6+v36G3GXg7Wrt2bZx00knx9NNPR0TEN7/5zfjqV7/a7uPNvXPcddddpe3mr6Q07x136623RrHl55Rv91fzFzs89thjpY83BYG5d44XXnghHn744YjY8vy4Aw88sPRYmpmX8we10rE+8IEPFBFRVFdXF0899dQ2j0+ZMmWXfwDw7mrZsmWl2Y0bN65dx/z2t78tqquri4gohg4dWqxdu3arx9euXVsMHTq09D373e9+1wkrz2PDhg3FCSecUJrz+eefv1PnMff2u+WWW4p169a1us+1115b+p4cfPDBxZtvvrnV4+bd8S655JLSzB977LEW9zH3HfOTn/xkmz+7zb300kvFkCFDSnP/7ne/u80+GWYu4hJ7+umniz59+hQRUdTW1hZXXnllMXfu3OLRRx8tzjvvvNIfzoEDBxZ/+ctfKr3cLu3xxx8vbrnlltKvq6++ujS/kSNHbvXYLbfcst3zfO1rXysdN2TIkOL2228vfvGLXxS33377Vv/D+PrXv16+L66LGjt2bGkeH/rQh4pnn322eO6557b767e//e12z2Xu7XPQQQcV/fr1K84999zitttuK5544onimWeeKR5//PHi//yf/1OMHDmyNKuePXsWDz/8cIvnMe+O1Z6IKwpz3xEHHXRQ0b9//2L8+PHFj370o+Kpp54qFi5cWDz88MPFN7/5zWKfffYpzeuYY44p1q9f3+J5uvrMRVxyP/nJT4p3vOMdpT9Ib/01cODA4vnnn6/0Mru8cePGbXeGLf3ans2bNxef/vSnWz32nHPOKTZv3lzGr65r2pF5R0Rx0EEHbfdc5t4+Bx10ULtmXVdXV/zsZz/b7nnMu2O1N+LMvf3a+2f94x//eLFq1artnqerz1zEvQ28+OKLxQUXXFAMHDiw2GOPPYq99tqrGDp0aDF58uRizZo1lV5eCh0VcU1mzpxZ1NfXF/379y969uxZ9O/fv6ivry8efPDBMnw1OXRkxDUx99YtWbKkuPHGG4t//ud/Lt7znvcUBxxwQFFdXV3U1tYWhx12WPHxj3+8uOWWW9r9/w3z7hjtjbgm5t62n//858W3v/3t4p/+6Z+KgQMHFv369Suqq6uLvfbaq3j3u99dfOYzn2nxaUjb01VnXlUU23l9LQAAXZZXpwIAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJPT/A2SCLGtb8KdDAAAAAElFTkSuQmCC", "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": "2023-12-10T18:30:13.243391Z", "iopub.status.busy": "2023-12-10T18:30:13.243080Z", "iopub.status.idle": "2023-12-10T18:30:13.255189Z", "shell.execute_reply": "2023-12-10T18:30:13.254386Z", "shell.execute_reply.started": "2023-12-10T18:30:13.243363Z" } }, "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.00000 0.00000, 10.00000 10.00000,...
11POLYGON ((20.00000 0.00000, 20.00000 10.00000,...
21POLYGON ((30.00000 0.00000, 30.00000 10.00000,...
31POLYGON ((40.00000 0.00000, 40.00000 10.00000,...
41POLYGON ((50.00000 0.00000, 50.00000 10.00000,...
\n", "
" ], "text/plain": [ " var geometry\n", "0 1 POLYGON ((10.00000 0.00000, 10.00000 10.00000,...\n", "1 1 POLYGON ((20.00000 0.00000, 20.00000 10.00000,...\n", "2 1 POLYGON ((30.00000 0.00000, 30.00000 10.00000,...\n", "3 1 POLYGON ((40.00000 0.00000, 40.00000 10.00000,...\n", "4 1 POLYGON ((50.00000 0.00000, 50.00000 10.00000,..." ] }, "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": "2023-12-10T18:30:13.259043Z", "iopub.status.busy": "2023-12-10T18:30:13.258769Z", "iopub.status.idle": "2023-12-10T18:30:13.270721Z", "shell.execute_reply": "2023-12-10T18:30:13.269893Z", "shell.execute_reply.started": "2023-12-10T18:30:13.259021Z" } }, "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": "2023-12-10T18:30:13.272547Z", "iopub.status.busy": "2023-12-10T18:30:13.271994Z", "iopub.status.idle": "2023-12-10T18:30:13.278387Z", "shell.execute_reply": "2023-12-10T18:30:13.277569Z", "shell.execute_reply.started": "2023-12-10T18:30:13.272518Z" } }, "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": "2023-12-10T18:30:13.280057Z", "iopub.status.busy": "2023-12-10T18:30:13.279536Z", "iopub.status.idle": "2023-12-10T18:30:13.284197Z", "shell.execute_reply": "2023-12-10T18:30:13.283175Z", "shell.execute_reply.started": "2023-12-10T18:30:13.280030Z" } }, "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": "2023-12-10T18:30:13.285937Z", "iopub.status.busy": "2023-12-10T18:30:13.285519Z", "iopub.status.idle": "2023-12-10T18:30:21.281332Z", "shell.execute_reply": "2023-12-10T18:30:21.280556Z", "shell.execute_reply.started": "2023-12-10T18:30:13.285911Z" } }, "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": "2023-12-10T18:30:21.282788Z", "iopub.status.busy": "2023-12-10T18:30:21.282097Z", "iopub.status.idle": "2023-12-10T18:30:21.291870Z", "shell.execute_reply": "2023-12-10T18:30:21.291189Z", "shell.execute_reply.started": "2023-12-10T18:30:21.282761Z" } }, "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:\n", " print(\"No feasible solution!\")" ] }, { "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": "2023-12-10T18:30:21.293541Z", "iopub.status.busy": "2023-12-10T18:30:21.292900Z", "iopub.status.idle": "2023-12-10T18:30:21.382198Z", "shell.execute_reply": "2023-12-10T18:30:21.381110Z", "shell.execute_reply.started": "2023-12-10T18:30:21.293515Z" }, "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": "2023-12-10T18:30:21.383798Z", "iopub.status.busy": "2023-12-10T18:30:21.383286Z", "iopub.status.idle": "2023-12-10T18:30:21.468270Z", "shell.execute_reply": "2023-12-10T18:30:21.467748Z", "shell.execute_reply.started": "2023-12-10T18:30:21.383773Z" }, "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": "2023-12-10T18:30:21.469504Z", "iopub.status.busy": "2023-12-10T18:30:21.469238Z", "iopub.status.idle": "2023-12-10T18:30:21.641463Z", "shell.execute_reply": "2023-12-10T18:30:21.640691Z", "shell.execute_reply.started": "2023-12-10T18:30:21.469482Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAUnCAYAAAAy/TFcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAB7CAAAewgFu0HU+AABTjklEQVR4nO3de5xXdb3o//fAcHOGVLz9wmF7R7BdJx4B4cZC22lbNCZJ3bo7hebWriSW0lW04yVBEztuT24eeOs8Ci9tLylKudXACxwiMW0HJYjFEG5vVHIXZv3+oPk2yDAzwMz3O295Ph8PHo/FfNda85n3oLxY6/v9TlVRFEUAAJBKt0ovAACAHSfiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJFRd6QXsLtavXx/PPfdcRETst99+UV1t9ACwu9i0aVO88sorERHx7ne/O3r37r3L51QSZfLcc8/F8OHDK70MAKDC5s+fH8OGDdvl87idCgCQkCtxZbLffvuVtufPnx/vfOc7K7gaAKCcVq5cWboj17wJdoWIK5Pmz4F75zvfGXV1dRVcDQBQKR31vHi3UwEAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIQ6NeJefvnleOCBB2LSpElx4oknxr777htVVVVRVVUVZ511VrvOsX79+rjvvvti/Pjx8f73vz/69esXPXr0iH79+sXRRx8dl156aaxcubLda1q7dm1cffXVMXz48OjXr1/U1tbG4MGD48ILL4w//OEPO/mVAgCUV3VnnvyAAw7YpeOfffbZOOaYY+KNN97Y5rFVq1bFvHnzYt68eXHttdfG9OnT4/TTT2/1fEuXLo2TTjopfvvb32718cWLF8fixYtj+vTp8aMf/ShGjx69S+sGAOhsZbudOmDAgDjhhBN26Ji//OUvpYAbOXJkfOc734mHH344nn766fjpT38an/nMZ6J79+7xxhtvxL/8y7/EQw89tN1zrV69Ok4++eRSwJ177rnxyCOPxFNPPRVXXHFF1NbWxp///Oc47bTT4tlnn935LxQAoAw69UrcpEmTYtiwYTFs2LA44IAD4sUXX4xDDjmk3cd369YtTj/99LjkkkviqKOO2ubxE044IU488cQ45ZRTYvPmzTF+/Ph4/vnno6qqapt9r7nmmli8eHFEREyZMiUuuuii0mNHH310HHfccfHBD34w1q5dGxMmTIhHH310J75iAIDyqCqKoijXJ2secePGjYtbb721Q8576qmnxn/8x39ERMTTTz8dQ4YM2erxN998M/bff//405/+FIMHD45f//rX0a3bthchP/vZz8a///u/R0TEggUL4n3ve1+HrC8ioqGhIQYMGBAREcuXL4+6uroOOzcA0LV1Rge8LV6detxxx5W2ly5dus3jP//5z+NPf/pTRGyJx5YCLiK2erHF3Xff3aFrBADoSG+LiNuwYUNpu6VAe/zxx0vbo0aN2u55hg4dGjU1NRER8cQTT3TgCgEAOtbbIuJmz55d2h40aNA2jy9atKjVx5tUV1fHYYcdts0xAABdTae+sKEcfvWrX8XMmTMjIuJd73pXiy+AWL58eURE1NTUxF577dXq+QYMGBDPPvtsvPLKK7Fhw4bo1atXu9bR0NDQ6uM78l52AABtSR1xGzZsiH/913+NzZs3R0TElVde2eJ+TW9TUltb2+Y5m26nRmx5W5L2RlzTkxUr6a677opJkya1+L56dI5169bF6tWro7a2Nvr06VPp5ewWzLwyzL38zLwy+vbtG5dddlmceuqplV5Km1JH3Be/+MVYsGBBRGx5wcKYMWNa3G/9+vUREdGzZ882z9k82tatW9cBqyyfSZMmld5GhfJ6/fXXK72E3Y6ZV4a5l5+Zl9/FF18s4jrTd77znZg+fXpERLzvfe+LG264Ybv79u7dOyIiNm7c2OZ5m79IYkf+5dN0y3Z7Vq5cGcOHD2/3+XZG6QpcVbfoXrN3p34utti8+rUtG2ZeNmZeGeZefmZefpvXrIooGtPc0UoZcf/+7/8e3/jGNyIi4sgjj4yHHnpoq9ugb9W3b9+I2HJ7tC1r1qwpbbfn9muTrvS+b91r9o66L9xW6WXsFn4/ZUxE0WjmZWTmlWHu5Wfm5ddww7i/xXMC6V6dOmPGjPj85z8fEREHHXRQ/Od//mfst99+rR7TFFhr1qwpvV/c9jRdUdtvv/3a/Xw4AIBySxVxP/nJT+JTn/pUNDY2xjvf+c545JFH2nUFrPkrVlt7ztimTZtKbxY8ePDgXV8wAEAnSRNxjzzySJx++umxadOm2GeffeLhhx8uvadbW4455pjSdvP3lHurBQsWlG6njhw5ctcWDADQiVJE3FNPPRX19fWxYcOGeMc73hE//elP413vele7jz/22GNjzz33jIiI2267Lbb342Kb/yzXU045ZZfWDADQmbp8xD3zzDNx0kknxZo1a6KmpiYefPDBHf7B9D179owvfelLEbHlJzFcc8012+wzd+7cuOmmmyJiy4/mGjZs2K4vHgCgk3Tqq1OfeOKJWLJkSen3r776aml7yZIlW135itj6B9BHbPlh9h/5yEdKL0a4/PLLY88994xf//rX2/2c+++/f+y///7bfPyiiy6KO+64I373u9/FxIkTY8mSJXHGGWdEnz594rHHHosrr7wyNm3aFH369Inrrrtuh79WAIBy6tSImz59etx2W8svi37yySfjySef3Opjb424xx9/PF5++eXS7y+44II2P+cll1wSl1566TYf79u3b8ycOTNGjx4dzz//fEybNi2mTZu21T7veMc74oc//GG8973vbfPzAABUUpe/ndqRDj/88Fi4cGFMnjw5hg4dGnvttVfsscceceSRR8YFF1wQzz77bJx88smVXiYAQJs69Urcrbfeus0t0x1x1llnbXN1blfV1NTExIkTY+LEiR16XgCActqtrsQBALxdiDgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEqoqiqKo9CJ2Bw0NDTFgwICIiFi+fHnU1dV1+OfYZ5994vXXX4+o6hbda/bu8POzrc2rX9uyYeZlY+aVYe7lZ+bl1zTzfv36xWuvvdah5+6MDqje5TPQZaxevXrLRtH4t//4KQ8zLz8zrwxzLz8zL7vq6hx5lGOVtEttbe2WK3EREVXulJdF0fi3bTMvDzOvDHMvPzMvv7/OvEePHhVeSPuIuLeRPn36RERE99p9ou4Lt1V4NbuH308ZE1E0mnkZmXllmHv5mXn5NdwwLtVVT2kPAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJBQp0bcyy+/HA888EBMmjQpTjzxxNh3332jqqoqqqqq4qyzztrh882aNSvGjh0bdXV10atXr6irq4uxY8fGrFmz2n2OtWvXxtVXXx3Dhw+Pfv36RW1tbQwePDguvPDC+MMf/rDDawIAqITqzjz5AQcc0CHnKYoiPvvZz8a0adO2+viKFSvinnvuiXvuuSfOO++8uPHGG6Oqqmq751m6dGmcdNJJ8dvf/narjy9evDgWL14c06dPjx/96EcxevToDlk3AEBnKdvt1AEDBsQJJ5ywU8d+61vfKgXckCFDYsaMGTF//vyYMWNGDBkyJCIipk2bFhdffPF2z7F69eo4+eSTSwF37rnnxiOPPBJPPfVUXHHFFVFbWxt//vOf47TTTotnn312p9YJAFAunXolbtKkSTFs2LAYNmxYHHDAAfHiiy/GIYccskPnWLJkSUyZMiUiIoYOHRpz5syJPn36RETEsGHDYsyYMTFq1KhYsGBBTJ48Oc4+++w47LDDtjnPNddcE4sXL46IiClTpsRFF11Ueuzoo4+O4447Lj74wQ/G2rVrY8KECfHoo4/u7JcNANDpOvVK3Le//e04+eSTd+m26tSpU2PTpk0REXH99deXAq7JHnvsEddff31ERGzatCmuu+66bc7x5ptvxve+972IiBg8eHB85Stf2Wafo48+Os4555yIiHjsscfil7/85U6vGQCgs3XpV6cWRRH33XdfREQMGjQoRowY0eJ+I0aMiCOPPDIiIu69994oimKrx3/+85/Hn/70p4iIGDduXHTr1vKX3fzFFnffffcurh4AoPN06YhbtmxZrFixIiIiRo0a1eq+TY83NDTEiy++uNVjjz/++Db7tWTo0KFRU1MTERFPPPHEziwZAKAsOvU5cbtq0aJFpe1Bgwa1um/zxxctWrTVc+/ae57q6uo47LDD4tlnn93qmPZoaGho9fGVK1fu0PkAAFrTpSNu+fLlpe26urpW9x0wYECLxzX/fU1NTey1115tnufZZ5+NV155JTZs2BC9evVq11qbf/5KWbduXUREbF6zKhpuGFfh1ewmisaIMPOyMvPKMPfyM/Oy27z6tYj429+nXV2Xjrg33nijtF1bW9vqvk23QSO2vJ1IS+dp6xwtnae9EdcVlL7uorH0B5EyMfPyM/PKMPfyM/Oyq67u0nlU0qVXuX79+tJ2z549W923eWy9taCbztPWOdo6T2veevXvrVauXBnDhw9v9/l2Rm1tbbz++utbflPVpZ/u+Pbx138pR4SZl4uZV4a5l5+Zl99fZ96jR48KL6R9unTE9e7du7S9cePGVvfdsGFDafutb0PSdJ62ztHWeVrT1u3ecmhab/fafaLuC7dVeDW7h99PGRNRNJp5GZl5ZZh7+Zl5+TXcMC7VVc8unfZ9+/Ytbb/1FulbrVmzprT91tumTedp6xxtnQcAoKvo0hHX/OpWW6/+bH47860vMmg6z5o1a0rvF9fWefbbb79Uz4cDAHYvXTrijjrqqNJ204/M2p7mjw8ePHinzrNp06ZYunRpi+cAAOhKunTEHXLIIdG/f/+IiJg9e3ar+86ZMyciIg488MA4+OCDt3rsmGOOKW23dp4FCxaUbqeOHDlyZ5YMAFAWXTriqqqqor6+PiK2XEGbN29ei/vNmzevdIWtvr4+qqqqtnr82GOPjT333DMiIm677bZtfixXk1tvvbW0fcopp+zq8gEAOk2XjriIiAkTJpTer2X8+PHbvO3HunXrYvz48RGx5X1dJkyYsM05evbsGV/60pciYstPb7jmmmu22Wfu3Llx0003RcSWH801bNiwjvwyAAA6VKe+xcgTTzwRS5YsKf3+1VdfLW0vWbJkqytfEVv/APomAwcOjAsvvDCuuuqqWLBgQYwcOTK++tWvxmGHHRZLly6NyZMnx8KFCyMi4qKLLoojjjiixbVcdNFFcccdd8Tvfve7mDhxYixZsiTOOOOM6NOnTzz22GNx5ZVXxqZNm6JPnz5x3XXX7fLXDgDQmTo14qZPnx633dbye9s8+eST8eSTT271sZYiLiLiiiuuiJdffjluvvnmWLhwYZxxxhnb7HPOOefE5Zdfvt219O3bN2bOnBmjR4+O559/PqZNmxbTpk3bap93vOMd8cMf/jDe+973tv6FAQBUWJe/nRoR0a1bt7jpppti5syZUV9fH/3794+ePXtG//79o76+Ph588MGYPn16dOvW+pdz+OGHx8KFC2Py5MkxdOjQ2GuvvWKPPfaII488Mi644IJ49tln4+STTy7TVwUAsPM69Urcrbfeus0t010xevToGD169C6do6amJiZOnBgTJ07soFUBAJRfiitxAABsTcQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkFBVURRFpRexO2hoaIgBAwZERMTy5cujrq6uwz/HPvvsE6+//npEVbfoXrN3h5+fbW1e/dqWDTMvGzOvDHMvPzMvv6aZ9+vXL1577bUOPXdndED1Lp+BLmP16tVbNorGv/3HT3mYefmZeWWYe/mZedlVV+fIoxyrpF1qa2u3XImLiKhyp7wsisa/bZt5eZh5ZZh7+Zl5+f115j169KjwQtpHxL2N9OnTJyIiutfuE3VfuK3Cq9k9/H7KmIii0czLyMwrw9zLz8zLr+GGcamuekp7AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACaWJuI0bN8ZNN90U//RP/xTvfOc7o1evXlFbWxtHHnlkfPrTn4558+a16zyzZs2KsWPHRl1dXfTq1Svq6upi7NixMWvWrE7+CgAAOk51pRfQHsuXL4+TTjopnnvuua0+vnHjxvjd734Xv/vd7+KWW26JCy64IL773e9GVVXVNucoiiI++9nPxrRp07b6+IoVK+Kee+6Je+65J84777y48cYbWzweAKAr6fJX4jZt2rRVwL3nPe+JW2+9NebOnRs/+9nPYtKkSVFTUxMREVOnTo1rrrmmxfN861vfKgXckCFDYsaMGTF//vyYMWNGDBkyJCIipk2bFhdffHEZvioAgF3T5a/E3XfffaWAO/roo+Pxxx+P7t27lx4//vjjY8yYMXH00UfHm2++Gd/5znfiggsuiOrqv31pS5YsiSlTpkRExNChQ2POnDnRp0+fiIgYNmxYjBkzJkaNGhULFiyIyZMnx9lnnx2HHXZYGb9KAIAd0+WvxD355JOl7a9//etbBVyT973vfXHyySdHRMSqVati8eLFWz0+derU2LRpU0REXH/99aWAa7LHHnvE9ddfHxFbrvxdd911HfklAAB0uC4fcRs3bixtH3roodvdr/mVsw0bNpS2i6KI++67LyIiBg0aFCNGjGjx+BEjRsSRRx4ZERH33ntvFEWxS+sGAOhMXT7iBg4cWNp+4YUXtrvf0qVLIyKiqqoqjjjiiNLHly1bFitWrIiIiFGjRrX6uZoeb2hoiBdffHFnlwwA0Om6fMSdeeaZ8Y53vCMiIiZPnhybN2/eZp+FCxfGzJkzIyLijDPOKO0fEbFo0aLS9qBBg1r9XM0fb34cAEBX0+Vf2LDffvvFrbfeGp/4xCfiySefjGHDhsWECRNi4MCBsXr16njyySfju9/9bmzcuDHe+973xrXXXrvV8cuXLy9t19XVtfq5BgwY0OJx7dHQ0NDq4ytXrtyh8wEAtKbLR1xExCmnnBILFiyIa6+9Nm6++eYYN27cVo8fcMAB8e1vfzvOO++80tuNNHnjjTdK27W1ta1+nubHrl69eofW2DwAK2XdunUREbF5zapouGFcG3vTIYrGiDDzsjLzyjD38jPzstu8+rWI+Nvfp11dioh7880340c/+lHcf//9Lb7g4L//+79jxowZMXDgwDjppJO2emz9+vWl7Z49e7b6eXr16lXazvINbK4UnkVj6Q8iZWLm5WfmlWHu5WfmZdf8bcq6si6/yjVr1sTo0aNjzpw50b1795g4cWKcffbZceihh8b69evj//2//xf/63/9r3jiiSfiox/9aEydOjXOP//80vG9e/cubTd/pWtLmr+q9a1vQ9KWtm6/rly5MoYPH75D59xRtbW18frrr2/5TVWXf7rj28Nf/6UcEWZeLmZeGeZefmZefn+deY8ePSq8kPbp8hF3ySWXxJw5cyIi4qabbtrqVmrPnj3j+OOPj+OOOy5OOOGEeOyxx+LLX/5yHHfccfGe97wnIiL69u1b2r+tW6Rr1qwpbbd16/Wt2nq+XTk0hWf32n2i7gu3VXg1u4ffTxkTUTSaeRmZeWWYe/mZefk13DAu1VXPLp32RVHELbfcEhFb3mrkrc+Fa1JdXR2XXXZZREQ0NjaWjonYOq7aevFB86tpXeE5bgAA29OlI+6///u/S7cHm36+6fa8733vK203/4kNRx11VIsfb0nzxwcPHrxDawUAKKcuHXHNn1jY9GOztufNN99s8bhDDjkk+vfvHxERs2fPbvUcTbdtDzzwwDj44IN3dLkAAGXTpSOuX79+pTfunTt3bqsh1zzQDjnkkNJ2VVVV1NfXR8SWK23z5s1r8fh58+aVrsTV19dHVVXVLq8fAKCzdOmI69atW+ktQ/74xz/GFVdc0eJ+q1atiq9+9aul35988slbPT5hwoTS1bnx48dv8/Yh69ati/Hjx0fElqt4EyZM6KgvAQCgU3TpiIuImDRpUuyxxx4REXHppZfGmDFj4j/+4z9i4cKFMXfu3Jg6dWq8973vjd/85jcREfGP//iPccIJJ2x1joEDB8aFF14YERELFiyIkSNHxh133BELFiyIO+64I0aOHBkLFiyIiIiLLrpoq5+9CgDQFXX5txgZNGhQ3HfffXHmmWfGq6++Gvfff3/cf//9Le77oQ99KO66664WH7viiivi5ZdfjptvvjkWLlwYZ5xxxjb7nHPOOXH55Zd36PoBADpDl78SFxHx4Q9/OBYvXhyTJ0+OY489Nvbbb7/o0aNH9OnTJw455JA4/fTT4957743//M//jL333rvFc3Tr1i1uuummmDlzZtTX10f//v2jZ8+e0b9//6ivr48HH3wwpk+fHt26pRgJALCb6/JX4prss88+MXHixJg4ceIunWf06NExevToDloVAEBluOwEAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgoaqiKIpKL2J30NDQEAMGDIiIiOXLl0ddXV2Hf4599tknXn/99YiqbtG9Zu8OPz/b2rz6tS0bZl42Zl4Z5l5+Zl5+TTPv169fvPbaax167s7ogOpdPgNdxurVq7dsFI1/+4+f8jDz8jPzyjD38jPzsquuzpFHOVZJu9TW1m65EhcRUeVOeVkUjX/bNvPyMPPKMPfyM/Py++vMe/ToUeGFtI+Iexvp06dPRER0r90n6r5wW4VXs3v4/ZQxEUWjmZeRmVeGuZefmZdfww3jUl31lPYAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmli7hXX301pkyZEiNHjoz/7//7/6JXr17Rv3//eP/73x8XXXRRzJ07t81zzJo1K8aOHRt1dXXRq1evqKuri7Fjx8asWbPK8BUAAOy66kovYEfcdddd8bnPfS5ee+21rT6+cuXKWLlyZcyfPz+ef/75uPfee1s8viiK+OxnPxvTpk3b6uMrVqyIe+65J+65554477zz4sYbb4yqqqrO+jIAAHZZmoj7wQ9+EGeffXY0NjbG/vvvH5/73OfimGOOiX79+sVLL70US5cujfvvvz969Oix3XN861vfKgXckCFDYuLEiXHYYYfF0qVLY8qUKbFw4cKYNm1a7LfffnH55ZeX60sDANhhKSJu0aJFcd5550VjY2N84AMfiPvvvz/23HPPbfYbP358bNy4scVzLFmyJKZMmRIREUOHDo05c+ZEnz59IiJi2LBhMWbMmBg1alQsWLAgJk+eHGeffXYcdthhnfdFAQDsghTPiRs/fnxs2LAh9t1337j77rtbDLgmPXv2bPHjU6dOjU2bNkVExPXXX18KuCZ77LFHXH/99RERsWnTprjuuus6ZvEAAJ2gy0fc4sWL45FHHomIiC9+8Yux77777vA5iqKI++67LyIiBg0aFCNGjGhxvxEjRsSRRx4ZERH33ntvFEWxk6sGAOhcXT7i7rrrrtL2aaedVtpetWpVPP/889u8yKEly5YtixUrVkRExKhRo1rdt+nxhoaGePHFF3dixQAAna/LR9y8efMiImLPPfeMwYMHxw9/+MP4H//jf0S/fv1i4MCBse+++8ahhx4a3/72t2P16tUtnmPRokWl7UGDBrX6+Zo/3vw4AICupMu/sOE3v/lNREQcfPDBMX78+Ljhhhu22WfZsmVx6aWXxo9//OP46U9/Gv3799/q8eXLl5e26+rqWv18AwYMaPG4tjQ0NLT6+MqVK9t9LgCAtnT5iHv99dcjYstz4371q1/FXnvtFVdddVWMHTs23vGOd8Rzzz0XkyZNioceeih+/etfx2mnnRaPP/54dOv2t4uMb7zxRmm7tra21c9XU1NT2t7elb2WNI+/Slm3bl1ERGxesyoabhhX4dXsJorGiDDzsjLzyjD38jPzstu8estTtJr+Pu3qunzErVmzJiIiNmzYEN27d4+HHnpoqxcmDB06NB544IE4+eST46GHHoqnnnoq7r777jj11FNL+6xfv760vb1Xrzbp1atXaTvLN7FJKTqLxtIfRMrEzMvPzCvD3MvPzMuuurrL51FEJIi43r17l0LutNNOa/GVpd26dYurr746HnrooYiImDFjxlYR17t379L29t5HrsmGDRtK2299G5LWtHXrdeXKlTF8+PB2n29n1NbWlq5cRlWXf7rj28Nf/6UcEWZeLmZeGeZefmZefn+deWs/OKAr6fIR17dv31LEnXjiidvd713velcceOCBsWLFivjFL36xzTmatHWLtOlzRbR967W5tp5rVw5N0dm9dp+o+8JtFV7N7uH3U8ZEFI1mXkZmXhnmXn5mXn4NN4xLddWzy6d98+eatfdFCS+//PJWH29+XFsvQGh+Ra0rPM8NAKAlXT7i3vWud5W2N2/e3Oq+TY+/9V72UUcdVdpevHhxq+do/vjgwYPbvU4AgHLq8hH3wQ9+sLS9dOnSVvd94YUXIiLiwAMP3OrjhxxySOltR2bPnt3qOebMmVM6x8EHH7yjywUAKIsuH3FjxowpPcHw7rvv3u5+s2fPLv30hg984ANbPVZVVRX19fURseVKW9MbCL/VvHnzSlfi6uvro6qqapfXDwDQGbp8xO2zzz7xr//6rxER8fDDD8ftt9++zT5vvPFGTJgwofT7z3zmM9vsM2HChNJt1vHjx2/z9iHr1q2L8ePHR8SW27HNzwcA0NV0+YiLiPj2t78df/d3fxcREZ/85Cdj/Pjx8dhjj8Uvf/nLuPXWW2P48OHxzDPPRETE5z73uRg2bNg25xg4cGBceOGFERGxYMGCGDlyZNxxxx2xYMGCuOOOO2LkyJGxYMGCiIi46KKL4ogjjijPFwcAsBO6/FuMRETst99+MWvWrBgzZkwsWbIk/u3f/i3+7d/+bZv9Pv3pT8f3vve97Z7niiuuiJdffjluvvnmWLhwYZxxxhnb7HPOOefE5Zdf3qHrBwDoaCmuxEVseaXoM888E1dffXW8//3vj379+kXPnj2jrq4u/vmf/zkeffTRuOmmm1p9g75u3brFTTfdFDNnzoz6+vro379/9OzZM/r37x/19fXx4IMPxvTp07f6kV0AAF1RiitxTWpqauLCCy8s3RbdWaNHj47Ro0d30KoAAMrPJScAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJVRVFUVR6EbuDhoaGGDBgQERELF++POrq6jr8c+yzzz7x+uuvR1R1i+41e3f4+dnW5tWvbdkw87Ix88ow9/Iz8/Jrmnm/fv3itdde69Bzd0YHVO/yGegyVq9evWWjaPzbf/yUh5mXn5lXhrmXn5mXXXV1jjzKsUrapba2dsuVuIiIKnfKy6Jo/Nu2mZeHmVeGuZefmZffX2feo0ePCi+kfUTc20ifPn0iIqJ77T5R94XbKrya3cPvp4yJKBrNvIzMvDLMvfzMvPwabhiX6qqntAcASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASCh1xE2cODGqqqpKv37+85+3ecysWbNi7NixUVdXF7169Yq6uroYO3ZszJo1q/MXDADQQdJG3K9+9auYOnVqu/cviiI+85nPxIknnhj33HNPrFixIjZu3BgrVqyIe+65J0488cT4zGc+E0VRdOKqAQA6RsqIa2xsjHPPPTc2bdoU+++/f7uO+da3vhXTpk2LiIghQ4bEjBkzYv78+TFjxowYMmRIRERMmzYtLr744k5bNwBAR0kZcf/7f//v+MUvfhGDBg2Kc845p839lyxZElOmTImIiKFDh8aTTz4ZZ5xxRgwbNizOOOOMeOKJJ2Lo0KERETF58uRYunRpp64fAGBXpYu45cuXl66Wff/734+ePXu2eczUqVNj06ZNERFx/fXXR58+fbZ6fI899ojrr78+IiI2bdoU1113XccuGgCgg6WLuM9//vOxevXqGDduXBx77LFt7l8URdx3330RETFo0KAYMWJEi/uNGDEijjzyyIiIuPfeez03DgDo0lJF3J133hkPPPBA9OvXL66++up2HbNs2bJYsWJFRESMGjWq1X2bHm9oaIgXX3xxl9YKANCZqiu9gPb605/+FOeff35EbHne2n777deu4xYtWlTaHjRoUKv7Nn980aJFccghh7R7fQ0NDa0+vnLlynafCwCgLWkibuLEifHSSy/FP/zDP7TrxQxNli9fXtquq6trdd8BAwa0eFx7ND+2UtatWxcREZvXrIqGG8ZVeDW7iaIxIsy8rMy8Msy9/My87Davfi0i/vb3aVeXIuKeeOKJmD59elRXV8eNN94YVVVV7T72jTfeKG3X1ta2um9NTU1pe/Xq1Tu+0AorrbloLP1BpEzMvPzMvDLMvfzMvOyqq1PkUdePuI0bN8Z5550XRVHEBRdcEO9+97t36Pj169eXttt6JWuvXr1K2zta4W1duVu5cmUMHz58h865o2pra+P111/f8puqVE93zOuv/1KOCDMvFzOvDHMvPzMvv7/OvEePHhVeSPt0+Yi78sorY9GiRfF3f/d3cckll+zw8b179y5tb9y4sdV9N2zYUNp+69uQtKWtW7Xl0LTm7rX7RN0XbqvwanYPv58yJqJoNPMyMvPKMPfyM/Pya7hhXKqrnl067RcvXhzf+c53ImLL+7s1v93ZXn379i1tt3WLdM2aNaXttm69AgBUUpe+Ejd16tTYuHFjHHroobF27dq4/fbbt9nn17/+dWn70UcfjZdeeikiIj760Y9GTU3NVlfI2noFafNbol3hhQoAANvTpSOu6fbmCy+8EGeeeWab+1922WWl7WXLlkVNTU0cddRRpY8tXry41eObPz548OAdXS4AQNl06dupHeGQQw6J/v37R0TE7NmzW913zpw5ERFx4IEHxsEHH9zZSwMA2GldOuJuvfXWKIqi1V/NX+zw2GOPlT7eFGFVVVVRX18fEVuutM2bN6/FzzVv3rzSlbj6+vodehsTAIBy69IR11EmTJhQes+X8ePHb/P2IevWrYvx48dHxJb3hpkwYUK5lwgAsEN2i4gbOHBgXHjhhRERsWDBghg5cmTccccdsWDBgrjjjjti5MiRsWDBgoiIuOiii+KII46o5HIBANrUpV/Y0JGuuOKKePnll+Pmm2+OhQsXxhlnnLHNPuecc05cfvnlFVgdAMCO2S2uxEVEdOvWLW666aaYOXNm1NfXR//+/aNnz57Rv3//qK+vjwcffDCmT58e3brtNiMBABJLfyXu0ksvjUsvvbTd+48ePTpGjx7deQsCACgDl50AABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkJOIAABIScQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAkVFUURVHpRewOGhoaYsCAARERsXz58qirq+vwz7HPPvvE66+/HlHVLbrX7N3h52dbm1e/tmXDzMvGzCvD3MvPzMuvaeb9+vWL1157rUPP3RkdUL3LZ6DLWL169ZaNovFv//FTHmZefmZeGeZefmZedtXVOfIoxyppl9ra2i1X4iIiqtwpL4ui8W/bZl4eZl4Z5l5+Zl5+f515jx49KryQ9hFxbyN9+vSJiIjutftE3Rduq/Bqdg+/nzImomg08zIy88ow9/Iz8/JruGFcqque0h4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCKSLu6aefjiuvvDJOPPHEGDBgQPTq1Stqa2tj4MCBcdZZZ8Xjjz++Q+ebNWtWjB07Nurq6qJXr15RV1cXY8eOjVmzZnXSVwAA0LGqK72AtowaNSrmzJmzzcc3btwYzz//fDz//PNx2223xSc/+cmYPn169OzZc7vnKooiPvvZz8a0adO2+viKFSvinnvuiXvuuSfOO++8uPHGG6OqqqrDvxYAgI7S5a/ErVixIiIi+vfvH+eff378+Mc/jvnz58fcuXPj2muvjQMPPDAiIv7v//2/cdZZZ7V6rm9961ulgBsyZEjMmDEj5s+fHzNmzIghQ4ZERMS0adPi4osv7rwvCACgA3T5K3GDBg2KK6+8Mj7+8Y9H9+7dt3psxIgR8clPfjJGjhwZv/vd72LGjBnxuc99Lj7wgQ9sc54lS5bElClTIiJi6NChMWfOnOjTp09ERAwbNizGjBkTo0aNigULFsTkyZPj7LPPjsMOO6zzv0AAgJ3Q5a/EPfDAA3H66advE3BN9t133/jud79b+v2Pf/zjFvebOnVqbNq0KSIirr/++lLANdljjz3i+uuvj4iITZs2xXXXXdcBqwcA6BxdPuLa49hjjy1tL126dJvHi6KI++67LyK2XNkbMWJEi+cZMWJEHHnkkRERce+990ZRFB2/WACADvC2iLiNGzeWtrt12/ZLWrZsWem5daNGjWr1XE2PNzQ0xIsvvthxiwQA6EBvi4ibPXt2aXvQoEHbPL5o0aJWH2+u+ePNjwMA6Eq6/Asb2tLY2BhXXXVV6fenn376NvssX768tF1XV9fq+QYMGNDicW1paGho9fGVK1e2+1wAAG1JH3FTp06N+fPnR0TEKaecEkOHDt1mnzfeeKO0XVtb2+r5ampqSturV69u9zqax1+lbV6zKhpuGFfpZeweisaIMPOyMvPKMPfyM/Oy27z6tUovYYekjrjZs2fH1772tYiI2H///eP73/9+i/utX7++tN3amwFHRPTq1au0vW7dug5YZfn07dt3y0bRmO4PYnpmXn5mXhnmXn5mXnZf/vKXK72Edkkbcf/1X/8Vp5xySmzatCl69eoVd955ZxxwwAEt7tu7d+/SdvMXQbRkw4YNpe23vg1Ja9q69bpy5coYPnx4u8+3My677LK4+OKLt7rySOdat25dVFdXR48ePSq9lN2GmVeGuZefmVfGl7/8ZRHXmZYtWxYnnHBCrFq1Krp37x4zZsxo9VWnpStU0fYt0jVr1pS227r12lxbz7Urh1NPPTVOPfXUSi8DACiDdK9O/eMf/xgf/vCH449//GNUVVXFzTffHKecckqrxzQPrLZegND8ilpXep4bAEBzqSLu1VdfjeOPPz5eeOGFiNjykxc+9alPtXncUUcdVdpevHhxq/s2f3zw4ME7uVIAgM6VJuL+/Oc/x0c+8pH4zW9+ExERV111VXzhC19o17GHHHJI9O/fPyK2fk+5lsyZMyciIg488MA4+OCDd37BAACdKEXErV27Nk466aR4+umnIyLim9/8Znz1q19t9/FVVVVRX18fEVuutM2bN6/F/ebNm1e6EldfXx9VVVW7uHIAgM7R5SNu48aNccopp8STTz4ZERHnn39+XH755Tt8ngkTJkR19ZbXcYwfP36btw9Zt25djB8/PiIiqqurY8KECbu2cACATtTlX5165plnxs9+9rOIiPjQhz4U55xzTvz617/e7v49e/aMgQMHbvPxgQMHxoUXXhhXXXVVLFiwIEaOHBlf/epX47DDDoulS5fG5MmTY+HChRERcdFFF8URRxzROV8QAEAHqCqKoqj0Ilqzo7c0DzrooO3+4PrGxsY499xz4+abb97u8eecc05MmzYtunXr2IuUDQ0NpVe7Ll++vEu8JQkAUB6d0QFd/nZqR+rWrVvcdNNNMXPmzKivr4/+/ftHz549o3///lFfXx8PPvhgTJ8+vcMDDgCgo3X526mdcaFw9OjRMXr06A4/LwBAubjkBACQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AIKHqSi+AjnPXXXfFpEmT4o033qj0UnYb69ati9WrV0dtbW306dOn0svZLZh5ZZh7+Zl5ZfTt2zcuu+yyOPXUUyu9lDZVFUVRVHoRu4OGhoYYMGBAREQsX7486urqOvxzDB48OBYvXtzh5wWA3cmgQYNi0aJFHXrOzugAV+LeRkpX4Kq6RfeavSu7mN3E5tWvbdkw87Ix88ow9/Iz8/LbvGZVRNGY5o6WiHsb6l6zd9R94bZKL2O38PspYyKKRjMvIzOvDHMvPzMvv4Ybxv0tnhPwwgYAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIR2y4j7wx/+EBdeeGEMHjw4ampqol+/fjF8+PC45pprYu3atZVeHgBAm6orvYBymzlzZnziE5+IP//5z6WPrV27Nn7xi1/EL37xi5g+fXo8+OCDceihh1ZwlQAArdutrsT96le/itNPPz3+/Oc/R21tbVxxxRXx1FNPxSOPPBLnnntuRET89re/jZNOOilWr15d4dUCAGzfbnUlbsKECbF27dqorq6On/3sZ3H00UeXHvvQhz4URxxxREycODEWL14c1157bUyaNKmCqwUA2L7d5krcL37xi/j5z38eERHnnHPOVgHX5Ctf+UoMHjw4IiKuu+66ePPNN8u5RACAdtttIu7ee+8tbZ999tkt7tOtW7f41Kc+FRERq1atKkUfAEBXs9tE3OOPPx4RETU1NfG+971vu/uNGjWqtP3EE090+roAAHbGbhNxixYtioiIww8/PKqrt/9UwEGDBm1zDABAV7NbvLBh/fr18eqrr0ZERF1dXav77r333lFTUxNr1qyJ5cuXt/tzNDQ0tPr4ypUr230uAIC27BYR98Ybb5S2a2tr29y/KeJ25G1GBgwYsFNr6wyb16yKhhvGVXoZu4eiMSLMvKzMvDLMvfzMvOw2r36t0kvYIbtFxK1fv7603bNnzzb379WrV0RErFu3rtPW1Bn69u27ZaNoTPcHMT0zLz8zrwxzLz8zL7svf/nLlV5Cu+wWEde7d+/S9saNG9vcf8OGDRER0adPn3Z/jrZuva5cuTKGDx/e7vPtjMsuuywuvvjira480rnWrVsX1dXV0aNHj0ovZbdh5pVh7uVn5pXx5S9/WcR1JaUrVBHtukW6Zs2aiGjfrdcmbT3XrhxOPfXUOPXUUyu9DACgDHaLV6f27t079t1334ho+wUIq1atKkVcV3qeGwBAc7tFxEVE6ScxLFmyJDZt2rTd/RYvXrzNMQAAXc1uE3HHHHNMRGy5VfrLX/5yu/vNnj27tD1y5MhOXxcAwM7YbSLuYx/7WGn7lltuaXGfxsbG+MEPfhAREXvttVccd9xx5VgaAMAO220ibvjw4fGBD3wgIiJuuummmDt37jb7fPe73y39lIbzzz/fK4IAgC5rt3h1apPvfe97MXLkyFi3bl2ccMIJ8Y1vfCOOO+64WLduXdx+++0xbdq0iIgYOHBgfOUrX6nwagEAtm+3irghQ4bEHXfcEf/zf/7P+Mtf/hLf+MY3ttln4MCBMXPmzK3elgQAoKvZbW6nNvnoRz8azz77bFxwwQUxcODA2GOPPWKvvfaKoUOHxuTJk2PhwoVx+OGHV3qZAACtqiqKoqj0InYHDQ0NpfedW758eZd4c2AAoDw6owN2uytxAABvByIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACCh6kovYHexadOm0vbKlSsruBIAoNya/93fvAl2hYgrk1deeaW0PXz48AquBACopFdeeSUOPvjgXT6P26kAAAlVFUVRVHoRu4P169fHc889FxER++23X1RXd+xF0JUrV5au8M2fPz/e+c53duj52ZaZl5+ZV4a5l5+Zl19nz3zTpk2lu3Lvfve7o3fv3rt8TrdTy6R3794xbNiwsnyud77znVFXV1eWz8UWZl5+Zl4Z5l5+Zl5+nTXzjriF2pzbqQAACYk4AICERBwAQEIiDgAgIREHAJCQiAMASEjEAQAk5M1+AQASciUOACAhEQcAkJCIAwBISMQBACQk4gAAEhJxAAAJiTgAgIREHABAQiIOACAhEfc28Ic//CEuvPDCGDx4cNTU1ES/fv1i+PDhcc0118TatWsrvbwUXn755XjggQdi0qRJceKJJ8a+++4bVVVVUVVVFWedddYOn2/WrFkxduzYqKuri169ekVdXV2MHTs2Zs2a1fGLT+zpp5+OK6+8Mk488cQYMGBA9OrVK2pra2PgwIFx1llnxeOPP75D5zP31v3lL3+J22+/Pb7yla/EqFGj4vDDD48999wzevbsGfvvv38ce+yxMWXKlHjttdfadT7z3nUTJ04s/b+mqqoqfv7zn7d5jLm3T/O5tvbr2GOPbfNcXXbmBak98MADxZ577llERIu/jjzyyGLp0qWVXmaXt735RUQxbty4dp+nsbGxOO+881o933nnnVc0NjZ23heTxAc/+MFW59T065Of/GSxYcOGVs9l7u3z8MMPt2vm++67bzFr1qztnse8O8YzzzxTVFdXbzW3xx57bLv7m/uOac+f9YgoRo0atd1zdPWZi7jEnnnmmWKPPfYoIqKora0trrjiiuKpp54qHnnkkeLcc88t/QEbNGhQ8cYbb1R6uV1a8/8gBwwYUJxwwgk7FXHf+MY3SscNGTKkmDFjRjF//vxixowZxZAhQ0qPffOb3+y8LyaJww47rIiIon///sX5559f/PjHPy7mz59fzJ07t7j22muLAw88sDSvM888s9VzmXv7PPzww8WAAQOKT33qU8X3vve94u677y7mzp1bPPnkk8Udd9xRnHbaaUX37t2LiCh69uxZ/OpXv2rxPOa96zZv3lwMGzasiIhi//33b1fEmfuOaZrH5z73ueK5557b7q8XXnhhu+fo6jMXcYkde+yxRUQU1dXVxVNPPbXN41OmTCn9Afv2t79dgRXmMWnSpOL+++8vXnrppaIoimLZsmU7HHHPP/986V/VQ4cOLdauXbvV42vWrCmGDh1a+p4tWbKko7+MVE466aTijjvuKDZt2tTi46+88koxcODA0vdhzpw5Le5n7u23vVk3d88995RmPnbs2G0eN++OMXXq1NI/sr/+9a+3GXHmvuOaZnrJJZfs1PEZZi7ikpo/f37pD+hnPvOZFvfZvHlzMXjw4CIiir333rvYuHFjmVeZ185E3Oc///nSMXPnzm1xn7lz55b2+eIXv9iBK357uv/++0vz+tKXvtTiPube8QYNGlS6rfpW5r3r/vCHPxS1tbWlaLvkkkvajDhz33G7GnEZZu6FDUnde++9pe2zzz67xX26desWn/rUpyIiYtWqVe16wiw7pyiKuO+++yIiYtCgQTFixIgW9xsxYkQceeSREbHle1gURdnWmFHzJxwvXbp0m8fNvXPU1NRERMT69eu3+rh5d4zPf/7zsXr16hg3bly7nlRv7uWXZeYiLqmmV+3V1NTE+973vu3uN2rUqNL2E0880enr2l0tW7YsVqxYERFbz7wlTY83NDTEiy++2NlLS23jxo2l7W7dtv3flbl3vEWLFsUzzzwTEVv+8mrOvHfdnXfeGQ888ED069cvrr766nYdY+7ll2XmIi6pRYsWRUTE4YcfHtXV1dvdr/n/hJuOoeM1n+1b/+J7K9+T9ps9e3Zpu6W5mnvHWLt2bTz//PNx7bXXxnHHHRebN2+OiIjzzz9/q/3Me9f86U9/Ks108uTJsd9++7XrOHPfNXfddVcceeSR0adPn+jbt28cccQRMW7cuHjssce2e0yWmW//b3+6rPXr18err74aERF1dXWt7rv33ntHTU1NrFmzJpYvX16O5e2Wms+2re/JgAEDWjyOrTU2NsZVV11V+v3pp5++zT7mvvNuvfXW7T4VIyLiwgsvjE984hNbfcy8d83EiRPjpZdein/4h3+Ic845p93Hmfuu+c1vfrPV75csWRJLliyJH/zgB/Gxj30sbr311thzzz232ifLzEVcQm+88UZpu7a2ts39myJu9erVnbms3dqOfE+anm8UEb4nrZg6dWrMnz8/IiJOOeWUGDp06Db7mHvHe+973xs33nhjvP/979/mMfPeeU888URMnz49qqur48Ybb4yqqqp2H2vuO2ePPfaIMWPGxD/+4z/GoEGDora2Nl555ZWYPXt23HjjjfHaa6/FvffeG/X19fHwww9Hjx49SsdmmbmIS6j5k4179uzZ5v69evWKiIh169Z12pp2dzvyPWn6fkT4nmzP7Nmz42tf+1pEROy///7x/e9/v8X9zH3nfexjHyuF8bp162Lp0qVx5513xj333BOf+MQn4rrrrouTTz55q2PMe+ds3LgxzjvvvCiKIi644IJ497vfvUPHm/vOWbFiRey1117bfPz444+P8ePHx4knnhgLFy6M2bNnx/e///340pe+VNony8w9Jy6h3r17l7abP/F7ezZs2BAREX369Om0Ne3uduR70vT9iPA9acl//dd/xSmnnBKbNm2KXr16xZ133hkHHHBAi/ua+87ba6+94u///u/j7//+72PYsGFxxhlnxN133x0/+MEP4oUXXoj6+vq49dZbtzrGvHfOlVdeGYsWLYq/+7u/i0suuWSHjzf3ndNSwDU54IAD4sc//nEp0K6//vqtHs8ycxGXUN++fUvb7bl0u2bNmoho361Xds6OfE+avh8RvidvtWzZsjjhhBNi1apV0b1795gxY0arrwwz9473yU9+Mk477bRobGyML37xi7Fq1arSY+a94xYvXhzf+c53ImJLKDS/9dZe5t45Dj300Dj++OMjYsvz5P74xz+WHssyc7dTE+rdu3fsu+++8eqrr0ZDQ0Or+65atar0B6z5ky/pWM2f+NrW96T5E199T/7mj3/8Y3z4wx+OP/7xj1FVVRU333xznHLKKa0eY+6do76+Pu68885Ys2ZNPPTQQ/Ev//IvEWHeO2Pq1KmxcePGOPTQQ2Pt2rVx++23b7PPr3/969L2o48+Gi+99FJERHz0ox+Nmpoac+9ERx11VMycOTMittx+7d+/f0Tk+bMu4pIaPHhwPP7447FkyZLYtGnTdt9mZPHixVsdQ+c46qijStvNZ94S35Ntvfrqq3H88cfHCy+8EBFbrlg0vVF1a8y9czR/64vf//73pW3z3nFNt9peeOGFOPPMM9vc/7LLLittL1u2LGpqasy9E23vzXmzzNzt1KSOOeaYiNhyGfeXv/zldvdr/j5bI0eO7PR17a4OOeSQ0r/gms+8JXPmzImIiAMPPDAOPvjgzl5al/fnP/85PvKRj5TeBuCqq66KL3zhC+061tw7R9ObnEZsfXvIvCvD3DtP87cfaZpxRJ6Zi7ikPvaxj5W2b7nllhb3aWxsjB/84AcRseUJnscdd1w5lrZbqqqqivr6+ojY8q+yefPmtbjfvHnzSv9qq6+v36G3GXg7Wrt2bZx00knx9NNPR0TEN7/5zfjqV7/a7uPNvXPcddddpe3mr6Q07x136623RrHl55Rv91fzFzs89thjpY83BYG5d44XXnghHn744YjY8vy4Aw88sPRYmpmX8we10rE+8IEPFBFRVFdXF0899dQ2j0+ZMmWXfwDw7mrZsmWl2Y0bN65dx/z2t78tqquri4gohg4dWqxdu3arx9euXVsMHTq09D373e9+1wkrz2PDhg3FCSecUJrz+eefv1PnMff2u+WWW4p169a1us+1115b+p4cfPDBxZtvvrnV4+bd8S655JLSzB977LEW9zH3HfOTn/xkmz+7zb300kvFkCFDSnP/7ne/u80+GWYu4hJ7+umniz59+hQRUdTW1hZXXnllMXfu3OLRRx8tzjvvvNIfzoEDBxZ/+ctfKr3cLu3xxx8vbrnlltKvq6++ujS/kSNHbvXYLbfcst3zfO1rXysdN2TIkOL2228vfvGLXxS33377Vv/D+PrXv16+L66LGjt2bGkeH/rQh4pnn322eO6557b767e//e12z2Xu7XPQQQcV/fr1K84999zitttuK5544onimWeeKR5//PHi//yf/1OMHDmyNKuePXsWDz/8cIvnMe+O1Z6IKwpz3xEHHXRQ0b9//2L8+PHFj370o+Kpp54qFi5cWDz88MPFN7/5zWKfffYpzeuYY44p1q9f3+J5uvrMRVxyP/nJT4p3vOMdpT9Ib/01cODA4vnnn6/0Mru8cePGbXeGLf3ans2bNxef/vSnWz32nHPOKTZv3lzGr65r2pF5R0Rx0EEHbfdc5t4+Bx10ULtmXVdXV/zsZz/b7nnMu2O1N+LMvf3a+2f94x//eLFq1artnqerz1zEvQ28+OKLxQUXXFAMHDiw2GOPPYq99tqrGDp0aDF58uRizZo1lV5eCh0VcU1mzpxZ1NfXF/379y969uxZ9O/fv6ivry8efPDBMnw1OXRkxDUx99YtWbKkuPHGG4t//ud/Lt7znvcUBxxwQFFdXV3U1tYWhx12WPHxj3+8uOWWW9r9/w3z7hjtjbgm5t62n//858W3v/3t4p/+6Z+KgQMHFv369Suqq6uLvfbaq3j3u99dfOYzn2nxaUjb01VnXlUU23l9LQAAXZZXpwIAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJCTiAAASEnEAAAmJOACAhEQcAEBCIg4AICERBwCQkIgDAEhIxAEAJPT/A2SCLGtb8KdDAAAAAElFTkSuQmCC", "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": "2023-12-10T18:30:21.643023Z", "iopub.status.busy": "2023-12-10T18:30:21.642562Z", "iopub.status.idle": "2023-12-10T18:30:21.827657Z", "shell.execute_reply": "2023-12-10T18:30:21.827018Z", "shell.execute_reply.started": "2023-12-10T18:30:21.642999Z" }, "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": "2023-12-10T18:30:21.829662Z", "iopub.status.busy": "2023-12-10T18:30:21.829254Z", "iopub.status.idle": "2023-12-10T18:30:21.839362Z", "shell.execute_reply": "2023-12-10T18:30:21.838823Z", "shell.execute_reply.started": "2023-12-10T18:30:21.829632Z" } }, "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 Exception as e:\n", " print(e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "-----------------" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:py312_spopt]", "language": "python", "name": "conda-env-py312_spopt-py" }, "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.0" } }, "nbformat": 4, "nbformat_minor": 4 }