{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "73f22d29",
   "metadata": {},
   "source": [
    "# `kl_divergence_profile` function walkthrough\n",
    "This function calculates the Kullback-Leiber (KL) divergence of different areas and the total population, in order to create a divergence profile across the whole area. This provides a multi-scalar measurement of the segregation, understanding how each area of a population interacts with those areas in close and distant proximity to it.\n",
    "\n",
    "For theoretical background to this methodology, please see the [Olteanu et al. (2019)](https://doi.org/10.1073/pnas.1900192116) article that inspired the creation of this function within PySAL's library.\n",
    "\n",
    "Many thanks to the work of [Cécile de Bézenac](https://github.com/ceciledebezenac/segregation_index), who's primary iteration of a wrapped code to compute this metric was incredibly helpful in directing this function's interpretation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73e3c286",
   "metadata": {},
   "source": [
    "## Contents\n",
    "1. **[The function](#The-function)** - the full function from the PySAL library\n",
    "\n",
    "\n",
    "2. **[Cincinnati example dataset](#Cincinnati-(cincin)-example-dataset)** - an overview of the example dataset used in the walkthrough, and creating the different forms of input acceptable for the data\n",
    "\n",
    "\n",
    "2. **[Making the function work](#Making-the-function-work)** - providing the functions with the correct inputs to ensure that this will run successfully\n",
    "\n",
    "\n",
    "3. **[Using the outputs](#Using-the-outputs)** - how the outputs can be deconstructed to create useful visualisations\n",
    "\n",
    "\n",
    "4. **[How the function works](#How-the-function-works)** - taking a step-by-step work through of the inner workings of the function to understand exactly what it is doing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "666242aa",
   "metadata": {},
   "source": [
    "## The function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "318b6aae",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "from scipy.special import rel_entr as relative_entropy\n",
    "\n",
    "\n",
    "def kl_divergence_profile(populations, coordinates = None, metric = 'euclidean'):\n",
    "    \"\"\"\n",
    "    A segregation metric, using Kullback-Leiber (KL) divergence to quantify the\n",
    "    difference in the population characteristics between (1) an area and (2) the total population.\n",
    "\n",
    "    This function utilises the methodology proposed in\n",
    "    Olteanu et al. (2019): 'Segregation through the multiscalar lens'. Which can be\n",
    "    found here: https://doi.org/10.1073/pnas.1900192116\n",
    "\n",
    "    Arguments\n",
    "    ----------\n",
    "    populations : GeoPandas GeoDataFrame object\n",
    "                  Pandas DataFrame object\n",
    "                  NumPy Array object\n",
    "                  Population information of raw group numbers (not percentages) to be\n",
    "                  included in the analysis.\n",
    "    coordinates : GeoPandas GeoSeries object\n",
    "                  NumPy Array object\n",
    "                  Spatial information relating to the areas to be included in the analysis.\n",
    "    metric : Acceptable inputs to `scipy.spatial.distance.pdist` - including:\n",
    "             ‘braycurtis’, ‘canberra’, ‘chebyshev’, ‘cityblock’, ‘correlation’,\n",
    "             ‘cosine’, ‘dice’, ‘euclidean’, ‘hamming’, ‘jaccard’, ‘jensenshannon’,\n",
    "             ‘kulsinski’, ‘mahalanobis’, ‘matching’, ‘minkowski’, ‘rogerstanimoto’,\n",
    "             ‘russellrao’, ‘seuclidean’, ‘sokalmichener’, ‘sokalsneath’, ‘sqeuclidean’, ‘yule’.\n",
    "             Distance metric for calculating pairwise distances,\n",
    "             using `scipy.spatial.distance.pdist` - 'euclidean' by default.\n",
    "\n",
    "    Returns\n",
    "    ----------\n",
    "    observation : an identifier of the area that forms the centre of the aggregation\n",
    "                  of population, from which the divergence is calculated.\n",
    "    distance : how far the most recently aggregated area is from the 'observation'\n",
    "               area, starting at zero for each observation to represent that only the\n",
    "               'observation' area being aggregated.\n",
    "    divergence : the KL divergence measure, between the aggregated population and the\n",
    "                 total population, will converge to zero for the final row of each\n",
    "                 observation to represent that the total population is covered.\n",
    "    population_covered : the population count within the aggregated population.\n",
    "    Returns a concatenated object of Pandas dataframes. Each dataframe contains a\n",
    "    set of divergence levels between an area and the total population. These areas\n",
    "    become consecutively larger, starting from a single location and aggregating\n",
    "    outward from this location, until the area represents the total population.\n",
    "    Thus, together the divergence levels within a dataframe represent a profile\n",
    "    of divergence from an area. The concatenated object is the collection of these\n",
    "    divergence profiles for every areas within the total population.\n",
    "\n",
    "    Example\n",
    "    ----------\n",
    "    from libpysal.examples import get_path\n",
    "    from libpysal.examples import load_example\n",
    "    cincin = load_example('Cincinnati')\n",
    "    cincin.get_file_list()\n",
    "    cincin_df = gpd.read_file(cincin.get_path('cincinnati.shp'))\n",
    "    cincin_ethnicity = cincin_df[[\"WHITE\", \"BLACK\", \"AMINDIAN\", \"ASIAN\", \"HAWAIIAN\", \"OTHER_RACE\", \"geometry\"]]\n",
    "    cincin_ethnicity.head()\n",
    "    kl_divergence_profile(cincin_ethnicity)\n",
    "    \"\"\"\n",
    "    # Store the observation index to return with the results\n",
    "    if hasattr(populations, 'index'):\n",
    "        indices = populations.index\n",
    "    else:\n",
    "        indices = np.arange(len(populations))\n",
    "\n",
    "    # Check for geometry present in populations argument\n",
    "    if hasattr(populations, 'geometry'):\n",
    "        if coordinates is None:\n",
    "            coordinates = populations.geometry\n",
    "        populations = populations.drop(populations.geometry.name, axis = 1).values\n",
    "    populations = np.asarray(populations)\n",
    "\n",
    "    #  Creating consistent coordinates - GeoSeries input\n",
    "    if hasattr(coordinates,'geometry'):\n",
    "        centroids = coordinates.geometry.centroid\n",
    "        coordinates = np.column_stack((centroids.x, centroids.y))\n",
    "    #  Creating consistent coordinates - Array input\n",
    "    else:\n",
    "        assert len(coordinates) == len(populations), \"Length of coordinates input needs to be of the same length as populations input\"\n",
    "\n",
    "    # Creating distance matrix using defined metric (default euclidean distance)\n",
    "    dist_matrix = squareform(pdist(coordinates, metric = metric))\n",
    "\n",
    "    # Preparing list for results\n",
    "    results = []\n",
    "\n",
    "    # Loop to calculate KL divergence profile\n",
    "    for (i, distances) in enumerate(dist_matrix):\n",
    "\n",
    "        # Creating the q and r objects\n",
    "        sorted_indices = np.argsort(distances)\n",
    "        cumul_pop_by_group = np.cumsum(populations[sorted_indices], axis = 0)\n",
    "        obs_cumul_pop = np.sum(cumul_pop_by_group, axis = 1)[:, np.newaxis]\n",
    "        q_cumul_proportions = cumul_pop_by_group / obs_cumul_pop\n",
    "        total_pop_by_group = np.sum(populations, axis = 0, keepdims = True)\n",
    "        total_pop = np.sum(populations)\n",
    "        r_total_proportions = total_pop_by_group / total_pop\n",
    "\n",
    "        # Input q and r objects into relative entropy (KL divergence) function\n",
    "        kl_divergence = relative_entropy(q_cumul_proportions,\n",
    "                                         r_total_proportions).sum(axis = 1)\n",
    "\n",
    "        # Creating an output dataframe\n",
    "        output = pd.DataFrame().from_dict(dict(\n",
    "            observation = indices[i],\n",
    "            distance = distances[sorted_indices],\n",
    "            divergence = kl_divergence,\n",
    "            population_covered = obs_cumul_pop.sum(axis=1)\n",
    "        ))\n",
    "\n",
    "        # Append (bring together) all outputs into results list\n",
    "        results.append(output)\n",
    "\n",
    "    return(pd.concat(results))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c64ff19",
   "metadata": {},
   "source": [
    "## Cincinnati (cincin) example dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fa2a609",
   "metadata": {},
   "source": [
    "For this example, we will be using an example dataset that is available via the [`examples`](https://pysal.org/notebooks/lib/libpysal/Example_Datasets.html) package of **`libpysal`** - it is a remote dataset, originating from the GeoDa program, an output of the [Center for Spatial Data Science at the Unversity of Chicago](https://spatial.uchicago.edu).\n",
    "\n",
    "This dataset is entitled '[**2008 Cincinnati Crime + Socio-Demographics**](https://geodacenter.github.io/data-and-lab/walnut_hills/)', and contains \"*Crime and socio-demographic data for the Clifton, Walnut Hills, Evanston, and Avondale neighborhoods in Cincinnati, OH for the last 6 months of 2008*\".\n",
    "\n",
    "Throughout this workbook, we will only be inputting the socio-demographic information it possesses. The total population of Cincinnati can be broken down within the data by six different ethnic groups: *AMINDIAN* (American-Indian), *ASIAN*, *BLACK*, *HAWAIIAN*, *WHITE*, and *OTHER_RACE*.\n",
    "\n",
    "This data will be referred to it as ***cincin*** throughout."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "834ae96c",
   "metadata": {},
   "source": [
    "You can read in this example dataset, and then repeat all the code in this notebook, by doing the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2d43890c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example not available: Cincinnati\n",
      "Example not downloaded: Chicago parcels\n",
      "Example not downloaded: Chile Migration\n",
      "Example not downloaded: Spirals\n"
     ]
    }
   ],
   "source": [
    "import geopandas as gpd\n",
    "from libpysal.examples import get_path\n",
    "from libpysal.examples import load_example\n",
    "\n",
    "cincin = load_example('Cincinnati')\n",
    "cincin = gpd.read_file(cincin.get_path('cincinnati.shp'))\n",
    "cincin = cincin.set_index('ID')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad24f596",
   "metadata": {},
   "source": [
    "Having read in this shapefile - as ***cincin*** - the data required for this analysis can be collected using either a *geodataframe* (gdf) or a combination of a *dataframe* (df) and *geoseries* data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2295eac4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>WHITE</th>\n",
       "      <th>BLACK</th>\n",
       "      <th>AMINDIAN</th>\n",
       "      <th>ASIAN</th>\n",
       "      <th>HAWAIIAN</th>\n",
       "      <th>OTHER_RACE</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>726907.0</th>\n",
       "      <td>433.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>POLYGON ((1407302.966 415693.734, 1407473.141 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695744.0</th>\n",
       "      <td>16.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1398841.243 416718.444, 1399605.382 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695762.0</th>\n",
       "      <td>15.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>POLYGON ((1398733.468 416975.853, 1398794.240 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695780.0</th>\n",
       "      <td>12.0</td>\n",
       "      <td>103.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>POLYGON ((1399564.078 416046.633, 1399605.382 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695798.0</th>\n",
       "      <td>52.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1398841.243 416718.444, 1398733.468 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703629.0</th>\n",
       "      <td>5.0</td>\n",
       "      <td>102.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>POLYGON ((1401900.256 422880.130, 1401964.217 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703648.0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>85.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1402856.266 422494.690, 1403591.085 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703666.0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1402861.897 422749.576, 1402925.846 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703686.0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>270.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1402709.472 423554.407, 1403020.559 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703708.0</th>\n",
       "      <td>10.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1402421.734 424653.671, 1402817.865 ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>457 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          WHITE  BLACK  AMINDIAN  ASIAN  HAWAIIAN  OTHER_RACE  \\\n",
       "ID                                                              \n",
       "726907.0  433.0   32.0       0.0    5.0       0.0         2.0   \n",
       "695744.0   16.0   66.0       0.0    0.0       0.0         0.0   \n",
       "695762.0   15.0   12.0       0.0    1.0       0.0         1.0   \n",
       "695780.0   12.0  103.0       0.0    0.0       0.0         2.0   \n",
       "695798.0   52.0   39.0       0.0    1.0       0.0         0.0   \n",
       "...         ...    ...       ...    ...       ...         ...   \n",
       "703629.0    5.0  102.0       0.0    0.0       1.0         2.0   \n",
       "703648.0    0.0   85.0       0.0    1.0       0.0         0.0   \n",
       "703666.0    0.0    0.0       0.0    0.0       0.0         0.0   \n",
       "703686.0    0.0  270.0       1.0    0.0       0.0         0.0   \n",
       "703708.0   10.0   43.0       1.0    0.0       0.0         0.0   \n",
       "\n",
       "                                                   geometry  \n",
       "ID                                                           \n",
       "726907.0  POLYGON ((1407302.966 415693.734, 1407473.141 ...  \n",
       "695744.0  POLYGON ((1398841.243 416718.444, 1399605.382 ...  \n",
       "695762.0  POLYGON ((1398733.468 416975.853, 1398794.240 ...  \n",
       "695780.0  POLYGON ((1399564.078 416046.633, 1399605.382 ...  \n",
       "695798.0  POLYGON ((1398841.243 416718.444, 1398733.468 ...  \n",
       "...                                                     ...  \n",
       "703629.0  POLYGON ((1401900.256 422880.130, 1401964.217 ...  \n",
       "703648.0  POLYGON ((1402856.266 422494.690, 1403591.085 ...  \n",
       "703666.0  POLYGON ((1402861.897 422749.576, 1402925.846 ...  \n",
       "703686.0  POLYGON ((1402709.472 423554.407, 1403020.559 ...  \n",
       "703708.0  POLYGON ((1402421.734 424653.671, 1402817.865 ...  \n",
       "\n",
       "[457 rows x 7 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cincin_gdf = cincin[[\"WHITE\", \"BLACK\", \"AMINDIAN\", \"ASIAN\", \"HAWAIIAN\", \"OTHER_RACE\", \"geometry\"]]\n",
    "cincin_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2a75d6e2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>WHITE</th>\n",
       "      <th>BLACK</th>\n",
       "      <th>AMINDIAN</th>\n",
       "      <th>ASIAN</th>\n",
       "      <th>HAWAIIAN</th>\n",
       "      <th>OTHER_RACE</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>726907.0</th>\n",
       "      <td>433.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695744.0</th>\n",
       "      <td>16.0</td>\n",
       "      <td>66.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695762.0</th>\n",
       "      <td>15.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695780.0</th>\n",
       "      <td>12.0</td>\n",
       "      <td>103.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695798.0</th>\n",
       "      <td>52.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703629.0</th>\n",
       "      <td>5.0</td>\n",
       "      <td>102.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703648.0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>85.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703666.0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703686.0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>270.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703708.0</th>\n",
       "      <td>10.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>457 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          WHITE  BLACK  AMINDIAN  ASIAN  HAWAIIAN  OTHER_RACE\n",
       "ID                                                           \n",
       "726907.0  433.0   32.0       0.0    5.0       0.0         2.0\n",
       "695744.0   16.0   66.0       0.0    0.0       0.0         0.0\n",
       "695762.0   15.0   12.0       0.0    1.0       0.0         1.0\n",
       "695780.0   12.0  103.0       0.0    0.0       0.0         2.0\n",
       "695798.0   52.0   39.0       0.0    1.0       0.0         0.0\n",
       "...         ...    ...       ...    ...       ...         ...\n",
       "703629.0    5.0  102.0       0.0    0.0       1.0         2.0\n",
       "703648.0    0.0   85.0       0.0    1.0       0.0         0.0\n",
       "703666.0    0.0    0.0       0.0    0.0       0.0         0.0\n",
       "703686.0    0.0  270.0       1.0    0.0       0.0         0.0\n",
       "703708.0   10.0   43.0       1.0    0.0       0.0         0.0\n",
       "\n",
       "[457 rows x 6 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cincin_df = cincin[[\"WHITE\", \"BLACK\", \"AMINDIAN\", \"ASIAN\", \"HAWAIIAN\", \"OTHER_RACE\"]]\n",
    "cincin_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f8dc7e9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>726907.0</th>\n",
       "      <td>POLYGON ((1407302.966 415693.734, 1407473.141 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695744.0</th>\n",
       "      <td>POLYGON ((1398841.243 416718.444, 1399605.382 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695762.0</th>\n",
       "      <td>POLYGON ((1398733.468 416975.853, 1398794.240 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695780.0</th>\n",
       "      <td>POLYGON ((1399564.078 416046.633, 1399605.382 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>695798.0</th>\n",
       "      <td>POLYGON ((1398841.243 416718.444, 1398733.468 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703629.0</th>\n",
       "      <td>POLYGON ((1401900.256 422880.130, 1401964.217 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703648.0</th>\n",
       "      <td>POLYGON ((1402856.266 422494.690, 1403591.085 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703666.0</th>\n",
       "      <td>POLYGON ((1402861.897 422749.576, 1402925.846 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703686.0</th>\n",
       "      <td>POLYGON ((1402709.472 423554.407, 1403020.559 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>703708.0</th>\n",
       "      <td>POLYGON ((1402421.734 424653.671, 1402817.865 ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>457 rows × 1 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                   geometry\n",
       "ID                                                         \n",
       "726907.0  POLYGON ((1407302.966 415693.734, 1407473.141 ...\n",
       "695744.0  POLYGON ((1398841.243 416718.444, 1399605.382 ...\n",
       "695762.0  POLYGON ((1398733.468 416975.853, 1398794.240 ...\n",
       "695780.0  POLYGON ((1399564.078 416046.633, 1399605.382 ...\n",
       "695798.0  POLYGON ((1398841.243 416718.444, 1398733.468 ...\n",
       "...                                                     ...\n",
       "703629.0  POLYGON ((1401900.256 422880.130, 1401964.217 ...\n",
       "703648.0  POLYGON ((1402856.266 422494.690, 1403591.085 ...\n",
       "703666.0  POLYGON ((1402861.897 422749.576, 1402925.846 ...\n",
       "703686.0  POLYGON ((1402709.472 423554.407, 1403020.559 ...\n",
       "703708.0  POLYGON ((1402421.734 424653.671, 1402817.865 ...\n",
       "\n",
       "[457 rows x 1 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cincin_geoseries = cincin[['geometry']]\n",
    "cincin_geoseries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d800c3c",
   "metadata": {},
   "source": [
    "## Making the function work"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75e34fb6",
   "metadata": {},
   "source": [
    "The function is designed for you to input data in one of three ways - that are outlined below:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9e9cbe5",
   "metadata": {},
   "source": [
    "### 1. `kl_divergence_profile(populations = geodataframe)`\n",
    "A geodataframe can be used as the sole input to the function. However, for this to work as desired, it **must only contain** columns of the raw data for group counts (for ***cincin*** this is the ethnic groups) along with the geometry column. This geometry must be in polygon form or point form (of the centroid of each area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "33036e31",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-1-8d5700952c7a>:100: RuntimeWarning: invalid value encountered in true_divide\n",
      "  q_cumul_proportions = cumul_pop_by_group / obs_cumul_pop\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>observation</th>\n",
       "      <th>distance</th>\n",
       "      <th>divergence</th>\n",
       "      <th>population_covered</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.964293e-01</td>\n",
       "      <td>472.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>674.319177</td>\n",
       "      <td>3.970494e-01</td>\n",
       "      <td>654.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>755.315782</td>\n",
       "      <td>3.966244e-01</td>\n",
       "      <td>728.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>830.744890</td>\n",
       "      <td>3.082997e-01</td>\n",
       "      <td>872.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>1067.478163</td>\n",
       "      <td>3.110239e-01</td>\n",
       "      <td>981.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>452</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13063.729732</td>\n",
       "      <td>1.881053e-06</td>\n",
       "      <td>37805.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>453</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13099.457827</td>\n",
       "      <td>1.377660e-06</td>\n",
       "      <td>37874.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>454</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13163.435867</td>\n",
       "      <td>1.492472e-06</td>\n",
       "      <td>37919.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>455</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13170.741000</td>\n",
       "      <td>1.595123e-07</td>\n",
       "      <td>37963.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>456</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13223.508518</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>38012.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>208849 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     observation      distance    divergence  population_covered\n",
       "0       726907.0      0.000000  8.964293e-01               472.0\n",
       "1       726907.0    674.319177  3.970494e-01               654.0\n",
       "2       726907.0    755.315782  3.966244e-01               728.0\n",
       "3       726907.0    830.744890  3.082997e-01               872.0\n",
       "4       726907.0   1067.478163  3.110239e-01               981.0\n",
       "..           ...           ...           ...                 ...\n",
       "452     703708.0  13063.729732  1.881053e-06             37805.0\n",
       "453     703708.0  13099.457827  1.377660e-06             37874.0\n",
       "454     703708.0  13163.435867  1.492472e-06             37919.0\n",
       "455     703708.0  13170.741000  1.595123e-07             37963.0\n",
       "456     703708.0  13223.508518  0.000000e+00             38012.0\n",
       "\n",
       "[208849 rows x 4 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kl_divergence_profile(cincin_gdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32ca3ac8",
   "metadata": {},
   "source": [
    "Here we have our first look at the output for the KL divergence profile function. This is four columns, containing:\n",
    "- **observation** - an identifier of the area that forms the centre of the aggregation of population, from which the divergence is calculated.\n",
    "- **distance** - how far the most recently aggregated area is from the 'observation' area, starting at zero for each observation to represent that the 'observation' area is the only one being aggregated.\n",
    "- **divergence** - the KL divergence measure, between the aggregated population and the total population, will converge to zero for the final row of each observation to represent that the total population is covered.\n",
    "- **nearby_population** - the population count within the aggregated population."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38b577b4",
   "metadata": {},
   "source": [
    "### 2. `kl_divergence_profile(populations = geodataframe[groups], coordinates = geodataframe.geometry.centroid)`\n",
    "If a geodataframe contains other columns, that are not desired to be used in the analysis - such as the original ***cincin*** dataset - this can be included as the input by subsetting only the columns that are desired.\n",
    "\n",
    "In order for this to work, the coordinates input must then be used to specify the spatial location of each observation. While point or polygon geoseries format can be passed, output from the geopandas `centroid` command will pass smoothly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e7efe4ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-1-8d5700952c7a>:100: RuntimeWarning: invalid value encountered in true_divide\n",
      "  q_cumul_proportions = cumul_pop_by_group / obs_cumul_pop\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>observation</th>\n",
       "      <th>distance</th>\n",
       "      <th>divergence</th>\n",
       "      <th>population_covered</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.964293e-01</td>\n",
       "      <td>472.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>674.319177</td>\n",
       "      <td>3.970494e-01</td>\n",
       "      <td>654.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>755.315782</td>\n",
       "      <td>3.966244e-01</td>\n",
       "      <td>728.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>830.744890</td>\n",
       "      <td>3.082997e-01</td>\n",
       "      <td>872.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>1067.478163</td>\n",
       "      <td>3.110239e-01</td>\n",
       "      <td>981.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>452</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13063.729732</td>\n",
       "      <td>1.881053e-06</td>\n",
       "      <td>37805.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>453</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13099.457827</td>\n",
       "      <td>1.377660e-06</td>\n",
       "      <td>37874.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>454</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13163.435867</td>\n",
       "      <td>1.492472e-06</td>\n",
       "      <td>37919.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>455</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13170.741000</td>\n",
       "      <td>1.595123e-07</td>\n",
       "      <td>37963.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>456</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13223.508518</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>38012.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>208849 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     observation      distance    divergence  population_covered\n",
       "0       726907.0      0.000000  8.964293e-01               472.0\n",
       "1       726907.0    674.319177  3.970494e-01               654.0\n",
       "2       726907.0    755.315782  3.966244e-01               728.0\n",
       "3       726907.0    830.744890  3.082997e-01               872.0\n",
       "4       726907.0   1067.478163  3.110239e-01               981.0\n",
       "..           ...           ...           ...                 ...\n",
       "452     703708.0  13063.729732  1.881053e-06             37805.0\n",
       "453     703708.0  13099.457827  1.377660e-06             37874.0\n",
       "454     703708.0  13163.435867  1.492472e-06             37919.0\n",
       "455     703708.0  13170.741000  1.595123e-07             37963.0\n",
       "456     703708.0  13223.508518  0.000000e+00             38012.0\n",
       "\n",
       "[208849 rows x 4 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kl_divergence_profile(populations = cincin_gdf[[\"WHITE\", \"BLACK\", \"AMINDIAN\", \"ASIAN\", \"HAWAIIAN\", \"OTHER_RACE\"]],\n",
    "                     coordinates = cincin_gdf[[\"geometry\"]].centroid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f6fdae3",
   "metadata": {},
   "source": [
    "The function will also run if the geometry polygon information is inputted into the *coordinates* argument - without being converted into centroid points. If this is done, the `centroid` command is ran as part of the function to distinguish a point against which distance is calculated.\n",
    "\n",
    "i.e: an alternative input is\n",
    "\n",
    "`kl_divergence_profile(populations = cincin_gdf[[\"WHITE\", \"BLACK\", \"AMINDIAN\", \"ASIAN\", \"HAWAIIAN\", \"OTHER_RACE\"]], coordinates = cincin_gdf[[\"geometry\"]])`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae40db45",
   "metadata": {},
   "source": [
    "### 3. `kl_divergence_profile(populations = dataframe, coordinates = geoseries)`\n",
    "If the population data is contained in a separate object to the spatial, geoseries information - as long as the indexing (order) of both datasets is consistent - they can be inputted separately into the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "303449ae",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-1-8d5700952c7a>:100: RuntimeWarning: invalid value encountered in true_divide\n",
      "  q_cumul_proportions = cumul_pop_by_group / obs_cumul_pop\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>observation</th>\n",
       "      <th>distance</th>\n",
       "      <th>divergence</th>\n",
       "      <th>population_covered</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>8.964293e-01</td>\n",
       "      <td>472.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>674.319177</td>\n",
       "      <td>3.970494e-01</td>\n",
       "      <td>654.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>755.315782</td>\n",
       "      <td>3.966244e-01</td>\n",
       "      <td>728.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>830.744890</td>\n",
       "      <td>3.082997e-01</td>\n",
       "      <td>872.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>726907.0</td>\n",
       "      <td>1067.478163</td>\n",
       "      <td>3.110239e-01</td>\n",
       "      <td>981.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>452</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13063.729732</td>\n",
       "      <td>1.881053e-06</td>\n",
       "      <td>37805.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>453</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13099.457827</td>\n",
       "      <td>1.377660e-06</td>\n",
       "      <td>37874.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>454</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13163.435867</td>\n",
       "      <td>1.492472e-06</td>\n",
       "      <td>37919.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>455</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13170.741000</td>\n",
       "      <td>1.595123e-07</td>\n",
       "      <td>37963.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>456</th>\n",
       "      <td>703708.0</td>\n",
       "      <td>13223.508518</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>38012.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>208849 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     observation      distance    divergence  population_covered\n",
       "0       726907.0      0.000000  8.964293e-01               472.0\n",
       "1       726907.0    674.319177  3.970494e-01               654.0\n",
       "2       726907.0    755.315782  3.966244e-01               728.0\n",
       "3       726907.0    830.744890  3.082997e-01               872.0\n",
       "4       726907.0   1067.478163  3.110239e-01               981.0\n",
       "..           ...           ...           ...                 ...\n",
       "452     703708.0  13063.729732  1.881053e-06             37805.0\n",
       "453     703708.0  13099.457827  1.377660e-06             37874.0\n",
       "454     703708.0  13163.435867  1.492472e-06             37919.0\n",
       "455     703708.0  13170.741000  1.595123e-07             37963.0\n",
       "456     703708.0  13223.508518  0.000000e+00             38012.0\n",
       "\n",
       "[208849 rows x 4 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kl_divergence_profile(populations = cincin_df, coordinates = cincin_geoseries)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b619f8",
   "metadata": {},
   "source": [
    "For those of you on-the-ball, you may notice that this input is rather similar to that suggested as an alternative input to step 2 - the data inputted is exactly the same in this case.\n",
    "\n",
    "However, due to the available format of your own data, or purely due to personal preference, you may prefer to input two separate objects to the *populations* and *coordinates* arguments."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a26c3964",
   "metadata": {},
   "source": [
    "## Using the outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1d97366a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-1-8d5700952c7a>:100: RuntimeWarning: invalid value encountered in true_divide\n",
      "  q_cumul_proportions = cumul_pop_by_group / obs_cumul_pop\n"
     ]
    }
   ],
   "source": [
    "cincin_kl_divergence = kl_divergence_profile(cincin_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "59a1ffc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='distance'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEGCAYAAACToKXdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsTUlEQVR4nO3deXhTZdrH8e/dNF2TttCNQoG2UCh7gbIWGRQcEWUVRVRQUUEBZ3ReR52Z12UcdRwdfR0FVEAFEQEBEVRU1HGjyNKyyL6DVtkRKTuF5/0jaQ20QKFJc9Len+viIpykJzeH8uuT5zzn3GKMQSmlVGAK8ncBSimlLp2GuFJKBTANcaWUCmAa4kopFcA0xJVSKoAF+2KncXFxJiUlxRe7VkqpSikvL2+vMSb+Yr/OJyGekpJCbm6uL3atlFKVkohsv5Sv0+kUpZQKYBriSikVwDTElVIqgPlkTlwpFXhOnjxJfn4+x44d83cplVpYWBjJycnY7Xav7E9DXCkFQH5+Pk6nk5SUFETE3+VUSsYY9u3bR35+PqmpqV7Zp06nKKUAOHbsGLGxsRrgPiQixMbGevXTjoa4UqqYBrjvefsYWyrE567cwfqdBf4uQymlAoal5sQfmvk9hacMz9/Qgh7NkvxdjlLKjx5//HEcDgcHDx6kc+fOdOvWzd8lWZKlRuLBQcKxwlMMn7yUFz/fwOnT2rBCqaruiSee8EqAnzp1ygvVWI+lQjw02EbvFjXp16oWL36+kZFTlnLkRKG/y1JKVZCnnnqKhg0b0q1bN9avXw/AbbfdxowZM/j444+54YYbil/71Vdf0bNnTwDmzZtHhw4daNWqFddffz2HDh0CXLcAeeKJJ+jUqRPTp09n7ty5ZGRk0KlTJ/7whz9w7bXXAnD48GGGDBlCmzZtaNmyJbNnzwZgwoQJ9OvXj+7du5Oens6DDz5Y/P6ffPIJrVq1okWLFnTt2vW8+/GlMk2niMg2oAA4BRQaY7J8UUyo3fUz5fnrW5BRw8k/P17H9n1HGDc4i5ox4b54S6VUKf7+wWrW/HzQq/tsXDOKx3o2OefzeXl5TJ06lWXLllFYWEirVq1o3bp18fNXXnklw4YN4/Dhw0RGRjJt2jQGDBjA3r17efLJJ/n888+JjIzkX//6Fy+88AKPPvoo4FqXPX/+fI4dO0Z6ejrffPMNqampDBw4sHjfTz31FFdccQVvvPEGBw4coG3btsWj/+XLl7Ns2TJCQ0Np2LAh9957L2FhYdx1113F+9q/f/959xMZGenVY+npYkbilxtjMn0V4AChwUEcLzyNiDC0cz1evzWL7fuO0GtUDnnbf/HV2yqlLODbb7+lb9++REREEBUVRa9evc54Pjg4mO7du/PBBx9QWFjIRx99RO/evVm4cCFr1qwhOzubzMxMJk6cyPbtv91LasCAAQCsW7eOtLS04vXZniE+b948nnnmGTIzM+nSpQvHjh3jhx9+AKBr165ER0cTFhZG48aN2b59OwsXLqRz587F+6pevfoF9+MrljqxGRps43jh6eI/X5GRyKzhHbnzrVwGjl3I0/2a0b91sh8rVKpqON+I2ZcutPxuwIABjB49murVq9OmTRucTifGGK688kqmTJlS6tcUjYLP1xTeGMPMmTNp2LDhGdsXLVpEaGho8Z9tNhuFhYUYY0qt9Vz78aWyjsQNME9E8kRkaGkvEJGhIpIrIrl79uy5pGJcI/EzTz6kJzp5f3g2retW44HpK3h67lpO6QlPpSqdzp07M2vWLI4ePUpBQQEffPBBidd06dKFpUuXMm7cuOIRdvv27cnJyWHTpk0AHDlyhA0bNpT42oyMDLZs2cK2bdsAmDZtWvFzV111FS+//HJx0C9btuy8tXbo0IGvv/6arVu3AhRPp1zsfryhrCGebYxpBVwNjBCRzme/wBgz1hiTZYzJio+/6PuaA6458eMnT5fYXi0yhLfuaMug9nUZ+80W7py4hIJjJy/pPZRS1tSqVSsGDBhAZmYm1113HZdddlmJ19hsNq699lo+/vjj4pOS8fHxTJgwgYEDB9K8eXPat2/PunXrSnxteHg4Y8aMoXv37nTq1InExESio6MBeOSRRzh58iTNmzenadOmPPLII+etNT4+nrFjx9KvXz9atGhR/APlYvfjDXK+jxilfoHI48AhY8y/z/WarKwscylNIYZMWMKeguN8cG+nc75m0sLtPD5nNalxkYwfnEVKnO9OGChVlaxdu5ZGjRr5uwyfOnToEA6HA2MMI0aMID09nfvvv7/C6yjtWItI3qWcc7zgSFxEIkXEWfQY+D2w6mLfqCxKm04526D2dZl0R1v2HjpOnzE5LNi01xelKKUqoXHjxpGZmUmTJk349ddfGTZsmL9LKreyTKckAvNFZAWwGPjIGPOJL4opWp1yIR3rxTF7RDZxjlAGvbGYSQsvqauRUqqKuf/++1m+fDlr1qxh8uTJRERE+Lukcrvg6hRjzBagRQXU4lqdUsqceGnqxkYya3hH/jh1OY+8v4oNOwt4tGdj7DZLXb+kVEA516oL5T0XO4V9IZZKvFD7hadTPDnD7IwbnMWwzmlMWridwa8v5pfDJ3xYoVKVV1hYGPv27fN6yKjfFN1PPCwszGv7tNQ68RBb2aZTPNmChL/0aESDRCd/eW8lfcbkMH5wFumJTh9VqVTllJycTH5+Ppe6RFiVTVFnH2+xVIi7RuIXF+JFrmudTEpcJMMm5dF3zAJeHtiSyzMSvFyhUpWX3W73WrcZVXGsNZ0SbOPUaUPhqUsL8tZ1qzFnZDZ1YyMYMnEJY7/ZrB8NlVKVmsVC3FXOpY7GAWrGhDP97g70aJrE03PX8T/TV3DsZOW8BaVSSlkyxE+UI8QBIkKCGXVTS+7v1oD3lv7EwHEL2V2gHbyVUpWPtULcbgPKNxIvIiL8sVs6r9zcinU7Cug9KodVP/1a7v0qpZSVWCvEi6dTvDf9cXWzJGbc0wEB+r+6gI++3+G1fSullL9ZLMS9NxL31KRmNLNHdqJxUhQj3lnK/32mrd+UUpWDxULcPRIv41WbFyPeGcqUoe3p3zqZ/3yxkRHvaOs3pVTgs1aI270/nXLG/oNtPNe/OX/r0YhPV++k/yvf8dOBoz55L6WUqgjWCnEfTad4EhHu6pzG67e14cf9R+g9aj552/f77P2UUsqXLBbivh2Je7q8YQKzRnQkMjSYgWMXMSMv3+fvqZRS3matELf7bk68NPUTnMwekU2bVFfrt6c+WqOt35RSAcVaIV4B0ylni4kIYcLtbRncoS7jvt3KnROXcFBbvymlAoTFQrziplM82W1BPNG7KU/1bcq3G/fSb8wCtu09XKE1KKXUpbBUiId44d4p5XFzu7pMuqMdew8dp/dobf2mlLI+S4W4L9eJl1WHerHMGdGJxCh367fvtvmtFqWUuhCLhXjRnLh/7zpYJzaCmfd0pEuDeB6ZvZq/zVrJyUu8Pa5SSvmSpULcbhNE/Ded4skZZmfs4Czu/l09Ji/6gUGvL9LWb0opy7FUiIsIocFB5b4VrbfYgoSHr87g/wa0YOkPB+g9OocNuwr8XZZSShWzVIiDu+O9RUK8SN+WyUwb2p6jJ0/Rb8wCvli7y98lKaUUYMkQv7iO9xWlZR1X67eUuAjufCuXV7/W1m9KKf+zXojbg/y6OuV8kqLDmT6sIz2aJfHMx+v4n3e19ZtSyr8s1e0erDmd4ik8xMaogS3JSHTy/Gcb2LL3MGMHtSYhKszfpSmlqiDrjcQtOp3iSUS4t2s6r97SivU7C+g1KoeV+dr6TSlV8Swa4tYdiXvq3tTV+s0WJFz/2gI+/P5nf5eklKpiLBjiNsvOiZfG1fotm6Y1oxn5zjJemLdeW78ppSqM9ULcbv3plLPFOUKZfFc7rm+dzEv/3cTwydr6TSlVMcoc4iJiE5FlIvKhLwsKpOkUT6HBNp7t35z/vaYR89bs5LpXviP/lyP+LkspVcldzEj8j8BaXxVSxOqrU85HRLjzsjTeuK0N+fuP0Gd0DrnbtPWbUsp3yhTiIpIMXAOM9205rtvRHg/wtdddGiYwa0Q2jtBgBo5byLu5P/q7JKVUJVXWkfiLwIPAOYfIIjJURHJFJHfPnj2XXFCgTqecrX6Cg/dHZNMuNZYHZ3zPkx9q6zellPddMMRF5FpgtzEm73yvM8aMNcZkGWOy4uPjL7mgQJ5OOZur9VsbbuuYwvj5WxkyQVu/KaW8qywj8Wygl4hsA6YCV4jI274qKBBXp5xPsC2Ix3s14em+zcjZtJe+o3PYqq3flFJecsEQN8b8xRiTbIxJAW4E/muMucVXBYUGB3HylKl0a61valeHt+9sx/7DJ+gzOof5G7X1m1Kq/Ky3Ttzd3edEJeyk0z4tljkjO1EjKoxb31zMxAXb9E6ISqlyuagQN8Z8ZYy51lfFgDX6bPpS7eoRzBzekcsbxvPYnNX8ddYqyzTBUEoFHuuNxO1FHe8rz7z42RyhwYwdlMXwLvWYstjV+m2/tn5TSl0C64V4cbPkyj06DQoSHuyewYsDMln24wF6j57P+p3a+k0pdXEsGOKVfyTuqU/LWrw7rAPHT56m35gcPlujrd+UUmVn2RA/VknnxEuTWTuGOSM7US/BwdBJuYz5apOe8FRKlYn1QtxeNaZTzlYjOox3h3Xg2uY1efaT9dw/bbm2flNKXZAF27NVrekUT2F2Gy/dmEnDRAf/nreBrfuOME5bvymlzsN6I/HiEK9aI/EiIsLIK9J59ZbWbNzlav32ff4Bf5ellLIoC4a4ezqlCs2Jl6Z70xrMvKejq/Xbq9/xwQpt/aaUKslyIR5ShadTztYoKYrZI7NpnhzNvVOW8by2flNKncVyIV7Vp1POFucIZfKd7RmQVZuX/7uJeybncfi4tn5TSrlYL8TtGuJnCwkO4pnrmvHotY35bM0urntlgbZ+U0oBVgzx4jlxnU7xJCIM6ZTKm7e35acDR+k9Kocl2vpNqSrPgiHuKqky3sXQG37XIJ73R2QTFW7npnELmbbkB3+XpJTyI8uGeFVfnXI+9eIdvD88m/ZpsTw0cyVPfLCGQv2hp1SVZLkQFxFXs2SdEz+v6Ag7b97WhtuzU3gjZytDJuby61Ft/aZUVWO5EIeiZsk6J34hwbYgHuvZhGf6NeO7zXvpOyaHLXsO+bsspVQFsmiIV55myRXhxrZ1ePuOdhw4cpI+o3P4duMef5eklKogFg3xIJ0Tv0jt0mKZPSKbmjHh3PbmEt7M2ap3QlSqCrBmiFeyjvcVpXb1CGbc05ErMhL4+wdr+Ousldr6TalKzpohrtMpl8wRGsxrt7RmxOX1mLL4R255fRH7Dh33d1lKKR+xaIjr6pTyCAoS/nxVBv+5MZMVPx6g9+gc1u086O+ylFI+YN0Q1ys2y613pqv124nC01w3ZgHzVu/0d0lKKS+zZojbdTrFW1rUjuGDeztRP8HBsLfzGP2ltn5TqjKxZIiH2HQ6xZsSo8KYNqwDPZvX5LlP13Oftn5TqtKwXHs20NUpvhBmt/GfGzNpWMPJc5+uZ9vew4wdnEWitn5TKqBZciSu68R9Q0QYcXl9xg5qzcbdh+g1aj4rfjzg77KUUuVg0RDXOXFf+n2TGrw3vCN2WxA3vPYds5f/5O+SlFKXyKIhHsQJnU7xqYwaUcwekU2L5Bj+OHU5z326Tlu/KRWArBnidj2xWRFiHaG8fWc7bmxTm9FfbmbY23kc0tZvSgWUC4a4iISJyGIRWSEiq0Xk774uqmg6RZfC+V5IcBD/7NeMx3o25ou1u+j/ygJ+3K+t35QKFGUZiR8HrjDGtAAyge4i0t6XRWl3n4olItyencqEotZvo3NYtGWfv8tSSpXBBUPcuBTdpNru/uXTIbJ2vPePzu7WbzHhdm55fRFTF2vrN6Wsrkxz4iJiE5HlwG7gM2PMIl8WFWovapasIV7R6sU7mOVu/fbweyv5+wertfWbUhZWphA3xpwyxmQCyUBbEWl69mtEZKiI5IpI7p495WtK8NtIXFeo+ENR67ch2am8mbON2ycs4dcj2vpNKSu6qNUpxpgDwFdA91KeG2uMyTLGZMXHx5erKJ1O8b9gWxCP9mzMv65rxsIt++g7JofN2vpNKcspy+qUeBGJcT8OB7oB63xZVGiwTqdYxYA2dZh8Z3sOHHW1fvtmg7Z+U8pKyjISTwK+FJHvgSW45sQ/9GVRoXadTrGStqnVmT0im1ox4dz25mLemK+t35SyigveAMsY8z3QsgJqKabTKdZTu3oEM+/pyH3TlvPEh2tYv7OAf/RpSkiwJa8XU6rKsOT/QA1xa4p0t34beXl9puX+yM3jF7JXW78p5VcWDfGiOXGdTrGaoCDhgasa8p8bM/k+/1d6j8ph7Q5t/aaUv1g0xHUkbnVFrd8KT5/mulcW8Km2flPKLywa4u6RuIa4pbWoHcOckZ1IT3AwbJK2flPKH6wZ4u7VKSc0xC2vqPVb70xX67c/TtXWb0pVJGu2Z9MrNgNKmN3GiwM8Wr/tO8zYQVnUiNbWb0r5mjVH4jqdEnBEhOFd6jN2UBab3a3flmvrN6V8zpIhXrT2WK/YDDxXNk5k5vCOhARr6zelKoIlQ9wWJNhtotMpAaqo9VtmbVfrt2c/0dZvSvmKJUMctFlyoIt1hPL2He0Y2LY2Y77S1m9K+YqFQzxIR+IBLiQ4iKf7NuPxno3577rdXDdGW78p5W3WDnGdEw94IsJt2alMuL0NO37V1m9KeZt1Q9yu0ymVyWXp7tZvEXZuHr+IKdr6TSmvsG6I63RKpZPmbv3WsX4cf3lvJY/P0dZvSpWXxUNc/4NXNtHhdt64NYs7OqUyYYG2flOqvCwc4jadE6+kgm1BPHJtY569rjkLt+yjj7Z+U+qSWTbEQ3Q6pdK7oU1t3rmrPQfdrd++1tZvSl00y4a4TqdUDW1SqjN7ZDbJ1SK4/c3FvK6t35S6KNYNcbuGeFWRXC2CGXd34MrGifzjwzU8NPN7/RSmVBlZN8SDbXor2iokMjSYV25uzR+uqM+7ufncMn6Rtn5TqgwsHOI6J17VBAUJf/p9Q14e2LK49duan7X1m1LnY/EQ15F4VdSzRU1m3N2RU6cN/V9dwCertPWbUudi3RC36xLDqqxZcjRzRmaTnujk7rfzePmLjXrCU6lSWDfE3dMp+h+36kqICmPa0Pb0bVmL5z/bwB+mLufoCZ1iU8qTJduzgSvETxsoPG2w28Tf5Sg/CbPbeOGGFjRIdPLsp+vYtvcw4wZr6zelilh4JK4t2pSLiHBPl3qMG5TFlj2H6DlqPst++MXfZSllCdYNcXtRizb9+KxcujVO5L3h2YTZgxgwdiHvL9PWb0pZN8SLO97rSFz9pmENJ7NHdKJl7Rjum7acf2nrN1XFWTjEdTpFla56ZAiT7mjHTe3q8MpXmxk6KVdbv6kqy8IhXjQS1+kUVVJIcBBP9WnKE72b8OX6PfQbk8MP+7T1m6p6LhjiIlJbRL4UkbUislpE/lgRhf02J64jcVU6EWFwhxTeGtKWXQeP03v0fBZq6zdVxZRlJF4I/I8xphHQHhghIo19WxaE2HQ6RZVNdv043h+RTfXIEG4Zv4h3FmnrN1V1XDDEjTE7jDFL3Y8LgLVALV8XVjwS1+kUVQapcZHMGpFNdv04/jprJY/NXqWt31SVcFFz4iKSArQEFpXy3FARyRWR3D17yn9z/+I5cZ1OUWUUFWbnjdvacNdlqUz8bju3vrmYA0dO+LsspXyqzCEuIg5gJnCfMabEreWMMWONMVnGmKz4+PhyF1a0OuWEjqbURbAFCX+7pjHP9W/Okq2/0Gd0Dpt2a+s3VXmVKcRFxI4rwCcbY97zbUkuujpFlcf1WbWZMrQdh44X0nd0Dl+u3+3vkpTyibKsThHgdWCtMeYF35fkoqtTVHm1rlud2SM7kVw9gjsmLGH8t1v0hmqq0inLSDwbGARcISLL3b96+LguvdhHeUWtmHBm3tOBq5rU4MmP1vLnGdr6TVUuF7yLoTFmPlDhtxHU6RTlLREhwYy+qRUvfrGRl77YyNa9h3n1ltbEO0P9XZpS5Wb9KzZ1OkV5QVCQ8KcrGzD6plas/vlXeo+az+qff/V3WUqVm2VDPNgWhC1IdDpFedU1zZOYcXdHDND/le/4ZNUOf5ekVLlYNsRBmyUr32haK5rZI7PJSHJy99tLeUlbv6kAFgAhriNx5X0JzjCm3NWefi1r8cJnGxg5ZZm2flMBybLt2cC1QkXnxJWvhNltPH9DCxrWcPLMJ+vYvs/V+i0pOtzfpSlVZtYeidt1OkX5logw7Hf1GD84i217j9BrVA5LtfWbCiDWDnGdTlEVpGujRN4b3pFwu40bxy7kvaX5/i5JqTKxdIiHaIirCtQg0cnsEdm0qhPDn95dwT8/Xsspbf2mLM7SIR4abNPpFFWhqrlbv93Svg6vfb2Fu97KpeDYSX+XpdQ5WTzEg/TEpqpwdlsQT/Zpxj96N+HrDXvoN2aBtn5TlmX5ENdb0Sp/GdQhhUlD2rK74Di9Rs9nwea9/i5JqRIsHuK6xFD5V8f6ccwekU2cI5TBry/m7YXb/V2SUmewdojrEkNlASlxkbw3vCOXpcfxv++v4pH3V3FSPyEqi7B2iOvqFGURUWF2xt/ahqGd05i0cDu3vqGt35Q1WDzEbRriyjJsQcJfezTi39e3IHfbL/QencOm3QX+LktVcRYP8SCOn9TpFGUt/VsnM2Voew4fP0Xf0Qv4cp22flP+Y+0Qt+t0irKm1nWrMWdkNnViIxgycQljv9msd0JUfmHtEA+2UXjaUKgnkZQF1YwJZ/rdHbi6aQ2enruOB6Zr6zdV8Swe4q7ydK24sqqIkGBGDWzFfd3Smbk0n4FjF7K74Ji/y1JVSECEuK4VV1YWFCTc160BY25uxZodB+kzKodVP2nrN1UxrB3idu14rwJHj2au1m8A17/6HXNXaus35XvWDnHteK8CTNNa0bw/MptGSU6GT17Ki59v4LTeCVH5kKU6+6Q8/FGp23UkrgJJgjOMKUPb89f3VvHi5xvZsKuAf1/fgogQS/13U5VEQHxX/f7/vil+3DuzJvdeUZ968Q5ExI9VKXVuocE2/n19czJqOHn647Vs33eEcYOzqBmjrd+Ud4kv1rZmZWWZ3Nzci/66c43Ey2Jo5zT6t06mfryDoCANd2UdX67bzR+mLCPUbuO1Qa1pXbeav0tSFiQiecaYrIv+usoS4qUJDhKGd6nH1c2SSE9wEGyz9CkAVYlt3FXAnW/lsuPAMZ7u14z+rZP9XZKyGA3xixDnCOXmdnW4snEiDRKdhARruCvf++XwCUa8s5QFm/cxtHMaD3XPwKafGpVblQjxq5ok8tmaXfjqZH/d2Ah6Z9aia0YCDRKdhIfYfPNGqso6eeo0//hwDW99t53LG8bz0sCWOMPs/i5LWUCVCHF/SYuP5JpmSVzuDndHaECcD1YW9vbC7Tw+ZzUpcZGMH5xFSlykv0tSfuazEBeRN4Brgd3GmKZl2WllC/FzSYuLpEezJDo3iKdBooOYiBB/l6QCyILNexk+eSkAY25qRcf6cX6uSPmTL0O8M3AIeEtDvGxS4yLp0awG2fXiSE90EucI0eWQqlTb9x3mzom5bNl7mMd7NmZQhxR/l6T8xKfTKSKSAnxYWUN8eJd6HDx2krcX/uDT96kbG8E1zZJonxZL/QQHSdFhGu6KgmMnuW/qcr5Yt5tb2tfhsZ5NsOtKqirH7yEuIkOBoQB16tRpvX37xTeUtdJIPDEqlFox4Sz94YDP36tubAQ9miWRVbca6QlOkquF61r3KubUacOzn67jta+30CEtljE3t6JapE7PVSV+D3FPgTYSt6ra1cPp0TSJlnWqkZ7ooG71CF3rXsm9tzSfh2eupEZ0GK/fmkV6otPfJakKcqkhrsssLOzH/Ud57ZstJbYnVwunR7MkmidHUz/BQWpcJKHBuhyyMujXKpmUuEiGvpVH3zELeGlgJldkJPq7LGVhGuIBKP+Xo4wtJdxrxYTTvWkNmtaKIj3BSb14h651D0Ct6rhavw2dlMsdE3N5uHsGQzun6fkTVaoLhriITAG6AHEikg88Zox53duF/HrkpLd3WeX8dOAor8/fWmJ7zegwft+kBo1rRlE/wUH9BAdReoGJpdWMCWf6sI48MGMF//x4Het3FvB0v2aE2fWHsjrTBUPcGDOwIgqJjrj0UGmbUp3F2/Z7sZrK5edfjzFhwbYS25Oiw+jWKJGMJCfpCU7SExx6Ms1CwkNsjBrYkoaJTl74bANb9x3mtUGtSXCG+bs0ZSF6xaYqoUZUGJdnJJBRwxXs9RMcxDtD9eO8H328cgd/encFMRF2xg3OommtaH+XpLxMT2wqr9l58BhTFpdcM58YFcrvGsTTINFJ/QQH6YlOaupa9wpxdbMk6sRGcNfEXPq/uoDnr8/kmuZJ/i5LWYCGuCqzXQeP825ufontCc5QOqXHFU/JpCc6SK4WoXfo87ImNaOZPbITd7+dx4h3lrJ+Vzr3dU3XawqqOA1xVW67C47z3tKfSmyPc4TSoV4sDdxTMumJDurGRurViOUQ7wzlnbva8b+zVvHSFxvZuKuA52/Q1m9Vmf7LK5/Ze+g4H6z4ucT2OEcIbVOrUz/eQf1E1+g9NS5SV16UUWiwjWf7N6dhDSdPz13L9leOMO7WLGpp67cqSUNcVbi9h04wd+XOEturR4bQyn11atEJ1XrxDiL11r8liAh3XpZGvQQHf3hnGb1HzXe3fqvu79JUBdPVKcryYiLsZNaOKQ72+gmuE6vR4brWHWDT7gLunJjLz9r6LaDp6hRVaR04cpKv1u/hq/V7ztgeFRZMs+Ro0t2hXhTysY5QP1XqH/UTnLw/IpsR7yzlgekrWL/zIA9f3UhPLFcRGuJVQK2YcHYdPEahr/ra+cnBY4XkbNpHzqZ9Z2x3hAbTuGaUa6WMe+SenuggoRKvdY+JCGHC7W158sM1jPt2K5t2H+I/A1vqlblVgIZ4FfDTgaP+LqFCHTpeyOKt+1m89cyreMPtNholFY3andR3z73XjK4ct/6124L4e++mNKjh5LHZq+k3ZoG2fqsCNMRVlXH05CmW/nCgxD3iQ4KDaOheJeMKdlfQ16kemGvdb25Xl7Q4B8Mn59F7dA6v3Kyt3yozPbGp1DnYbUK9eNeVqfXjHcWrZurGRhISbP217j/sO8Kdby1h8x5t/RYI9MSmUl528pRh3c4C1u0sOGO7LUhIiY1wXaGa6CienkmLt9Za9zqxEcy8pyP3T1vOI7NXs25nAY/30tZvlY2GuFIX6dRpw+Y9h9m85zCfrP5tuwjUqR5RfDK1aMVMvQQHDj+tdXeG2XltUBb/nreeV77azOY9h3jl5tZ6t8pKRENcKS8xBrbvO8L2fUf4fO3uM56rFRNOPXeoF91fpn68s1y3YC4rW5DwUPcMGiQ6eGjmSnqPzmH8rVk00NZvlYKGuFIV4KcDR/npwFG+2XDmWvd4Z6jHUsjflkPGRoZ4fTlk35bJpMRGMnRSHv3GLOA/N2bStZG2fgt0GuJK+dGeguPsKTjOgs1nrnWvFmF3tdjzGLmnJzhJjCrfWveWRa3f3srjzrdyeah7BsO09VtA0xBXyoJ+OXKSxdv2l+hY5QwNPiPYi06q1oop+1r3pOhw3h3WgT/PWMEzH69jg7Z+C2ga4koFkILjhSz/8QDLfzxwxvYwe9BvFzEV3frXvdY9uJTVKOEhNl52t357/rMNbNl7mLGDWpMQpa3fAo2GuFKVwLGTp1n100FW/XTwjO0htiBS4yKLr04tCvmUuAhCg23c2zWd9EQnf3p3Ob1G5TBucBbNkrX1WyDREFeqEjtx6jTrdxWwfldBiefS4iOLg31Qh7q89vUWeo6az6ibWnJt85p+qFZdCr1iUylVqrsuS6VBotN1xaof17pXFXrFplLKq8Z9u7XEttZ1q9EiOcbjpKqDmAi9cMifNMSVUmWWt/0X8rb/UmJ79cgQumYk0Cw5uvgEa5zD+2vdVUka4kqpctt/+ATT8/KZnpdf4rk2KdW4IiORJjWjqJ/gICk6TMPdizTElVI+tWTbLyzZVnL0HiRwQ1ZtujRMoHFSFMnVKsd93SuahrhSyi9OG5i65EemLvmxxHNxjlBuz07h940TSY2LLHWtu3LREFdKWc7eQ8d57tP1PPfp+hLPZdRwMiQ7lV6ZNfUqUzTElVIBZt3OAh6c+T0Pzvy+xHPR4Xb+2iODni1qEhFSNeKtavwtlVJVwq9HT/LQzJU8NHNlqc/3aFaDq5rUoEVyDDWiwyrFSL5MIS4i3YH/ADZgvDHmGZ9WpZRSPjB35U7mrtx5zuevapJI29RY0uIjqREVRo2oMGIi7JZeTXPBEBcRGzAauBLIB5aIyBxjzBpfF6eUUhXp09W7+HT1rvO+Ji0uklZ1q9EiOZqGNaKoGRNGgjPMb31XyzISbwtsMsZsARCRqUBvQENcKVXlbNl7mC17DzOjlDXxANueuaZC6ynLj45agOcaoHz3tjOIyFARyRWR3D179pz9dJlMuqMtaXGRtE2tTrNa0aTFR17SfpRSqqooy0i8tMmgEnfNMsaMBcaC6wZYl1LMZenx/PeBLpfypUopVSWVZSSeD9T2+HMy8LNvylFKKXUxyhLiS4B0EUkVkRDgRmCOb8tSSilVFhecTjHGFIrISOBTXEsM3zDGrPZ5ZUoppS6oTOvEjTFzgbk+rkUppdRF0rvKKKVUANMQV0qpAKYhrpRSAUxDXCmlAphPut2LyB5g+yV+eRyw14vl+Fqg1QuBV7PW63uBVnOg1QsXrrmuMSb+YnfqkxAvDxHJNcZk+buOsgq0eiHwatZ6fS/Qag60esF3Net0ilJKBTANcaWUCmBWDPGx/i7gIgVavRB4NWu9vhdoNQdaveCjmi03J66UUqrsrDgSV0opVUYa4kopFcAsE+Ii0l1E1ovIJhF52I911BaRL0VkrYisFpE/urc/LiI/ichy968eHl/zF3fd60XkKo/trUVkpfu5l8SH3VZFZJv7vZaLSK57W3UR+UxENrp/r2aFmkWkocdxXC4iB0XkPqsdYxF5Q0R2i8gqj21eO6YiEioi09zbF4lIig/qfU5E1onI9yIyS0Ri3NtTROSox7F+1SL1eu17wNv1nqfmaR71bhOR5e7tFXOMjTF+/4XrFrebgTQgBFgBNPZTLUlAK/djJ7ABaAw8DjxQyusbu+sNBVLdfw+b+7nFQAdc3ZE+Bq72Yd3bgLiztj0LPOx+/DDwLyvV7PFvvxOoa7VjDHQGWgGrfHFMgeHAq+7HNwLTfFDv74Fg9+N/edSb4vm6s/bjz3q99j3g7XrPVfNZzz8PPFqRx9gqI/HiZszGmBNAUTPmCmeM2WGMWep+XACspZSeoh56A1ONMceNMVuBTUBbEUkCoowx3xnXv8hbQB/fVl9qbRPdjyd6vL+Vau4KbDbGnO8KX7/Ua4z5BthfSi3eOqae+5oBdC3PJ4nS6jXGzDPGFLr/uBBXZ65z8ne95+H343uhmt37vgGYcr59eLtmq4R4mZoxVzT3R5mWwCL3ppHuj6VveHyMPlfttdyPz97uKwaYJyJ5IjLUvS3RGLMDXD+cgASL1Qyu0YbnN72VjzF495gWf407aH8FYn1WOQzBNeorkioiy0TkaxG5zKMmf9frre+Bij6+lwG7jDEbPbb5/BhbJcTL1Iy5IomIA5gJ3GeMOQi8AtQDMoEduD42wblrr+i/U7YxphVwNTBCRDqf57WWqFlc7f56AdPdm6x+jM/nUmqssPpF5G9AITDZvWkHUMcY0xL4E/COiERdoKaKqNeb3wMV/f0xkDMHJBVyjK0S4pZqxiwidlwBPtkY8x6AMWaXMeaUMeY0MA7XFBCcu/Z8zvzo6tO/kzHmZ/fvu4FZ7vp2uT+6FX2E222lmnH9wFlqjNnlrt3Sx9jNm8e0+GtEJBiIpuzTC2UmIrcC1wI3uz++456W2Od+nIdrjrmBv+v18vdAhRxfj/33A6YVbauoY2yVELdMM2b3/NPrwFpjzAse25M8XtYXKDo7PQe40X1WORVIBxa7P2oXiEh79z4HA7N9VHOkiDiLHuM6mbXKXdut7pfd6vH+fq/Z7YyRi5WPsQdvHlPPffUH/lsUst4iIt2Bh4BexpgjHtvjRcTmfpzmrneLBer15veAz+v10A1YZ4wpniapsGN8sWdnffUL6IFrJchm4G9+rKMTro8v3wPL3b96AJOAle7tc4Akj6/5m7vu9XisjgCycH0TbgZG4b5C1gc1p+E6c78CWF10/HDNpX0BbHT/Xt1CNUcA+4Boj22WOsa4fsDsAE7iGiHd4c1jCoThmkrahGu1QpoP6t2Ea4616Hu5aOXDde7vlRXAUqCnRer12veAt+s9V83u7ROAu896bYUcY73sXimlAphVplOUUkpdAg1xpZQKYBriSikVwDTElVIqgGmIK6VUAAv2dwFKlZWIPA4cAqKAb4wxn5/jdX2ADcaYNRVXnVL+oSNxFXCMMY+eK8Dd+uC6651SlZ6GuLI0EfmbuO4f/TnQ0L1tgoj0dz9+RkTWuG+Y9G8R6YjrfizPiesezvVE5C4RWSIiK0RkpohEeOznJRFZICJbivbpfu5Bcd3veYWIPOPeVk9EPnHfZOxbEcmo8AOi1Fl0OkVZloi0xnULhpa4vleXAnkez1fHdWl2hjHGiEiMMeaAiMwBPjTGzHC/7oAxZpz78ZO4rgx82b2bJFxX6WbgukJwhohcjWs0384Yc8T9PuBqdHu3MWajiLQDxgBX+O4IKHVhGuLKyi4DZhn3PT/c4ezpIHAMGC8iHwEfnmM/Td3hHQM4gE89nnvfuG62tEZEEt3bugFvFr2vMWa/+66WHYHp8tvtnUPL85dTyhs0xJXVnfO+EMaYQhFpi6uxxI3ASEofGU8A+hhjVojIbUAXj+eOezwWj9/Pft8g4IAxJvMialfK53ROXFnZN0BfEQl336Wxp+eT7tFxtDFmLnAfrntQAxTgaq1XxAnscN9i+OYyvO88YIjH3Hl147qn/FYRud69TUSkxSX/zZTyEg1xZVnG1SZvGq67780Evj3rJU7gQxH5HvgauN+9fSrwZ3F1VKkHPIKrO9NnwLoyvO8nuObHc8XV9PYB91M3A3eISNHdIv3SQlApT3oXQ6WUCmA6EldKqQCmIa6UUgFMQ1wppQKYhrhSSgUwDXGllApgGuJKKRXANMSVUiqA/T+pXff+JXOk0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cincin_kl_divergence.plot(\"distance\", \"divergence\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "48134eda",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cincin_kl_divergence = cincin_kl_divergence.groupby('observation').first().reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "72614e90",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>observation</th>\n",
       "      <th>distance</th>\n",
       "      <th>divergence</th>\n",
       "      <th>population_covered</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>654795.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.142210</td>\n",
       "      <td>69.0</td>\n",
       "      <td>POLYGON ((1398269.926 411448.786, 1398333.955 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>654816.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.060205</td>\n",
       "      <td>49.0</td>\n",
       "      <td>POLYGON ((1397564.804 411646.644, 1397600.475 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>654853.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.210239</td>\n",
       "      <td>0.0</td>\n",
       "      <td>POLYGON ((1397634.522 412227.972, 1397580.224 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>654869.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.189339</td>\n",
       "      <td>44.0</td>\n",
       "      <td>POLYGON ((1397788.495 411495.938, 1397826.597 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>654889.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.210239</td>\n",
       "      <td>142.0</td>\n",
       "      <td>POLYGON ((1397634.522 412227.972, 1397695.312 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>452</th>\n",
       "      <td>739686.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.244836</td>\n",
       "      <td>45.0</td>\n",
       "      <td>POLYGON ((1410127.237 417672.072, 1410268.460 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>453</th>\n",
       "      <td>741823.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.934698</td>\n",
       "      <td>29.0</td>\n",
       "      <td>POLYGON ((1410140.124 416441.192, 1410495.720 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>454</th>\n",
       "      <td>741843.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.945546</td>\n",
       "      <td>256.0</td>\n",
       "      <td>POLYGON ((1408692.225 416938.374, 1408734.116 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>455</th>\n",
       "      <td>741889.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.244836</td>\n",
       "      <td>10.0</td>\n",
       "      <td>POLYGON ((1411000.065 417361.617, 1411163.889 ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>456</th>\n",
       "      <td>741903.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.276574</td>\n",
       "      <td>357.0</td>\n",
       "      <td>POLYGON ((1410268.460 418944.037, 1410337.877 ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>457 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     observation  distance  divergence  population_covered  \\\n",
       "0       654795.0       0.0    0.142210                69.0   \n",
       "1       654816.0       0.0    0.060205                49.0   \n",
       "2       654853.0       0.0    0.210239                 0.0   \n",
       "3       654869.0       0.0    0.189339                44.0   \n",
       "4       654889.0       0.0    0.210239               142.0   \n",
       "..           ...       ...         ...                 ...   \n",
       "452     739686.0       0.0    1.244836                45.0   \n",
       "453     741823.0       0.0    0.934698                29.0   \n",
       "454     741843.0       0.0    0.945546               256.0   \n",
       "455     741889.0       0.0    1.244836                10.0   \n",
       "456     741903.0       0.0    0.276574               357.0   \n",
       "\n",
       "                                              geometry  \n",
       "0    POLYGON ((1398269.926 411448.786, 1398333.955 ...  \n",
       "1    POLYGON ((1397564.804 411646.644, 1397600.475 ...  \n",
       "2    POLYGON ((1397634.522 412227.972, 1397580.224 ...  \n",
       "3    POLYGON ((1397788.495 411495.938, 1397826.597 ...  \n",
       "4    POLYGON ((1397634.522 412227.972, 1397695.312 ...  \n",
       "..                                                 ...  \n",
       "452  POLYGON ((1410127.237 417672.072, 1410268.460 ...  \n",
       "453  POLYGON ((1410140.124 416441.192, 1410495.720 ...  \n",
       "454  POLYGON ((1408692.225 416938.374, 1408734.116 ...  \n",
       "455  POLYGON ((1411000.065 417361.617, 1411163.889 ...  \n",
       "456  POLYGON ((1410268.460 418944.037, 1410337.877 ...  \n",
       "\n",
       "[457 rows x 5 columns]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cincin_kl_divergence = cincin_kl_divergence.merge(cincin_geoseries, left_on = 'observation', right_on = 'ID')\n",
    "cincin_kl_divergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7dec92b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuBUlEQVR4nO3deXxV1b3//9cng4QpQAApJSj4E0XAJEAYBFQqCuktBQdQ7FWgWPGqON3ba8HWorT2ob1eqdSqF78gSFFAWxV7pQo4lRbBIDgwySjmioCEWRmSfH5/nJ2Tk3DYCTSQkLyfj8d5nH0+e621116E/Tl7OHubuyMiInIsCVXdARERqd6UKEREJJQShYiIhFKiEBGRUEoUIiISKqmqO1DZmjVr5m3atKnqboiInFaWLVv2tbs3jzevxiWKNm3akJubW9XdEBE5rZjZ58eap0NPIiISSolCRERCKVGIiEgoJQoREQmlRCEiIqGUKEREJJQShYiIhKpxv6MQEalNlm7KZ/GGnRS6853UFH7U46xKX4YShYhINfT+xp28vWY7Re7s+fYIuZt3cbiwCHcoLHIK3XF3dh44TPFjhTqf1ViJQkSktnhg7krWbd9PnaQEEhOM7LOb0KTeGSQkGAkGiQmGmdGs/hmMvvT/o/4ZiZjZSemLEoWISDWz/1ABn23bxx2XteOeK86r6u7oZLaISHXz8Re7KfLIoaTqQIlCRKSa+eT/9gCQ1bpx1XYkoEQhIlLNrNu+nzMb1qFxvTOquiuAEoWISLWzbvt+2rVoUNXdiNLJbBGRKvTRF7vJ/XwX7k6RO0UO67bt49rs1lXdtahyE4WZnQ/MjgmdA/wSeC6ItwE2A9e6+66gzjjgJqAQuNPd3wjiXYFpQF3gdeAud3czqxO01xXYCVzn7puDOiOAXwTL/rW7Tz/htRWRWu1QQSGf/t8ejhRGNsruRDfORe5Q5rPHTBc5JRvzIqL1DxUUMn/1dnYdOFxSt6h4o3/0Mko+R2I79x+iyEv30wx6npNWNYMUR7mJwt3XAlkAZpYI/B/wMjAWWOjuD5vZ2ODzz8ysAzAM6Ah8F1hgZue5eyHwFDAaeJ9IosgB5hFJKrvc/VwzGwY8AlxnZmnAeCAbcGCZmc0tTkgiUr0UFBaxeecBCopKb0xjN4ylN75Hb0gjG9mSDfXBI0UsWL2NPd8eOaqtkg1ypN1CP3q5hTHTuw4cZu/Bgkpf7+YN69Dpu6nR3zYkmpGQAGZGgkV+95BghgXvsb+DaFr/DG7seTYpZyRG6pmRmGCckVR9zgwc76GnfsAGd//czAYDfYP4dOAd4GfAYGCWux8CNpnZeqC7mW0GUt19MYCZPQdcSSRRDAYeCNp6CXjCIr8cGQDMd/f8oM58IsnlheNdUZHT2c79h/j2SGGcb6exG9j4G2aI3aiWbFBzN+/i850HjvlNtyLfpsuW37b3EF/tPVjp65+aksQ5zRuQGPzYzIKNbVLwYzSL2RjHbpiLN9aJQSw5MYFLzmtO0wZnlCprZTbipTfqpZdZHDODhATDgGYN6lSrDXtlO95EMYySjXQLd98K4O5bzezMIN6KyB5DsbwgdiSYLhsvrvNF0FaBme0BmsbG49SJMrPRRPZUOOusyv/5utQeB48UsvarfdHbI5TdwDpxNtRFJbGv9nzLx3l74h5q8Hgb3WOUid0If3ukiNVb956U9U1vUpfkxIQ4G8nYjWL8jWRigpGcYKXKt0hN4c7z29GkXnLpOqHfruMvIyGhpNxZafWod4ZOqVaVCo+8mZ0BDALGlVc0TsxD4idapyTgPhmYDJCdnX3UfJGKevC1VbywdMs/1UbzhnWod0ZiqY2hUWbjmHD0Bjm2TOxGuEGK0a/9uZzVtF7FN7TBhrm4zaPKJBhnNqzD2U3rV8awSQ13PCn6+8CH7r4t+LzNzFoGexMtge1BPA+IPV2fDnwZxNPjxGPr5JlZEtAIyA/ifcvUeec4+ixSIUcKizhUUMQbK7/i0vOa8+PebUIPTcQ9DGFGneQEzmlW/6Tdc0ekKhxPorie0ucG5gIjgIeD91dj4s+b2WNETma3A5a6e6GZ7TOznsASYDjw+zJtLQaGAG8FV0O9AfzGzJoE5fpT/h6NSLncnQOHC1m4ehsbdhzg6Xc3cLigCIBrs1vT9/wzy2lBpPaoUKIws3rAFcAtMeGHgTlmdhOwBRgK4O4rzWwOsAooAG4PrngCuJWSy2PnBS+AKcCM4MR3PpFzIbh7vpn9CvggKDeh+MS21BwFhUUsXLOdvd8eiVylUkRw9YpTGFzVUnxb5aIiZ9/BAt79bAcHjxRGr3KJXvEScyVMYVHpK1/KXiUTq8+5zejTrhkNU5IY0LFF1QyESDVl7jXrkH52drbn5uZWdTfkODz6xlqeeHv9cdXpenYT0pvUjXOVS8mlh6XmJZRMx17C2K1tE85t3oDmDevocJHUama2zN2z483TZQRSadydjV8f4GDYZZxFpS+9zP/mME+/u4Ers77Lf/Q/n4SEkmvQE4PryYtjxRv/xITIS0RODSWK01RBYRE7DxyOHpo5kevr3Z3DhUW899nXbN93sMwPnYgevom9nLMwzvX0xe3t+uYwG3ccOO51adagDuN/2JEm9avHDdBEpDQlikp24FABhwuKSo6xF29si6CgqKjkGHyRs+ubw/x9/dccLiiKbqALo8fQI3WK2yneSBf/cnXll3vZkv9NpfQ5waBlo7rRHzPFHqY5+pr24PBNzNU+SQmRHz3Vb1KPYd1ac3bT+hW7jDNY3jnNGihJiFRjShSV5N3PdvDcPzazcM328gvHSDBISY78dL/4+vnYY+iR6SAeU6ZFah1G9W5DSnLp6/VjN+iJZTfMx/jh09lN65HepN5JGhkROd0pUVSCgsIi7n3pI7btPcTIXm04K61edIOfGD3mbiQmQGJCQpAAIrcT6NY2jdSU5KpeBRGRY1KiCBwqKGTRuq9jDv1Q6jBQYVHJYZ/YSy0PFxTxwgdb2Lb3EE/f0IWcTi2relVERCqVEkVg38ECbpp+YpfVtmyUwvXdz+Ky9rr+XkRqHiWKQKO6ycwd0zt60rXsid3EmOP/iQmlr8lvmJKsyzVFpMZSoggkJyaQkd64qrshIlLt1NwbqIuISKVQohARkVBKFCIiEkqJQkREQilRiIhIKCUKEREJpUQhIiKhlChERCSUEoWIiISqUKIws8Zm9pKZrTGz1WZ2kZmlmdl8M1sXvDeJKT/OzNab2VozGxAT72pmnwTzJlnw7Ekzq2Nms4P4EjNrE1NnRLCMdWY2ohLXXUREKqCiexSPA3919/ZAJrAaGAssdPd2wMLgM2bWARgGdARygCfNLDFo5ylgNNAueOUE8ZuAXe5+LjAReCRoKw0YD/QAugPjYxOSiIicfOUmCjNLBS4BpgC4+2F33w0MBqYHxaYDVwbTg4FZ7n7I3TcB64HuZtYSSHX3xe7uwHNl6hS39RLQL9jbGADMd/d8d98FzKckuYiIyClQkT2Kc4AdwLNmttzM/p+Z1QdauPtWgOD9zKB8K+CLmPp5QaxVMF02XqqOuxcAe4CmIW2VYmajzSzXzHJ37NhRgVUSEZGKqkiiSAK6AE+5e2fgAMFhpmOId79tD4mfaJ2SgPtkd8929+zmzZuHdE1ERI5XRRJFHpDn7kuCzy8RSRzbgsNJBO/bY8q3jqmfDnwZxNPjxEvVMbMkoBGQH9KWiIicIuUmCnf/CvjCzM4PQv2AVcBcoPgqpBHAq8H0XGBYcCVTWyInrZcGh6f2mVnP4PzD8DJ1itsaArwVnMd4A+hvZk2Ck9j9g5iIiJwiFX1w0R3ATDM7A9gI/JhIkpljZjcBW4ChAO6+0szmEEkmBcDt7l4YtHMrMA2oC8wLXhA5UT7DzNYT2ZMYFrSVb2a/Aj4Iyk1w9/wTXFcRETkBFvniXnNkZ2d7bu6JPftaRKS2MrNl7p4db55+mS0iIqGUKEREJJQShYiIhFKiEBGRUEoUIiISSolCRERCKVGIiEgoJQoREQmlRCEiIqGUKEREJJQShYiIhFKiEBGRUEoUIiISSolCRERCKVGIiEgoJQoREQmlRCEiIqGUKEREJFSFEoWZbTazT8xshZnlBrE0M5tvZuuC9yYx5ceZ2XozW2tmA2LiXYN21pvZJDOzIF7HzGYH8SVm1iamzohgGevMbESlrbmIiFTI8exRfM/ds2KeqToWWOju7YCFwWfMrAMwDOgI5ABPmlliUOcpYDTQLnjlBPGbgF3ufi4wEXgkaCsNGA/0ALoD42MTkoiInHz/zKGnwcD0YHo6cGVMfJa7H3L3TcB6oLuZtQRS3X2xuzvwXJk6xW29BPQL9jYGAPPdPd/ddwHzKUkuIiJyClQ0UTjwppktM7PRQayFu28FCN7PDOKtgC9i6uYFsVbBdNl4qTruXgDsAZqGtFWKmY02s1wzy92xY0cFV0lERCoiqYLlerv7l2Z2JjDfzNaElLU4MQ+Jn2idkoD7ZGAyQHZ29lHzRUTkxFUoUbj7l8H7djN7mcj5gm1m1tLdtwaHlbYHxfOA1jHV04Evg3h6nHhsnTwzSwIaAflBvG+ZOu9UdOVEarsjR46Ql5fHwYMHq7orUk2kpKSQnp5OcnJyheuUmyjMrD6Q4O77gun+wARgLjACeDh4fzWoMhd43sweA75L5KT1UncvNLN9ZtYTWAIMB34fU2cEsBgYArzl7m5mbwC/iTmB3R8YV+G1E6nl8vLyaNiwIW3atCG4yFBqMXdn586d5OXl0bZt2wrXq8geRQvg5eCPLAl43t3/amYfAHPM7CZgCzA06MhKM5sDrAIKgNvdvTBo61ZgGlAXmBe8AKYAM8xsPZE9iWFBW/lm9ivgg6DcBHfPr/DaidRyBw8eVJKQKDOjadOmHO+53HIThbtvBDLjxHcC/Y5R5yHgoTjxXKBTnPhBgkQTZ95UYGp5/RSR+JQkJNaJ/D3ol9kiIhJKiUJETrnNmzfTqdNRBxdOuRUrVvD6669HP8+dO5eHH364CntUPSlRiEiNUFBQcNx1yiaKQYMGMXbs2MrsVo2gRCEiJ91jjz1Gp06d6NSpE7/73e+AyIZ9xIgRZGRkMGTIEL755hsAxo4dS4cOHcjIyOCnP/0pADt27OCaa66hW7dudOvWjb///e8APPDAA4wePZr+/fszfPhwevTowcqVK6PL7du3L8uWLWPp0qX06tWLzp0706tXL9auXcvhw4f55S9/yezZs8nKymL27NlMmzaNMWPGAPD555/Tr18/MjIy6NevH1u2bAFg5MiR3HnnnfTq1YtzzjmHl1566VQNY5Wp6A/uROQ09+BrK1n15d5KbbPDd1MZ/8OOoWWWLVvGs88+y5IlS3B3evTowaWXXsratWuZMmUKvXv3ZtSoUTz55JOMGjWKl19+mTVr1mBm7N69G4C77rqLe+65hz59+rBlyxYGDBjA6tWro+0vWrSIunXrMnHiRObMmcODDz7I1q1b+fLLL+natSt79+7lvffeIykpiQULFnDffffxpz/9iQkTJpCbm8sTTzwBwLRp06L9HjNmDMOHD2fEiBFMnTqVO++8k1deeQWArVu3smjRItasWcOgQYMYMmRIpY5rdaM9ChE5qRYtWsRVV11F/fr1adCgAVdffTV/+9vfaN26Nb179wbghhtuYNGiRaSmppKSksJPfvIT/vznP1OvXj0AFixYwJgxY8jKymLQoEHs3buXffv2AZHDRXXr1gXg2muv5cUXXwRgzpw5DB0auZhyz549DB06lE6dOnHPPfeU2us4lsWLF/OjH/0IgBtvvJFFixZF51155ZUkJCTQoUMHtm3bVkkjVX1pj0Kklijvm//JErkH6NHKXqZpZiQlJbF06VIWLlzIrFmzeOKJJ3jrrbcoKipi8eLF0YQQq379+tHpVq1a0bRpUz7++GNmz57N//zP/wBw//33873vfY+XX36ZzZs307dv3+Nej9j+1qlTp9z1q0m0RyEiJ9Ull1zCK6+8wjfffMOBAwd4+eWXufjii9myZQuLFy8G4IUXXqBPnz7s37+fPXv28C//8i/87ne/Y8WKFQD0798/engIiMbjGTZsGL/97W/Zs2cPF154IRDZo2jVKnI/0djDSw0bNozumZTVq1cvZs2aBcDMmTPp06fPiQ7BaU+JQkROqi5dujBy5Ei6d+9Ojx49+MlPfkKTJk244IILmD59OhkZGeTn53Prrbeyb98+Bg4cSEZGBpdeeikTJ04EYNKkSeTm5pKRkUGHDh14+umnj7m8IUOGMGvWLK699tpo7N5772XcuHH07t2bwsLCaPx73/seq1atip7MjjVp0iSeffZZMjIymDFjBo8//nglj8zpw2rablN2drbn5uZWdTdEqoXVq1dzwQUXVHU3pJqJ93dhZstiHkxXivYoREQklBKFiIiEUqIQEZFQShQiIhJKiUJEREIpUYiISCj9MltETqkHHniABg0asHfvXi655BIuv/zyuOVeeeUVzjvvPDp06HCKeyhlVXiPwswSzWy5mf0l+JxmZvPNbF3w3iSm7DgzW29ma81sQEy8q5l9EsybZMFv4s2sjpnNDuJLzKxNTJ0RwTLWmdmISllrEalyEyZMOGaSgEiiWLVq1SnskRzL8Rx6ugtYHfN5LLDQ3dsBC4PPmFkHIs+87gjkAE+aWWJQ5ylgNNAueOUE8ZuAXe5+LjAReCRoKw0YD/QAugPjYxOSiJweHnroIc4//3wuv/xy1q5dC0Ru1118i+6ytxb/xz/+wdy5c/nP//xPsrKy2LBhA8888wzdunUjMzOTa665Jnpb8rDbfv/2t7/lwgsvJDMzM/qciQ0bNpCTk0PXrl25+OKLWbNmzSkejdNPhQ49mVk68AMiz8H+9yA8GOgbTE8H3gF+FsRnufshYJOZrQe6m9lmINXdFwdtPgdcCcwL6jwQtPUS8ESwtzEAmO/u+UGd+USSywsnsrIitdq8sfDVJ5Xb5ncuhO+HPxFu2bJlzJo1i+XLl1NQUECXLl3o2rVrdH5+fv5RtxZv3LgxgwYNYuDAgdFbeDdu3Jibb74ZgF/84hdMmTKFO+64A4h/2+958+bxyiuvsGTJEurVq0d+fj4Ao0eP5umnn6Zdu3YsWbKE2267jbfeeqtyx6WGqeg5it8B9wINY2It3H0rgLtvNbMzg3gr4P2YcnlB7EgwXTZeXOeLoK0CM9sDNI2Nx6kTZWajieypcNZZZ1VwlUTkVPjb3/7GVVddFb1l+KBBg0rNj721+A9+8AMGDhwYt51PP/2UX/ziF+zevZv9+/czYED0qHbc234vWLCAH//4x9HlpqWlsX//fv7xj39Ebz8OcOjQoUpd35qo3ERhZgOB7e6+zMz6VqBNixPzkPiJ1ikJuE8GJkPkXk8V6KNI7VPON/+TqewtxWMd69biZY0cOZJXXnmFzMxMpk2bxjvvvBOdF++23+5+1HKLiopo3Lhx6N1n5WgVOUfRGxgUHDqaBVxmZn8EtplZS4DgfXtQPg9oHVM/HfgyiKfHiZeqY2ZJQCMgP6QtETlNXHLJJbz88st8++237Nu3j9dee63U/GPdWrzsLcD37dtHy5YtOXLkCDNnzix3uf3792fq1KnRcxn5+fmkpqbStm3b6MON3J2PPvqokta05io3Ubj7OHdPd/c2RE5Sv+XuNwBzgeKrkEYArwbTc4FhwZVMbYmctF4aHKbaZ2Y9g/MPw8vUKW5rSLAMB94A+ptZk+Akdv8gJiKniS5dunDdddeRlZXFNddcw8UXX1xq/rFuLT5s2DD+67/+i86dO7NhwwZ+9atf0aNHD6644grat29f7nJzcnIYNGgQ2dnZZGVl8eijjwKRZ0tMmTKFzMxMOnbsyKuvvlpOS3JctxkPDj391N0HmllTYA5wFrAFGBpz0vnnwCigALjb3ecF8WxgGlCXyEnsO9zdzSwFmAF0JrInMczdNwZ1RgH3BV14yN2fDeujbjMuUkK3GZd4jvc248f1gzt3f4fI1U24+06g3zHKPUTkCqmy8VygU5z4QWBo2Xgwbyow9Xj6KSIilUe38BARkVBKFCIiEkqJQkREQilRiIhIKCUKEREJpduMi8gpVdHbjEv1oUQhIlViwoQJldJOYWEhiYmJ5ReUE6ZDTyJy0oXdZnzevHlce+210bLvvPMOP/zhDwF48803ueiii+jSpQtDhw5l//79ALRp04YJEybQp08fXnzxRV5//XXat29Pnz59uPPOO6M3Fjxw4ACjRo2iW7dudO7cOfor7GnTpnH11VeTk5NDu3btuPfee6PL/+tf/0qXLl3IzMykX79+oe3UFtqjEKklHln6CGvyK/fZC+3T2vOz7j8LLVPebcavuOIKbrnlFg4cOED9+vWZPXs21113HV9//TW//vWvWbBgAfXr1+eRRx7hscce45e//CUAKSkpLFq0iIMHD9KuXTvee+892rZty/XXXx9t+6GHHuKyyy5j6tSp7N69m+7du0cPda1YsYLly5dTp04dzj//fO644w5SUlK4+eabo20V35r8WO3Ur1+/UsezutIehYicVLG3GU9NTT3qNuNJSUnk5OTw2muvUVBQwP/+7/8yePBg3n//fVatWkXv3r3Jyspi+vTpfP7559F61113HQBr1qzhnHPOoW3btgClEsWbb77Jww8/TFZWFn379uXgwYNs2bIFgH79+tGoUSNSUlLo0KEDn3/+Oe+//z6XXHJJtK20tLRy26kNtEchUkuU983/ZAq7zThENvp/+MMfSEtLo1u3bjRs2BB354orruCFF+I/p6z423zY/ercnT/96U+cf/75peJLliwpdWvyxMRECgoK4t6aPKyd2kJ7FCJyUpV3m3GAvn378uGHH/LMM89E9xR69uzJ3//+d9avXw/AN998w2effXZU3fbt27Nx40Y2b94MwOzZs6PzBgwYwO9///toMlm+fHloXy+66CLeffddNm3aBBA99HS87dQ0ShQiclKVd5txiHyjHzhwIPPmzYueiG7evDnTpk3j+uuvJyMjg549e8Z9vnXdunV58sknycnJoU+fPrRo0YJGjRoBcP/993PkyBEyMjLo1KkT999/f2hfmzdvzuTJk7n66qvJzMyMJq3jbaemOa7bjJ8OdJtxkRK15Tbj+/fvp0GDBrg7t99+O+3ateOee+6p6m5VW8d7m3HtUYjIae+ZZ54hKyuLjh07smfPHm655Zaq7lKNopPZInLau+eee7QHcRJpj0JEREIpUYiISKhyE4WZpZjZUjP7yMxWmtmDQTzNzOab2brgvUlMnXFmtt7M1prZgJh4VzP7JJg3yYILls2sjpnNDuJLzKxNTJ0RwTLWmdmISl17EREpV0X2KA4Bl7l7JpAF5JhZT2AssNDd2wELg8+YWQdgGNARyAGeNLPiO3Y9BYwG2gWvnCB+E7DL3c8FJgKPBG2lAeOBHkB3YHxsQhIRkZOv3EThEfuDj8nBy4HBwPQgPh24MpgeDMxy90PuvglYD3Q3s5ZAqrsv9sg1uc+VqVPc1ktAv2BvYwAw393z3X0XMJ+S5CIiIqdAhc5RmFmima0AthPZcC8BWrj7VoDg/cygeCvgi5jqeUGsVTBdNl6qjrsXAHuApiFtle3faDPLNbPcHTt2VGSVRKSG2Lx5M506dSq3zPPPPx/9nJuby5133nmyu1blGjRoUCntVChRuHuhu2cB6UT2DsL+VeLd1MVD4idaJ7Z/k909292zmzdvHtI1EamNyiaK7OxsJk2aVIU9OnEFBQWnfJnH9TsKd99tZu8QOfyzzcxauvvW4LDS9qBYHtA6plo68GUQT48Tj62TZ2ZJQCMgP4j3LVPnnePps4hEfPWb33BodeXeZrzOBe35zn33hZbZvHkzOTk59OjRg+XLl3Peeefx3HPPsXjxYn76059SUFBAt27deOqpp6hTpw5t2rThuuuu4+233wbg+eef59xzz2XkyJEMHDiQIUOGAJFvy8XPp4hd1o033siBAwcAeOKJJ+jVqxdjx45l9erVZGVlMWLECDp37syjjz7KX/7yF/Lz8xk1ahQbN26kXr16TJ48mYyMDB544AG2bNnCxo0b2bJlC3fffXfoXshzzz3Ho48+ipmRkZHBjBkz+Pzzzxk1ahQ7duygefPmPPvsszRq1IjMzEw2btxIQkIC33zzDeeff350Obfffjs7duygXr16PPPMM7Rv356RI0eSlpbG8uXL6dKlC7fddlvccps2beJHP/oRBQUF5ORU3lH6ilz11NzMGgfTdYHLgTXAXKD4KqQRQPGTPOYCw4IrmdoSOWm9NDg8tc/MegbnH4aXqVPc1hDgreA8xhtAfzNrEpzE7h/EROQ0snbtWkaPHs3HH39Mamoqjz32GCNHjmT27Nl88sknFBQU8NRTT0XLp6amsnTpUsaMGcPdd99d4eWceeaZzJ8/nw8//JDZs2dHN+wPP/wwF198MStWrDjqh3njx4+nc+fOfPzxx/zmN79h+PDh0Xlr1qzhjTfeYOnSpTz44IMcOXIk7nJXrlzJQw89xFtvvcVHH33E448/DsCYMWMYPnw4H3/8Mf/6r//KnXfeGU0U7777LgCvvfYaAwYMIDk5mdGjR/P73/+eZcuW8eijj3LbbbdFl/HZZ5+xYMEC/vu///uY5e666y5uvfVWPvjgA77zne9UeNzK5e6hLyADWA58DHwK/DKINyVytdO64D0tps7PgQ3AWuD7MfHsoI0NwBOU3GsqBXiRyInvpcA5MXVGBfH1wI/L62/Xrl1dRCJWrVpV1V3wTZs2eevWraOfFy5c6H379vWLL744GluwYIFfddVV7u5+9tln+4YNG9zd/fDhw56Wlubu7iNGjPAXX3wxWqd+/frR9jt27Oju7rt37/YbbrjBO3Xq5JmZmV63bl13d3/77bf9Bz/4QbRu7OesrKzo8tzd09PTfffu3T5+/Hj/9a9/HY23b9/ev/jii7jrOGnSJL/vvvuOijdt2tQPHz4cXZemTZu6u/vMmTP9lltucXf3K6+80t98803ft2+fp6SkeGZmZvTVvn376LpPmzbN3T20XFpaWnR5e/bsiY5RWfH+LoBcP8Z2tdxDT+7+MdA5Tnwn0O8YdR4CHooTzwWOOr/h7geBocdoayowtbx+ikj1Vd7zKMLKF08nJSVRVFQERL7gHj58+Kh6EydOpEWLFnz00UcUFRWRkpJS7rI8zo1Ri5cZ75kVx2qjIutYXGbQoEGMGzeO/Px8li1bxmWXXcaBAwdo3LgxK1asiFu3+PkbRUVFoeWOd6wrQr/MFpGTbsuWLSxevBiAF154gcsvv5zNmzdHnzUxY8YMLr300mj54mdKzJ49m4suugiIPCd72bJlALz66qtxDwPt2bOHli1bkpCQwIwZMygsLASgYcOG7Nu3L27fLrnkEmbOnAlEntfdrFkzUlNTj2v9+vXrx5w5c9i5cydQ8hyLXr16MWvWLABmzpxJnz59gMj5le7du3PXXXcxcOBAEhMTSU1NpW3btrz44otAJPl89NFHRy0rrFzv3r1LLa+yKFGIyEl3wQUXMH36dDIyMsjPz+eee+7h2WefZejQoVx44YUkJCTwb//2b9Hyhw4dokePHjz++ONMnDgRgJtvvpl3332X7t27s2TJkrjPq77tttuYPn06PXv25LPPPouWycjIICkpiczMzGh7xR544AFyc3PJyMhg7NixTJ8+/ah2y9OxY0d+/vOfc+mll5KZmcm///u/AzBp0iSeffbZ6Mnt4nMXEHmq3x//+MfoMy8gsnGfMmUKmZmZdOzYkVdfffWoZYWVe/zxx/nDH/5At27d2LNnz3Gvx7HoeRQiNVh1eB7F5s2bGThwIJ9++mmFyrdp04bc3FyaNWt2kntWe+l5FCIiUqn0PAoROanatGlT4b0JIPrs6+po586d9Ot39DU8CxcupGnTplXQo1NDiUKkhqvoFTlSvqZNmx7zaqPTxYmcbtChJ5EaLCUlhZ07d57QxkFqHndn586dFbpsOJb2KERqsPT0dPLy8tDNMqVYSkoK6enp5ReMoUQhUoMlJyfTtm3bqu6GnOZ06ElEREIpUYiISCglChERCaVEISIioZQoREQklBKFiIiEUqIQEZFQShQiIhKqIs/Mbm1mb5vZajNbaWZ3BfE0M5tvZuuC9yYxdcaZ2XozW2tmA2LiXc3sk2DepODZ2QTP154dxJeYWZuYOiOCZawzsxGIiMgpVZE9igLgP9z9AqAncLuZdQDGAgvdvR2RZ2aPBQjmDQM6AjnAk2aWGLT1FDAaaBe8coL4TcAudz8XmAg8ErSVBowHegDdgfGxCUlERE6+chOFu2919w+D6X3AaqAVMBgofhTUdODKYHowMMvdD7n7JmA90N3MWgKp7r44eJD3c2XqFLf1EtAv2NsYAMx393x33wXMpyS5iIjIKXBc5yiCQ0KdgSVAC3ffCpFkApwZFGsFfBFTLS+ItQqmy8ZL1XH3AmAP0DSkLREROUUqnCjMrAHwJ+Bud98bVjROzEPiJ1ontm+jzSzXzHJ1l0wRkcpVoURhZslEksRMd/9zEN4WHE4ieN8exPOA1jHV04Evg3h6nHipOmaWBDQC8kPaKsXdJ7t7trtnN2/evCKrJCIiFVSRq54MmAKsdvfHYmbNBYqvQhoBvBoTHxZcydSWyEnrpcHhqX1m1jNoc3iZOsVtDQHeCs5jvAH0N7MmwUns/kFMREROkYo8j6I3cCPwiZmtCGL3AQ8Dc8zsJmALMBTA3Vea2RxgFZErpm5398Kg3q3ANKAuMC94QSQRzTCz9UT2JIYFbeWb2a+AD4JyE9w9/8RWVUREToTVtEckZmdne25ublV3Q0TktGJmy9w9O948/TJbRERCKVGIiEgoJQoREQmlRCEiIqGUKEREJJQShYiIhFKiEBGRUEoUIiISSolCRERCKVGIiEgoJQoREQmlRCEiIqGUKEREJJQShYiIhFKiEBGRUEoUIiISSolCRERCKVGIiEiochOFmU01s+1m9mlMLM3M5pvZuuC9Scy8cWa23szWmtmAmHhXM/skmDfJzCyI1zGz2UF8iZm1iakzIljGOjMbUWlrLSIiFVaRPYppQE6Z2Fhgobu3AxYGnzGzDsAwoGNQ50kzSwzqPAWMBtoFr+I2bwJ2ufu5wETgkaCtNGA80APoDoyPTUgiInJqlJso3P09IL9MeDAwPZieDlwZE5/l7ofcfROwHuhuZi2BVHdf7O4OPFemTnFbLwH9gr2NAcB8d893913AfI5OWCIicpKd6DmKFu6+FSB4PzOItwK+iCmXF8RaBdNl46XquHsBsAdoGtLWUcxstJnlmlnujh07TnCVREQknso+mW1xYh4SP9E6pYPuk909292zmzdvXqGOiohIxZxootgWHE4ieN8exPOA1jHl0oEvg3h6nHipOmaWBDQicqjrWG2JiMgpdKKJYi5QfBXSCODVmPiw4EqmtkROWi8NDk/tM7OewfmH4WXqFLc1BHgrOI/xBtDfzJoEJ7H7BzERETmFksorYGYvAH2BZmaWR+RKpIeBOWZ2E7AFGArg7ivNbA6wCigAbnf3wqCpW4lcQVUXmBe8AKYAM8xsPZE9iWFBW/lm9ivgg6DcBHcve1JdREROMot8ea85srOzPTc3t6q7ISJyWjGzZe6eHW+efpktIiKhlChERCSUEoWIiIRSohARkVBKFCIiEkqJQkREQilRiIhIKCUKEREJpUQhIiKhlChERCSUEoWIiIRSohARkVBKFCIiEkqJQkREQilRiIhIKCUKEREJpUQhIiKhlChERCTUaZEozCzHzNaa2XozG1vV/RERqU2qfaIws0TgD8D3gQ7A9WbWoWp7JdWJFxRwIs9+L/r2W7yg4CT0SCrqyFdf4YWFVd0NKUdSVXegAroD6919I4CZzQIGA6sqcyG7tm/h08HfB469wTGHhKLIdGEiFCWUFLeYMpQzfdS84omYmFvJqyghrFfVQLACHkwmFkZeFmy83QwsMt+Ly1rJdIMDzqE6cCTJSuYF7TqQctixIihIssiYl9For/NtCiQVwqEzjMPJFety2m6nMBG+STEKE4/+NzqWUutwjLGIli0zK94yYkOJRZB8xDl8hlFUpq2ybVeaSvrjMqDOYSf5SOTfIaHIKUo0ChIj/2+K17343z6xCBrvdfbVh2/rWKl2SrVbTv8q8u/2z7YRb+gLEyLbgeK/e/z4tgNHbQMcrPgfw0tvF8KmY9vZ9d0ULn99efjKnIDTIVG0Ar6I+ZwH9IgtYGajgdEAZ5111gktJCEpmd3pdaDwSGi54g1VQlHkVbxBg5j/bzF/VR7nL8zD5puBO+aUelVbZRIllCRRNyuZF6wTUPrdYUcdI7HAI0nYS/7TFNfdnQxFiUbSEY87FtvbGMmHIxulaDsV8PVZhhskHXESgy+18f69yiru91FFy34ZiFO3bHIpW8YNCpLjrGuctiv1z6KSktCRZKMgGZIPR3pXnCBivxjE/vtvPcNIPuRHj8Nx9q8i/27HW6e8+QlFRP9uoPSXG4LpY7V1zG2AWan65bVVtt2kMxuGd/oEnQ6JIu7/t1If3CcDkwGys7NP6P9Po7SW/PDFD0+kqohIjVbtz1EQ2YNoHfM5HfiyivoiIlLrnA6J4gOgnZm1NbMzgGHA3Cruk4hIrVHtDz25e4GZjQHeABKBqe6+soq7JSJSa1T7RAHg7q8Dr1d1P0REaqPT4dCTiIhUISUKEREJpUQhIiKhlChERCSUncg9cqozM9sBfP5PNNEM+LqSunM60ziU0FiU0FiUqGljcba7N483o8Ylin+WmeW6e3ZV96OqaRxKaCxKaCxK1Kax0KEnEREJpUQhIiKhlCiONrmqO1BNaBxKaCxKaCxK1Jqx0DkKEREJpT0KEREJpUQhIiKhlCgCZpZjZmvNbL2Zja3q/pxsZjbVzLab2acxsTQzm29m64L3JjHzxgVjs9bMBlRNryufmbU2s7fNbLWZrTSzu4J4bRyLFDNbamYfBWPxYBCvdWNRzMwSzWy5mf0l+Fwrx0KJgsgfA/AH4PtAB+B6M+tQtb066aYBOWViY4GF7t4OWBh8JhiLYUDHoM6TwZjVBAXAf7j7BUBP4PZgfWvjWBwCLnP3TCALyDGzntTOsSh2F7A65nOtHAsliojuwHp33+juh4FZwOAq7tNJ5e7vAfllwoOB6cH0dODKmPgsdz/k7puA9UTG7LTn7lvd/cNgeh+RjUIraudYuLvvDz4mBy+nFo4FgJmlAz8A/l9MuFaOhRJFRCvgi5jPeUGstmnh7lshsgEFzgzitWJ8zKwN0BlYQi0di+BQywpgOzDf3WvtWAC/A+4FimJitXIslCgiLE5M1w2XqPHjY2YNgD8Bd7v73rCicWI1ZizcvdDds4g8m767mXUKKV5jx8LMBgLb3X1ZRavEidWIsQAlimJ5QOuYz+nAl1XUl6q0zcxaAgTv24N4jR4fM0smkiRmuvufg3CtHIti7r4beIfI8fbaOBa9gUFmtpnIoejLzOyP1M6xUKIIfAC0M7O2ZnYGkZNSc6u4T1VhLjAimB4BvBoTH2ZmdcysLdAOWFoF/at0ZmbAFGC1uz8WM6s2jkVzM2scTNcFLgfWUAvHwt3HuXu6u7chsj14y91voBaOBZwmz8w+2dy9wMzGAG8AicBUd19Zxd06qczsBaAv0MzM8oDxwMPAHDO7CdgCDAVw95VmNgdYReQqodvdvbBKOl75egM3Ap8Ex+YB7qN2jkVLYHpwtU4CMMfd/2Jmi6l9Y3EstfHvQrfwEBGRcDr0JCIioZQoREQklBKFiIiEUqIQEZFQShQiIhJKiUJEREIpUYiISKj/H+MPa1GUuUafAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.pyplot import plot\n",
    "\n",
    "cincin_kl_divergence.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ccda2a",
   "metadata": {},
   "source": [
    "## How the function works\n",
    "Here the code will be broken down into its constituent parts, in order to explain how it works to create the complete `kl_divergence_profile` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2f0afa1",
   "metadata": {},
   "source": [
    "### The bigger picture\n",
    "\n",
    "The methodology suggested by [Olteanu et al. (2019)](https://doi.org/10.1073/pnas.1900192116) inventively calculates consecutive divergence levels in relation to a chosen area within a city.\n",
    "\n",
    "To do this, the KL divergence calculations consectively aggregate more and more population. This aggregation starts by including solely the population of this chosen area, then including populations of the closest neighbour to the chosen area, followed by the next closest neighbour, and so on until the population of the whole city is analysed.\n",
    "\n",
    "The result is a profile of divergence levels that represents an aggregation of population gradually from this chosen area, showing how the dynamics of a population change as you consider areas covering a larger and larger distance from this chosen area.\n",
    "\n",
    "Once this has been undertaken for a single chosen area, it can be repeated for every other area within the city.\n",
    "\n",
    "When analysed on this scale, this gives a more nuanced picture of which areas of a city closest represent the total population of the city, and vice versa which areas are most different to the total population of the city; along with their spatial relationships."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84822e09",
   "metadata": {},
   "source": [
    "### Managing inputs\n",
    "First, as there are different forms of input possible, they need to be managed into a consistent form that will eventually be used to calculate the KL divergence."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d030cd1b",
   "metadata": {},
   "source": [
    "As a reminder, this function currently has three arguments that you can use to tailor it to your own data:\n",
    "\n",
    "```python\n",
    "def kl_divergence_profile(populations, coordinates = None, metric = 'euclidean'):\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b63e90f",
   "metadata": {},
   "source": [
    "Where:\n",
    "- ***populations*** - is data regarding each group within the population information to be included in the analysis (in raw numbers, not percentages)\n",
    "- ***coordinates*** - is the spatial information, relating to the population information  inputted within the *populations* argument\n",
    "- ***metric*** - the style of distance metric to be used in calculating the distances between areas, such that they can be ordered"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d50ce74d",
   "metadata": {},
   "source": [
    "#### Storing indices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b632388d",
   "metadata": {},
   "source": [
    "If the *populations* input contains an [index](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.set_index.html), then it is stored, such that the results can be related to the inputted data more easily. If not, an array is created to take this place."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e65b77f8",
   "metadata": {},
   "source": [
    "```python\n",
    "    # Store the observation index to return with the results\n",
    "    if hasattr(populations, 'index'):\n",
    "        indices = populations.index\n",
    "    else:\n",
    "        indices = np.arange(len(populations))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81443f75",
   "metadata": {},
   "source": [
    "#### Creating *coordinates* from a geodataframe input to *populations*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9de9ddd2",
   "metadata": {},
   "source": [
    "If the *populations* input has a geometry attribute, and there is no *coordinates* input, then the geometry of the *populations* input is passed as the *coordinates* input. The *populations* input is then stripped of its geometry attribute, and converted into an array."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f77e2bd6",
   "metadata": {},
   "source": [
    "```python\n",
    "    # Check for geometry present in populations argument\n",
    "    if hasattr(populations, 'geometry'):\n",
    "        if coordinates is None:\n",
    "            coordinates = populations.geometry\n",
    "        populations = populations.drop(populations.geometry.name, axis = 1).values\n",
    "    populations = np.asarray(populations)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afcdc4dc",
   "metadata": {},
   "source": [
    "#### Ensuring *coordinates* is centroid's coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1760f7d2",
   "metadata": {},
   "source": [
    "Regardless whether the *coordinates* input has a point or polygon geometry attribute, the `centroid` of each observation is calculated and re-stored as *coordinates* in a two column array, representing the *x* and *y* values.\n",
    "\n",
    "A test is inputted into this using a *else* clause, whereby if the *coordinates* and *populations* inputs are not of equal length, then the user is warned of this."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbeee169",
   "metadata": {},
   "source": [
    "```python\n",
    "    #  Creating consistent coordinates - GeoSeries input\n",
    "    if hasattr(coordinates,'geometry'):\n",
    "        centroids = coordinates.geometry.centroid\n",
    "        coordinates = np.column_stack((centroids.x, centroids.y))\n",
    "    #  Creating consistent coordinates - Array input\n",
    "    else:\n",
    "        assert len(coordinates) == len(populations), \"Length of coordinates input needs to be of the same length as populations input\"\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11b1773",
   "metadata": {},
   "source": [
    "### Creating a distance matrix\n",
    "Having managed the inputs, we have to think about how we want the KL divergence levels to be calculated.\n",
    "\n",
    "To create this a divergence profile over space, we first need to be able to sort the inputted data in terms of their proximity to the chosen area of the analysis.\n",
    "\n",
    "A distance matrix is therefore created, to store the distances between each observation and every other observation. This uses the [`pdist`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html) function from *scipy.spatial.distance* to compute this, using the optional *metric* argument from the input of the function.\n",
    "\n",
    "This calculates the Euclidean (also known as Pythagorean) distance between all the points by default, but could be altered as desired by the user using the *metric* input."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d283b7dc",
   "metadata": {},
   "source": [
    "``` python\n",
    "    # Creating distance matrix using defined metric (default euclidean distance)\n",
    "    dist_matrix = squareform(pdist(coordinates, metric = metric))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54fb420c",
   "metadata": {},
   "source": [
    "### Creating a list for the results\n",
    "This is purely a piece of admin to open a list, such that the results can be filled into this (currently empty) list later in the function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32e594c9",
   "metadata": {},
   "source": [
    "``` python\n",
    "    # Preparing list for results\n",
    "    results = []\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "484498ab",
   "metadata": {},
   "source": [
    "### Undertaking the loop"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f701536",
   "metadata": {},
   "source": [
    "#### Part A - Sorting the population information by the distance matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d54b76b",
   "metadata": {},
   "source": [
    "In order to undertake all of the statistical work with regards to the *populations* input, the data must first be ordered by the distance matrix, which has been created above.\n",
    "\n",
    "For instance - when the *metric* is set to Euclidean - for an example observation '*A*', sorting by the distance matrix will return the *populations* information ordered from closest observation to the furthest observation, all with respect to observation *A*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7dadb4f",
   "metadata": {},
   "source": [
    "``` python\n",
    "    # Loop to calculate KL divergence\n",
    "    for (i, distances) in enumerate(dist_matrix):\n",
    "\n",
    "        # Creating the q and r objects\n",
    "        sorted_indices = np.argsort(distances)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e756902",
   "metadata": {},
   "source": [
    "The `sorted_indices` can then be used as a method of ordering the *populations* object, which we will see in the next part of the loop."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41ae92b9",
   "metadata": {},
   "source": [
    "#### Part B - Preparing the inputs to calculate the KL divergence\n",
    "With the end goal of this function being is to calculate the KL divergence - specific inputs for this calculation must be created."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0db6235",
   "metadata": {},
   "source": [
    "In the article, Olteanu et al. state that\n",
    "> $\\text{\"}$ The KL divergence between two multinomial distributions (*$\\space q_1, \\space ... \\space , q_k$* ) , (*$\\space r_1, \\space ... \\space , r_k$* ) is defined as:\n",
    "> $$ d_{KL}(q \\space|\\space r) = \\sum_{i=1}^{k} q_i \\space ln \\space{{q_i} \\over {r_i}} \\space \\space \\text{\"}  $$\n",
    "> [Olteanu et al. (pg. 12250, 2019)](https://doi.org/10.1073/pnas.1900192116)\n",
    "\n",
    "When applied in practice, we create statistics from the distributions such that\n",
    "\n",
    "$$ ln \\space {{q_i} \\over {r_i}} \\space = \\space ln \\space {{{{cumulative \\space proportions}\\space_i}} \\over {{total \\space proportions}\\space_i}} $$ "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "230dbde8",
   "metadata": {},
   "source": [
    "The *culmulative proportions* and *total proportions* objects are the inputs that will be used to calculate the KL divergence for each area.\n",
    "\n",
    "- **cumulative proportions** - culmulatively aggregating with respect to the chosen area, the proportion (percentage) that each group represents\n",
    "- **total proportions** - the proportions (percentages) each group represents in the total population, placed into an ordered list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16bba0c8",
   "metadata": {},
   "source": [
    "Therefore, to calculate *cumulative proportions* (`q_cumul_proportions`), the objects representing the culmulative group populations (`cumul_pop_by_group`) and culmulative total populations (`obs_cumul_pop`) are created."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "855fdcc7",
   "metadata": {},
   "source": [
    "``` python\n",
    "        cumul_pop_by_group = np.cumsum(populations[sorted_indices], axis = 0)\n",
    "        obs_cumul_pop = np.sum(cumul_pop_by_group, axis = 1)[:, np.newaxis]\n",
    "        q_cumul_proportions = cumul_pop_by_group / obs_cumul_pop\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a44f7f2",
   "metadata": {},
   "source": [
    "To calculate *total proportions* (`r_total_proportions`), the objects representing the total group populations (`total_pop_by_group`) and total populations (`total_pop`) are created"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bdb98ef",
   "metadata": {},
   "source": [
    "``` python\n",
    "        total_pop_by_group = np.sum(populations, axis = 0, keepdims = True)\n",
    "        total_pop = np.sum(populations)\n",
    "        r_total_proportions = total_pop_by_group / total_pop\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92d7e123",
   "metadata": {},
   "source": [
    "#### Part C - Calculating the KL divergence\n",
    "To calculate the KL divergence profile, the two objects created to form the *q* and *r* elements of the formula are inputted into the [`scipy.special.rel_entr`](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.special.rel_entr.html) function.\n",
    "\n",
    "While a [`scipy.special.kl_div`](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.special.kl_div.html#scipy.special.kl_div) also exists - and might be thought to be the obvious choice - the relative entropy function contains beneficial functionality for this use-case as a result of its handling of $ln(0)$, where it bends the rule of mathematics to return $0$ instead of $- \\infty$. For measuring population statistics, this makes a lot more logical sense."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb1d45c8",
   "metadata": {},
   "source": [
    "``` python\n",
    "        # Input q and r objects into relative entropy (KL divergence) function\n",
    "        kl_divergence = relative_entropy(q_cumul_proportions,\n",
    "                                         r_total_proportions).sum(axis = 1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e43b888",
   "metadata": {},
   "source": [
    "####  Part D - Placing the outputs into a dataframe\n",
    "Now we have undertaken all the calculations, a dataframe for the outputs can be constructed, returning the four outputs we have [already discussed](#1.-kl_divergence_profile(populations-=-geodataframe))."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81465ab2",
   "metadata": {},
   "source": [
    "``` python\n",
    "        # Creating an output dataframe\n",
    "        output = pd.DataFrame().from_dict(dict(\n",
    "            observation = indices[i],\n",
    "            distance = distances[sorted_indices],\n",
    "            divergence = kl_divergence,\n",
    "            population_covered = obs_cumul_pop.sum(axis=1)\n",
    "        ))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e106d52b",
   "metadata": {},
   "source": [
    "### Bringing all the dataframes together\n",
    "Finally, the *output* dataframes relating to each observation are brought together and concatenated (stacked) using the [`pandas.concat`](https://pandas.pydata.org/docs/reference/api/pandas.concat.html) function to return our full output. Designed to be deconstructed, broken down and used for visualisations and further analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd1a082d",
   "metadata": {},
   "source": [
    "``` python\n",
    "        # Append (bring together) all outputs into results list\n",
    "        results.append(output)\n",
    "\n",
    "    return(pd.concat(results))\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}