safegraph-comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mjUR2XYBtSO3"
   },
   "source": [
    "![alt text](https://global-uploads.webflow.com/5baafc2653bd67278f206724/5be267a03f7813daf821b31e_safegraph-logo-hidpi%403x-p-500.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZE97HH76SArR"
   },
   "source": [
    "# This notebook is based on a notebook provided by SafeGraph and modified by Alex Berke\n",
    "\n",
    "It was modified to measure the correlation / reporting bias of devices at the census tract level, constrained to 3 census tracts in MA. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xsNNli6GTN6s"
   },
   "source": [
    "# \"What about bias in your dataset?\" \n",
    "\n",
    "## Quantifying Sampling Bias in SafeGraph Patterns\n",
    "\n",
    "--------------\n",
    "**[Ryan Fox Squire](https://www.linkedin.com/in/ryanfoxsquire/) | Data Scientist, [SafeGraph](https://safegraph.com/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias)**\n",
    "\n",
    "\n",
    "Oct 2019\n",
    "\n",
    "--------------\n",
    "*Share this notebook: [Shareable Link](https://colab.research.google.com/drive/1u15afRytJMsizySFqA2EPlXSh3KTmNTQ#offline=true&sandboxMode=true\n",
    ")*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oBADlZgHTsWg"
   },
   "source": [
    "## Introduction: \"What about bias in your dataset?\"\n",
    "\n",
    "SafeGraph Patterns measures foot-traffic patterns to 3.6MM commercial points-of-interest from over 45MM mobile devices in the United States and provides a monumental window into American commerce. SafeGraph data users look through this window to ask detailed questions about consumer behavior. \n",
    "<br/><br/>\n",
    "\n",
    "![alt text](https://user-images.githubusercontent.com/6432906/67054597-f33e9a80-f0f9-11e9-81e2-777d78ef36d9.png)\n",
    "\n",
    "###### Image Credit: [123RF Stock Photo ](https://www.123rf.com/profile_elenabsl)\n",
    "<br/><br/>\n",
    "\n",
    "A common question we hear from SafeGraph Patterns customers is “What about bias in your dataset?”. \"Does your panel really represent the true American public?\". \"How do we know that your panel isn't oversampling wealthier people?\".  This is the kind of sophisticated data skepticism we ***love*** to hear.  A key part of [SafeGraph’s vision](https://blog.safegraph.com/safegraph-vision-and-values-circa-2017-d188d4185df2/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias) is to *“Seek the Truth About the World...Of course, data can never be 100% true, but we should strive to make it 100% true.”* \n",
    "\n",
    "And although SafeGraph Patterns aggregates data from ~ 10% of devices in the United States (a very impressive sample, if we don't say so ourselves!) this sample is not a perfect representative subset of the population. Like all samples, the SafeGraph dataset has [sampling error](https://en.wikipedia.org/wiki/Sampling_error). \n",
    "<br/><br/>\n",
    "\n",
    "### This is a technical blog post designed for data scientists and analysts to answer the following questions:\n",
    "* What is Sampling Bias and why is it a problem?\n",
    "* How biased is the SafeGraph Patterns dataset? How do we quantify sampling bias?\n",
    "\n",
    "We also include a short preview of a future post:\n",
    "* How do we correct the sampling bias to answer questions about consumer behavior? (e.g., What is a brand's true customer demographic? Or how far do people travel to visit my stores?)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YDKLdlIk79-m"
   },
   "source": [
    "## What is Sampling Bias?\n",
    "\n",
    "\n",
    "\n",
    "[Sampling error](https://en.wikipedia.org/wiki/Sampling_error) is the difference between a sample and the population. Here, the “sample” is the panel of devices aggregated into SafeGraph Patterns and the “population” is the US Population. There are many possible causes of sampling error, but usually the most worrisome is **sampling bias**. Sampling bias is problematic because, unlike other sources of errors which diminish as your sample size increases, sampling bias is a systematic error, not helped by increasing sample size.\n",
    "<br/><br/>\n",
    "\n",
    "[Sampling bias](https://en.wikipedia.org/wiki/Sampling_bias) goes by multiple names (e.g., selection bias, ascertainment bias... these terms have nuanced technical differences, but are generally used synonymously). **Sampling bias is sampling error that arises from the collection of data from sub-groups of the population disproportionately.**  The intuition is simple: you want to sample from sub-groups of the population proportional to their presence in the population. If your population is 80% women and 20% men, then a true random sample should also approximately be 80% women and 20% men. If your sample is actually 50% men and 50% women, then your sample is biased in favor of men. \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MXRIRdBE-Puu"
   },
   "source": [
    "\n",
    "## Why is Sampling Bias (sometimes) a Huge Problem?\n",
    "\n",
    "Sampling bias can be insidious and is one of the most common types of mistakes in statistical analysis of survey or panel data. \n",
    "<br/><br/>\n",
    "\n",
    "**Example Scenario:** \n",
    "\n",
    " <img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAT4AAACfCAMAAABX0UX9AAABL1BMVEXo2u7///81h8ezs7P7rxqZmZnAwMDp6ena2tq6urrV1dXw8PDk5OTZzt5lZWVpaWmBgYGUlZP7uDgqrF//+Oz7rgDNzc3Sx9b5+fkngcWmx+T5kA3E1+qhoaH7rBn5lg94q9elnpP6ohVzc3PHx8e4z+f8x2/92qr7syX8zXz94bb8ukD+6svl2e7x3/DN1t6ltsRFjsrAxuRlmtAOp1O2w8+ApNSetdzV0emNrNiuveCu277U3OLAy9U+jMnIwspdls60r7Zgi67ZdGLFgXaHlaOsj4v6XD3sZ061hn6lrLGMuN/h7PZrwYo8s22c1LB9yJnP6tjp9u7X7uBhvYSi17V6x5c7smq+48vS4vOgmn+wnm/905DPpVX90YjjqUVSjLpyiZ+pkXz+8t3Fo2A9zkQKAAAJeklEQVR4nO2d/V/jthnAo0CcV9ZCiM+mB2U7097KJZDEQICDS9muWwvkjrTrwdqt3br//2+YLPlFSmJbipWz5Or7i7BjG/z9PI8eWbZDqaTRaDQajUaj0Wg0Go1Go9FoNBqNRpOZXq9XKnW7Xf8nDTu93mhwdjocj/v9/ng8PD0bjEraIRO90uhsWLYh5QC0MDyDCvP+42SnNzobE+JIbLt/PtICkxgMY9wFBscXOgRj6F2Mk9z5Bvtnef+dUtIb9NPlIYHXFzoAZ+h1T9nkIYFD3QdS9C7K7PY8gTqDSThCLwjAbt5/szSMGHs9yl95oBPYozfidocF6gpS8iouf+hpfwG9wZLyPH9nv3d/y8eejj/Ikv1e6O/3XT+619n0lcujvE+Bhw3BDLOkLqKvzPhv4+tNwdxktgfHz4qk70Zls8LJX/6a+PHbzPLKypSPDW55lW/+9vfEzzN3fBglur+vue1VKt9+k/TpTbaoCyZW7VMFwq+LznhNIN9lcVce3t7u+1PT9iBvOemg3N0FArldum7Y4+f4ENYBmmwYyx9+z0Tru1ve3gGYXt5P2pP7SwB2bCXCT7g+KvhsHu7Au7bb9nDbj+DJViH8EvU1WxFVvGaX6ufmd7FIe9c7z9m5A/dYHhL4Hjy3FQi/JH1WjaTlrarSw5z5fXZIfXdccRvZm7ht950Xx9IX3yR9zZoVLbRqTA76RObuc9l7dKPIu4T+Hqyy/GM/Zn1VJn1k4bD3mzz6Jr689uTRc+lOYPjZF3n7SYE3+pqJWFTh2D8+Yrf3gILPnVxO0eL3bXcKq8cwbz8pcOqju8N5iNz1kvfkxFp86HneIX2X4fJ79xKWIVvyiRex0UcN+ry+7+i4xajvR8/eNFqewOIh/9BPbOl4PqsPVFkT+B5mrmfPgsDmBxdWEHiMc7lrL68+K4mZMTOqvBZjAr9zYbiBx/u267a91G27/wDyj5zF9n3j8pw+AF4fv2bQB2vtFPyI68ej1xG6D6grkLvz49VXTaBZvV6kDzSPr9LkWc2p+8Eb7uHBC2rxEFzukZ/QcR89XUAMm6+Ok4eAMMPB+w+gTQBz1yvjkteOZH3ESbOUjqc4fSkJDOvLSRO4pL32BNchyZ+54rzmTe77jmP1gebJVWwF8UY31WM4UiaZWvggcpfexBkXa27GhVwzR3UnXp+XwNXFvwSX5qMTMP0QXvZOpmCsvD5eEvXFDAFbwdojaPFy4nr88Aju/IfbJL/j8RH1LRwCHkUxifrH6cPD93D0vBMeRAV922JI0Td/DUcPaayr46PXzbung/3oNQYV9P2xI4Z/pumbSeC5emy1jq5uyReP1NC3HgOhJm4TcuufUvXBCDtpzv9IckAfRGF9nS8iDrHAPQp+fVHIxV2LzOhToPLG6Hu5TZ7Wz56/vc8oZvT9i0Wf581KGEjP6JP7SdMEfZ2fwZd/Dvg32EbB9yeKmT1+YdKHhoBXsfMwM/rknq5P0vcl2Aq7vr1X4GVMikdsLZ4ymOd1wizgjD4Fpgxi9R36Px0edlj0dX5l1AcSpgBpfSrMuCTo8zq8zhYAvr6Zvm+W/7DqS4DW15e660vT14EFowNLyG+Bvs8T+e/8bHNGfXKPW9L6vq/AK2yPLXl/Ea1P8sqRpu8QgFe/QXvrTPrWO6KTV/LKkZq80B/YhuLYkvfzt7ZYff28/aSQWjq2ttGAj610fPaTWH2SX3MwDFw6h1uEvjRekiO/4VL6yGkb2XM3NXm5hs1wHyp7n5awZ5HPeUh+lzdZH7xQ+yoEbKNr3pS+b4+uvbc73JDhK3vdTZlxeUWGxRcLpgzm2OvcUMWXH3J3ye+Rl9ImrLYi1hdNWM2x3vlf9heyQvlyz7Z4iJwuxbvcpGphRvrg8/VtCURY+CkQfOJfTABgX5A/2YfMHivQt/yLMRSSP92CWYG+mdu9y9qTfK4Fw/xaTAvNcDYruyQx+2V/mbxcvpa/bpTSHhEi7OFHhJq7LPqaAoJPhdQV/1oM5ilr+KlQdT0EvxYTcJDNnxodX4n/4dxW0tO51egJluXf6i0r9D0Qol+LiXbI4M+W/V2iCKEvJljk/celR8/2OG8p7KymdCCWHP7Zp3k74UDsS1n0E1NL1Q/p5+cpVtb3edxd8wtU49tvAlYZfdA2Zwdo9xX77twVjftCDri+OlepxPVgfS3G4q68UQAyClTxa5vZX4upzq9K6/swd0wZbI9V/MrD5NdimsT7fnhN4kVHdfGLL+Do1zR3UF7eJpZiFfN9c1Q3Kzf92E5Q5f848bH0VSpvz/uz9yHRnczr80FXUXm+vkrKcCQrNfRLSr3u4Ox0fB3d0h0Pzwcjpf/bjqdv8w8r5wXU500fQ1fd7mjgMRp11f9HRVDf7uaLkE9EEh32xae7WF/ReEZ/d+GnIqG/VXG3iPrqIr84MpmNvM9VPBsfz97am7xPVjxh8BmGiW6oGQ5qTAPdZasFa43l1xrh7yhc9obBVzcMx4Q4UKNHuLh4LbWYuCtcNAobfkTwGegmkNNAd3zwgmmihXrNaxoGuioz0EINLzj+Bgu2NoPLOJMIv4L1fmTw4cyr13Di4TRsoPSr42R0iHUNvIFD7uSvwweqGeHt9eKGHxl8ZkDdzz0/9ch19QXrqO38vJ35KNRXrN4vKruwvwqSzc9LnIpVp0YkM/2RwfxRlL2FCr/wrLzC0fAxnBpqTLRUR41poAXH8LcgG3+L+qIt8DEKGn5k8BnRjKeJf3Yc1BgmbhqoqaOmUac2TNje36CQ4WcQwRcWSph5DSIBgzzExdVv/Ez1G78kNxwq6YkDksVjLe+TFgYVfE5YObzu33GCxgkaM6Yx4z+LjkjoK0z4vSH10VPuQQD6g8D4paoftgYVc2TozQ79ihJ+UfCFg75wqIYHckGDr75a9QYx8GsFV2i4cYJLPbKJjlfA8KOCz6+vfkNV1nA5BnN2OwNDLRQu/Kjgc8hSEJQJv3YEy7GEqerns4lKsH/VVjNxjS5a+NHB5w89yLGbH4vhoC8WJwhOg9rf3ysMwrUihV83Ops6JPwhjrV44IdJe84foQDh9yZeR4yjJeUuIu+Tz0w3/Rxpe0YinEdTPvx4gy/ZHq8+1cOPN/jS9PFm77O8BWSDN/hE61M7/LiDT7g+pcOPO/hSKge/PqXDj/tkxaPyTSOHZZi7WvJWkIUVP83HgupPVmk0Go1Go9FoNBqNRqPRaDQajUaj0Wg0Go1m5fwffToFJ+D06WMAAAAASUVORK5CYII=\" alt=\"Email Survey\" ALIGN=\"right\">\n",
    "\n",
    "> You manage a national retail business, and you email a survey to all your customers. When you analyze the results, you notice that most of the respondents live in California. In reality, most of your stores are in the Midwest. For whatever reason, your survey has significant sampling bias in favor of California. Customers in California may be different from your customers in the midwest in important ways (e.g. demographics, income, retail alternatives, product preferences, etc.). This is a problem. If you applied the results of your survey to your Midwest locations, then you would have a skewed view which could lead to bad (and expensive) business decisions (like stocking the wrong product mix or launching the wrong marketing campaign or buying the wrong real estate).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "How do you know whether sampling bias is a problem in your study? In general, if an independent variable (e.g. a consumer's home state) shows a significant correlation with your outcome variable (e.g. product preferences), then the distribution of this variable in your sample should be similar to the distribution of this variable in the whole population. If it isn’t, then you have a problematic sampling bias.\n",
    "\n",
    "\n",
    "A [simple random sample](https://en.wikipedia.org/wiki/Simple_random_sample) generally solves sampling bias (and is why random assignment to groups is a critical feature of AB tests). But often in the real world (like the email survey example above) random sampling during data collection isn't possible. You have to work with the data you get, and sometimes that data has sampling bias. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pvPYDci4Dnf-"
   },
   "source": [
    "# How to Quantify Sampling Bias in SafeGraph Patterns\n",
    "\n",
    "We want to measure the (sampling) bias in the SafeGraph Patterns dataset so that we can answer questions like \"How many people visited this store?\" and \"Who are the customers of my competitor?\" without getting a biased or incomplete picture. \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2PgWXrESeVvE"
   },
   "source": [
    "\n",
    "## What is SafeGraph Patterns?\n",
    "SafeGraph Patterns is a dataset about commercial brick-and-mortar points-of-interest (POI) and  includes anonymized counts of how many people visit these POI each month. The counts of visitors are derived from a panel of anonymized mobile device data (see **What is the Panel Sample?** below). These counts are broken down by different dimensions, including by the **home census block group of the visitor**. For example, the dataset may report that 85 total devices visited a specific POI. It also reports that 50 of those visitors live in census block group X, 25 live in census block Y, and 10 live in census block group Z. \n",
    "\n",
    "\n",
    "\n",
    "The aggregated form of SafeGraph Patterns helps to ensure the protection of individuals' privacy, while also providing actionable data for statistical analysis and data science. For all the details on SafeGraph Patterns, see the [SafeGraph Patterns docs](https://docs.safegraph.com/docs/places-schema#section-patterns/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias). \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WUWFYLLXowUJ"
   },
   "source": [
    "### What is the Panel Sample?\n",
    "\n",
    "But what is the sample from which SafeGraph Patterns is aggregated? SafeGraph Patterns is aggregated from a [panel](https://en.wikipedia.org/wiki/Panel_data) (sample of population that is measured longitudinally) of ~40MM mobile devices (e.g., smart phones) in the USA. \n",
    "\n",
    "Protecting individual consumer privacy is at the core of the SafeGraph mission:  \n",
    "\n",
    "\n",
    "\n",
    "> \"SafeGraph’s mission is to make the world’s data open for innovation while protecting individuals privacy.\" - [SafeGraph Vision and Values](https://blog.safegraph.com/safegraph-vision-and-values-circa-2017-d188d4185df2/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias)\n",
    "\n",
    "The devices in the panel are fully anonymized; **no device-level demographic data exists for devices in the panel**. However, for every device in the panel, SafeGraph *estimates* the home census block group (CBG) of that device, based on the most common location during night time hours. Importantly, none of the anonymized device-specific data is available to us as customers of SafeGraph data; SafeGraph Patterns is an aggregated view. SafeGraph provides an aggregated summary of the panel, [Home Location Distributions by State/Census Block Group](https://docs.safegraph.com/docs/places-schema#section-home-location-distributions-by-state-census-block-group/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias), as part of the Panel Overview Data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OakYTjSBom-z"
   },
   "source": [
    "\n",
    "## We can quantify Sampling Bias by comparing  the Panel Sample to the US Census Population\n",
    "\n",
    "Since the US Census publishes population and demographic data for every census block group (CBG), we can  use the US Census data as a groundtruth against which to benchmark sampling bias in the panel. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5PHoLRDdXk2H"
   },
   "source": [
    "\n",
    "We will explore sampling bias on several dimensions: \n",
    "\n",
    "* Geographic\n",
    "  * State \n",
    "  * County\n",
    "  * Census Block Group\n",
    "  \n",
    "* Demographic (inferred from census block group)\n",
    "  * Race\n",
    "  * Educational Attainment\n",
    "  * Household Income\n",
    " \n",
    "<br/><br/>\n",
    "\n",
    "To start, let's load up the datasets. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "t86CDrDHfi7a"
   },
   "source": [
    "## Notebook and Data Setup\n",
    "\n",
    "**Quick Start -- How to use this notebook**\n",
    "* This is an interactive [Google CoLab notebook](https://colab.research.google.com/notebooks/welcome.ipynb). \n",
    "* If you want, you can simply read through and see my code and my results below. Enjoy!\n",
    "* Moreover, you are viewing an ephemeral copy of this notebook **that you can edit and re-run yourself**. You are highly encouraged to interact with the code and try things--only you can see your ephemeral copy (try double-clicking into this cell right now!)\n",
    "  * If you want to save your changes, then you save a copy to your own Google Drive.\n",
    "* As-is, this notebook reads data from a public Google Drive owned by SafeGraph.\n",
    "  * For this to work, Google will prompt you to login, and confirm that you want to run a notebook written by SafeGraph (you should click \"Run Anyway\"), and authenticate your user. \n",
    "  * Be assured, this notebook does not read any data from your personal drive and your data is not shared with SafeGraph. That's just how Google CoLab reads from Google Drive. \n",
    "\n",
    "Pro-Tip for Power Users:\n",
    "* Alternatively, you can upload your own copy of these datasets into *your* Google Drive and [mount your Google drive to this notebook (see \"Bonus Method\")](hhttps://towardsdatascience.com/3-ways-to-load-csv-files-into-colab-7c14fcbdcb92). This is how I worked with the data during the development of the notebook. It slightly changes how you read the data, so I have left those commands inline below, but commented out with the tag `#ForMountedDrive` in case you prefer this option.\n",
    "  * All of the datasets used in this notebook can be obtained for free from [SafeGraph.com](https://shop.safegraph.com/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias) with the coupon code at the bottom of the notebook.\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "If you want to edit the code interactively, make sure you can run the cell directly below. Otherwise, just read on!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "OqkjFI5pYwHm",
    "outputId": "8d7d0be8-db18-4b7b-daa3-222d53c8128e"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "If your name is A, then you've successfuly run this cell. \n",
      "You are authenticated and can edit and re-run anything in the notebook. Enjoy!\n"
     ]
    }
   ],
   "source": [
    "your_name = \"A\" # << You can edit this\n",
    "\n",
    "################################\n",
    "\n",
    "# These libraries are required to read directly from SafeGraph's public GoogleDrive\n",
    "# This is not necessary if you mount your own Drive directly. \n",
    "\n",
    "from pydrive.auth import GoogleAuth\n",
    "from pydrive.drive import GoogleDrive\n",
    "from google.colab import auth\n",
    "from oauth2client.client import GoogleCredentials\n",
    "\n",
    "# Authenticate and create the PyDrive client.\n",
    "auth.authenticate_user()\n",
    "gauth = GoogleAuth()\n",
    "gauth.credentials = GoogleCredentials.get_application_default()\n",
    "drive = GoogleDrive(gauth)\n",
    "\n",
    "################################\n",
    "\n",
    "print(\"If your name is {0}, then you've successfuly run this cell. \\nYou are authenticated and can edit and re-run anything in the notebook. Enjoy!\".format(your_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WC76yxp30Mxl"
   },
   "outputs": [],
   "source": [
    "# #ForMountedDrive -- mount your google drive with the below command\n",
    "# from google.colab import drive\n",
    "# drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lNGAVwMbm_lo"
   },
   "source": [
    "**Load libraries:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "p2PqcEYpVnXp",
    "outputId": "461db318-76c1-478f-c998-02bedf3fa23b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook is running python version: 3.7.12\n",
      "This notebook is running pandas version: 1.1.5\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import os\n",
    "\n",
    "\n",
    "#~~~~~~ Display Settings\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [8, 8]\n",
    "pd.set_option('display.max_rows', 5000)\n",
    "pd.set_option('display.max_columns', 500)\n",
    "pd.set_option('display.width', 10000)\n",
    "\n",
    "#~~~~~~\n",
    "\n",
    "import platform\n",
    "print(\"This notebook is running python version: {0}\".format(platform.python_version()))\n",
    "print(\"This notebook is running pandas version: {0}\".format(pd.__version__))\n",
    "\n",
    "# This is what the original author is running:\n",
    "# This notebook is running python version: 3.6.8\n",
    "# This notebook is running pandas version: 0.24.2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "JNr78OyMHpn_"
   },
   "outputs": [],
   "source": [
    "################################\n",
    "\n",
    "# These functions are required to read directly from SafeGraph's Public Drive\n",
    "# These are not necessary if you've directly mounted your own google drive\n",
    "\n",
    "def pd_read_csv_drive(id, drive=drive, dtype=None):\n",
    "  downloaded = drive.CreateFile({'id':id}) \n",
    "  downloaded.GetContentFile('Filename.csv')  \n",
    "  return(pd.read_csv('Filename.csv',dtype=dtype))\n",
    "\n",
    "\n",
    "# Note: OpenCensusData public GDrive folder: https://drive.google.com/open?id=1btSS6zo7_wJCCXAigkbhnaoeU-Voa9pG\n",
    "\n",
    "drive_ids = {'cbg_b19.csv' : '1d9GscpWbrnP2xNLKKlgd6xLcFTzJydY4',\n",
    "              'cbg_b01.csv' : '1QqttoDRoKpZM2TyyRwJ8B9c5bYZrHysB',\n",
    "              'cbg_b02.csv' : '1Zqqf3iLDkDWPl2theLlUm_cAbvznj-Kx',\n",
    "              'cbg_b15.csv' : '1xeSZShcX3egZFsalGOFD6Ze2jTof6ri-',\n",
    "              'cbg_field_descriptions.csv' : '1a7_7WxY6eaUIObkVwfknl9C7nPltYxPd',\n",
    "              'cbg_fips_codes.csv' : '1dB_HeAw11TmsZ8MATMedC9j2csTRAiVm',\n",
    "             'home_panel_summary.csv': '1aiwhO6Pw1ZfUOoqUf6mS70s9tJgsAVp5'\n",
    "}\n",
    "\n",
    "################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VHjCUWtF6SaG"
   },
   "source": [
    "### Load Open Census Data\n",
    "\n",
    "Census data can be be difficult to navigate, so SafeGraph [packaged all of it in a convenient single download called Open Census Data](https://blog.safegraph.com/download-open-census-data-visualize-neighborhood-insights-9348145bb5f1/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias). \n",
    "\n",
    "This notebook is pulling Open Census Data from a public Google Drive created by SafeGraph (See Notebook and Data Setup, above).  \n",
    "\n",
    "If you want your own copy of Open Census Data, you can [download Open Census Data here](https://www.safegraph.com/open-census-data/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias). See the [SafeGraph Open Census Data docs](https://docs.safegraph.com/docs/open-census-data/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias) for documentation on what is included and how the files are organized.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "faR7Yxi31FyR"
   },
   "outputs": [],
   "source": [
    "# #ForMountedDrive --  After Mounting Drive we can access files in your google drive\n",
    "# sg_ocd_dir = \"/content/drive/My Drive/publicly_shared/openCensusData/\" # copy/paste path for SafeGraph Open Census Data\n",
    "# print(\"expected output is ['metadata', 'geometry', 'data']\")\n",
    "# print(os.listdir(sg_ocd_dir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aDwGk1sS8dgQ"
   },
   "source": [
    "The total population of each `census_block_group` is stored in the file `cbg_b01.csv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "id": "id9Y2g6A3QQO",
    "outputId": "33c3c3e6-e8c9-4704-a1e1-3a072e261b71"
   },
   "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>census_block_group</th>\n",
       "      <th>B01001e1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>010010201001</td>\n",
       "      <td>745</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>010010201002</td>\n",
       "      <td>1265</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>010010202001</td>\n",
       "      <td>960</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>010010202002</td>\n",
       "      <td>1236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>010010203001</td>\n",
       "      <td>2364</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  census_block_group  B01001e1\n",
       "0       010010201001       745\n",
       "1       010010201002      1265\n",
       "2       010010202001       960\n",
       "3       010010202002      1236\n",
       "4       010010203001      2364"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# cbg_b01 = pd.read_csv(sg_ocd_dir + \"data/cbg_b01.csv\",dtype = {'census_block_group': str}) #ForMountedDrive\n",
    "cbg_b01 = pd_read_csv_drive(drive_ids['cbg_b01.csv'],dtype = {'census_block_group': str})\n",
    "columns = ['census_block_group', 'B01001e1']\n",
    "cbg_b01 = cbg_b01[columns]\n",
    "cbg_b01.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "WUDg_XLy3qY9",
    "outputId": "9c2794b5-812c-4b43-9214-63a11f8aa0a8"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "B01001e1\n",
      "4    SEX BY AGE: Total: Total population -- (Estimate)\n",
      "Name: field_full_name, dtype: object\n"
     ]
    }
   ],
   "source": [
    "# cbg_field_desc = pd.read_csv(sg_ocd_dir + \"metadata/cbg_field_descriptions.csv\") #ForMountedDrive\n",
    "cbg_field_desc = pd_read_csv_drive(drive_ids['cbg_field_descriptions.csv'])\n",
    "print('B01001e1')\n",
    "print(cbg_field_desc[cbg_field_desc.table_id == 'B01001e1'].field_full_name)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hoIIvgHv9Fif"
   },
   "source": [
    "So, `B01001e1` is the `table_id` code for the  `Total Population`. This dataset tells us the Total Population for each census block group. \n",
    "\n",
    "For the sake of clarity, we rename this column.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "1rBmxern3jGq",
    "outputId": "1d96103c-75d6-40bc-91c6-b8b71bb81469"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total USA population according to US Census: 322087547\n"
     ]
    }
   ],
   "source": [
    "cbg_b01.rename({'B01001e1':'census_population'}, axis='columns', inplace=True)\n",
    "\n",
    "total_usa_pop = cbg_b01['census_population'].sum()\n",
    "print(\"Total USA population according to US Census: {0}\".format(total_usa_pop))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IJW1dfsK_RZf"
   },
   "source": [
    "### Load SafeGraph Panel Summary\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TBqeRr3GTMza"
   },
   "outputs": [],
   "source": [
    "# #ForMountedDrive - Upload your patterns data to Google drive and copy the path\n",
    "# patterns_path = \"/content/drive/My Drive/publicly_shared/Major-Brands-2019-Oct-CORE-PATTERNS\"\n",
    "# # Confirm correct structure\n",
    "# print(\"expected files: ['home_panel_summary.csv','core_poi-patterns.csv','brand_info.csv','visit_panel_summary.csv']\")\n",
    "# os.listdir(patterns_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3yF9cR_KACSm"
   },
   "source": [
    "For reference, see [schema of home_panel_summary.csv ](https://docs.safegraph.com/docs/places-schema#section-home-location-distributions-by-state-census-block-group/?utm_source=content&utm_medium=referral&utm_campaign=colabnotebook&utm_content=panel_bias)\n",
    "\n",
    "This dataset tells us how many devices are residing in each CBG in the SafeGraph dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "id": "8wsO2TG-ABbx",
    "outputId": "60c6bcab-c08f-4f4e-b2ec-510291cd9522"
   },
   "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>year</th>\n",
       "      <th>month</th>\n",
       "      <th>state</th>\n",
       "      <th>census_block_group</th>\n",
       "      <th>number_devices_residing</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2019</td>\n",
       "      <td>9</td>\n",
       "      <td>ak</td>\n",
       "      <td>021580001001</td>\n",
       "      <td>31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2019</td>\n",
       "      <td>9</td>\n",
       "      <td>ak</td>\n",
       "      <td>020200023033</td>\n",
       "      <td>114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2019</td>\n",
       "      <td>9</td>\n",
       "      <td>ak</td>\n",
       "      <td>021050003001</td>\n",
       "      <td>105</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2019</td>\n",
       "      <td>9</td>\n",
       "      <td>ak</td>\n",
       "      <td>021700006042</td>\n",
       "      <td>67</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2019</td>\n",
       "      <td>9</td>\n",
       "      <td>ak</td>\n",
       "      <td>021880001003</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   year  month state census_block_group  number_devices_residing\n",
       "0  2019      9    ak       021580001001                       31\n",
       "1  2019      9    ak       020200023033                      114\n",
       "2  2019      9    ak       021050003001                      105\n",
       "3  2019      9    ak       021700006042                       67\n",
       "4  2019      9    ak       021880001003                        8"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# home_panel = pd.read_csv(patterns_path + '/home_panel_summary.csv', dtype = {'census_block_group': str}) #ForMountedDrive\n",
    "home_panel = pd_read_csv_drive(drive_ids['home_panel_summary.csv'], dtype = {'census_block_group': str})\n",
    "home_panel.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IiDR7_7kA2D0"
   },
   "source": [
    "I want to compute the correlation between number of devices residing and census populations at the census __tract__ level in specific tracts in MA.\n",
    "\n",
    "counties = ['middlesex', 'norfolk', 'suffolk']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "H7IT9SfbIE74"
   },
   "outputs": [],
   "source": [
    "counties = ['middlesex', 'norfolk', 'suffolk']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "id": "vMA3FjQqJkfF",
    "outputId": "a374ef8e-fe0e-4338-cf22-111b72df76bd"
   },
   "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>census_tract</th>\n",
       "      <th>census_population</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>01001020100</td>\n",
       "      <td>2010</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>01001020200</td>\n",
       "      <td>2196</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>01001020300</td>\n",
       "      <td>3136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>01001020400</td>\n",
       "      <td>4563</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>01001020500</td>\n",
       "      <td>10529</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  census_tract  census_population\n",
       "0  01001020100               2010\n",
       "1  01001020200               2196\n",
       "2  01001020300               3136\n",
       "3  01001020400               4563\n",
       "4  01001020500              10529"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cbg_b01['census_tract'] = [fips[0:11] for fips in cbg_b01.census_block_group]\n",
    "ct = cbg_b01.groupby([\"census_tract\"]).sum().reset_index()\n",
    "ct.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "id": "LB_T5rIEK_ip",
    "outputId": "d70df11a-80d9-4f89-e2b3-572a87dcbd47"
   },
   "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>census_tract</th>\n",
       "      <th>number_devices_residing</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>25001010100</td>\n",
       "      <td>1843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>25001010206</td>\n",
       "      <td>1339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>25001010208</td>\n",
       "      <td>392</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>25001010304</td>\n",
       "      <td>400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>25001010306</td>\n",
       "      <td>244</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  census_tract  number_devices_residing\n",
       "0  25001010100                     1843\n",
       "1  25001010206                     1339\n",
       "2  25001010208                      392\n",
       "3  25001010304                      400\n",
       "4  25001010306                      244"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "home_panel['census_tract'] = [fips[0:11] for fips in home_panel.census_block_group]\n",
    "panel = home_panel[(home_panel.state==\"ma\")].groupby([\"census_tract\"]).sum()[[\"number_devices_residing\"]].reset_index()\n",
    "panel.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EEPQVGlfGtdH"
   },
   "source": [
    "### Join Census and Home Panel data to compare relative population distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "id": "MOLhli-kLZiz",
    "outputId": "cdca0144-8343-4c61-b8e5-bdea903015aa"
   },
   "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>census_tract</th>\n",
       "      <th>number_devices_residing</th>\n",
       "      <th>census_population</th>\n",
       "      <th>state_fips</th>\n",
       "      <th>county_fips</th>\n",
       "      <th>state</th>\n",
       "      <th>county</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>25001010100</td>\n",
       "      <td>1843</td>\n",
       "      <td>2962</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>25001010206</td>\n",
       "      <td>1339</td>\n",
       "      <td>3168</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>25001010208</td>\n",
       "      <td>392</td>\n",
       "      <td>1589</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>25001010304</td>\n",
       "      <td>400</td>\n",
       "      <td>2107</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>25001010306</td>\n",
       "      <td>244</td>\n",
       "      <td>2817</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  census_tract  number_devices_residing  census_population state_fips county_fips state             county\n",
       "0  25001010100                     1843               2962         25         001    MA  Barnstable County\n",
       "1  25001010206                     1339               3168         25         001    MA  Barnstable County\n",
       "2  25001010208                      392               1589         25         001    MA  Barnstable County\n",
       "3  25001010304                      400               2107         25         001    MA  Barnstable County\n",
       "4  25001010306                      244               2817         25         001    MA  Barnstable County"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "panel = pd.merge(panel,ct)\n",
    "panel['state_fips'] = [fips[0:2] for fips in panel.census_tract]\n",
    "panel['county_fips'] = [fips[2:5] for fips in panel.census_tract]\n",
    "# Also load in FIPS codes for county names\n",
    "# fips_codes = pd.read_csv(sg_ocd_dir + \"metadata/cbg_fips_codes.csv\", dtype = {'state_fips': str, 'county_fips':str}) #ForMountedDrive\n",
    "fips_codes = pd_read_csv_drive(drive_ids['cbg_fips_codes.csv'], dtype = {'state_fips': str, 'county_fips':str})\n",
    "county_codes = fips_codes[['state_fips','county_fips', 'state', 'county']].drop_duplicates()\n",
    "panel = pd.merge(panel, county_codes)\n",
    "panel.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 222
    },
    "id": "wynazbO1_jBW",
    "outputId": "326e8b2f-8768-4fd4-b334-1ad018ee64d1"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "counties in MA: 14\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>census_tract</th>\n",
       "      <th>number_devices_residing</th>\n",
       "      <th>census_population</th>\n",
       "      <th>state_fips</th>\n",
       "      <th>county_fips</th>\n",
       "      <th>state</th>\n",
       "      <th>county</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>25001010100</td>\n",
       "      <td>1843</td>\n",
       "      <td>2962</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>25001010206</td>\n",
       "      <td>1339</td>\n",
       "      <td>3168</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>25001010208</td>\n",
       "      <td>392</td>\n",
       "      <td>1589</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>25001010304</td>\n",
       "      <td>400</td>\n",
       "      <td>2107</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>25001010306</td>\n",
       "      <td>244</td>\n",
       "      <td>2817</td>\n",
       "      <td>25</td>\n",
       "      <td>001</td>\n",
       "      <td>MA</td>\n",
       "      <td>Barnstable County</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  census_tract  number_devices_residing  census_population state_fips county_fips state             county\n",
       "0  25001010100                     1843               2962         25         001    MA  Barnstable County\n",
       "1  25001010206                     1339               3168         25         001    MA  Barnstable County\n",
       "2  25001010208                      392               1589         25         001    MA  Barnstable County\n",
       "3  25001010304                      400               2107         25         001    MA  Barnstable County\n",
       "4  25001010306                      244               2817         25         001    MA  Barnstable County"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('counties in MA: %s' % len(fips_codes[fips_codes.state == 'MA']))\n",
    "panel[panel.state == 'MA'].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YS5Hl_kENa5D"
   },
   "source": [
    "Middlesex county (017)\n",
    "\n",
    "Norfolk county (021)\n",
    "\n",
    "Suffolk county (025)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 258
    },
    "id": "13ETs8kq72kQ",
    "outputId": "8a9291e9-0ef8-492d-84f8-4e4a75b83cd0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total population for counties (middlesex, norfolk, suffolk): 3026547\n",
      "total devices for counties (middlesex, norfolk, suffolk): 396061\n",
      "13.09 %\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>census_tract</th>\n",
       "      <th>number_devices_residing</th>\n",
       "      <th>census_population</th>\n",
       "      <th>state_fips</th>\n",
       "      <th>county_fips</th>\n",
       "      <th>state</th>\n",
       "      <th>county</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>547</th>\n",
       "      <td>25017300100</td>\n",
       "      <td>235</td>\n",
       "      <td>3196</td>\n",
       "      <td>25</td>\n",
       "      <td>017</td>\n",
       "      <td>MA</td>\n",
       "      <td>Middlesex County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>548</th>\n",
       "      <td>25017301101</td>\n",
       "      <td>300</td>\n",
       "      <td>4092</td>\n",
       "      <td>25</td>\n",
       "      <td>017</td>\n",
       "      <td>MA</td>\n",
       "      <td>Middlesex County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>549</th>\n",
       "      <td>25017301102</td>\n",
       "      <td>361</td>\n",
       "      <td>5277</td>\n",
       "      <td>25</td>\n",
       "      <td>017</td>\n",
       "      <td>MA</td>\n",
       "      <td>Middlesex County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>550</th>\n",
       "      <td>25017310100</td>\n",
       "      <td>904</td>\n",
       "      <td>5962</td>\n",
       "      <td>25</td>\n",
       "      <td>017</td>\n",
       "      <td>MA</td>\n",
       "      <td>Middlesex County</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>551</th>\n",
       "      <td>25017310200</td>\n",
       "      <td>533</td>\n",
       "      <td>6735</td>\n",
       "      <td>25</td>\n",
       "      <td>017</td>\n",
       "      <td>MA</td>\n",
       "      <td>Middlesex County</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    census_tract  number_devices_residing  census_population state_fips county_fips state            county\n",
       "547  25017300100                      235               3196         25         017    MA  Middlesex County\n",
       "548  25017301101                      300               4092         25         017    MA  Middlesex County\n",
       "549  25017301102                      361               5277         25         017    MA  Middlesex County\n",
       "550  25017310100                      904               5962         25         017    MA  Middlesex County\n",
       "551  25017310200                      533               6735         25         017    MA  Middlesex County"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "panel = panel[(panel.county_fips=='017') | (panel.county_fips=='021') | (panel.county_fips=='025')]\n",
    "print('total population for counties (middlesex, norfolk, suffolk): %s' % panel.census_population.sum())\n",
    "print('total devices for counties (middlesex, norfolk, suffolk): %s' % panel.number_devices_residing.sum())\n",
    "print('%0.2f %%' % (100*panel.number_devices_residing.sum()/panel.census_population.sum()))\n",
    "panel.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zq5uJVTaHgqU"
   },
   "source": [
    "Check if we are missing any CBGs in the Home Panel dataset, and set the value to 0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "id": "pBAbiovm0NBc",
    "outputId": "bc1fb7da-c67a-4176-9858-7287ac64a407"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We are missing SG data for 0 tracts containing 0 residents\n",
      "Set number_devices_residing to 0 for missing CBGs\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>census_tract</th>\n",
       "      <th>number_devices_residing</th>\n",
       "      <th>census_population</th>\n",
       "      <th>state_fips</th>\n",
       "      <th>county_fips</th>\n",
       "      <th>state</th>\n",
       "      <th>county</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [census_tract, number_devices_residing, census_population, state_fips, county_fips, state, county]\n",
       "Index: []"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "missing_cts = panel[panel.number_devices_residing.isnull()].sort_values('census_population', ascending = False)\n",
    "print(\"We are missing SG data for {0} tracts containing {1} residents\".format(missing_cts.shape[0],missing_cts.census_population.sum()))\n",
    "panel = panel.fillna(0)\n",
    "print(\"Set number_devices_residing to 0 for missing CBGs\") # For these, we will set the value to 0, since the panel has 0 devices in these CBGs. \n",
    "panel[(panel.number_devices_residing == 0) & (panel.census_population > 0)].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "qUAHU3UBIj1S",
    "outputId": "537b110a-d42f-4396-d810-b64c45d939d7"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "census_population          3026547\n",
       "number_devices_residing     396061\n",
       "dtype: int64"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_population = panel[['census_population','number_devices_residing']].sum(axis = 0)\n",
    "total_population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nUrAZqaTxzU8"
   },
   "source": [
    "Briefly, we define some functions given the above framework to aid in exploratory data analysis and visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NPPbHRkeLmwN"
   },
   "outputs": [],
   "source": [
    "# functions for exploratory data analysis and visualization of Sampling Bias\n",
    "\n",
    "\n",
    "def compute_bias(df, sample_var_col='number_devices_residing', pop_var_col='census_population'):\n",
    "  # compute bias for a dataframe of observed and expected values\n",
    "  # returns a pd.Series\n",
    "  observed_proportion = df[sample_var_col] / df[sample_var_col].sum()*100\n",
    "  expected_proportion = df[pop_var_col] / df[pop_var_col].sum()*100\n",
    "  bias = observed_proportion - expected_proportion\n",
    "  return(bias)\n",
    "  \n",
    "\n",
    "def make_center_line_coords(df,x,y):\n",
    "  # This creates two coords to define the hypothetical line in which x and y are perfectly correlatd\n",
    "  expected_slope = df[y].sum()/df[x].sum()\n",
    "  xcoord = np.array([df[x].min(), df[x].max()])\n",
    "  ycoord = expected_slope*xcoord\n",
    "  return((xcoord,ycoord))\n",
    "  \n",
    "def plot_sg_vs_census(df, xvar, yvar,bias_col='bias',xlim_=None,ylim_=None,text_col=None, num_to_label=5,alpha=0.9,labelfontsize=14):\n",
    "  # `text_col` is the column used to annotate the points\n",
    "  # `num_to_label` sets the top num_to_label records sorted by yvar to be annotated\n",
    "  \n",
    "  plt.rcParams['figure.figsize'] = [7, 7]\n",
    "  \n",
    "  # Report Stats\n",
    "  pearson_r = df.corr().loc[xvar, yvar]\n",
    "  print(\"pearson r: {0:.3f}\".format(pearson_r))\n",
    "  \n",
    "  # scatter plot\n",
    "  ax = df.plot.scatter(xvar, yvar, alpha = alpha, edgecolor = 'k', xlim = xlim_, ylim = ylim_)\n",
    "  # annotate \n",
    "  xadj = df[xvar].max()*0.02\n",
    "  if(text_col is not None):\n",
    "    for i, point in df.sort_values(by=xvar,ascending=False).iloc[0:num_to_label].iterrows():\n",
    "      ax.text(point[xvar] + xadj, point[yvar], str(point[text_col]),fontsize=labelfontsize, verticalalignment='center')\n",
    "  # plot expectation line\n",
    "  line_coords = make_center_line_coords(df,xvar,yvar)\n",
    "  plt.plot(line_coords[0],line_coords[1], color = 'k',linestyle='-', linewidth=1)\n",
    "  return(ax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9V0gb9B8QjZg",
    "outputId": "45f8be07-0d29-4e10-d867-cba79f005b7e"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "651"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(panel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xz0kAF1ePjjL"
   },
   "source": [
    "## Sampling Bias at the Tract Level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 479
    },
    "id": "iOBSaYSYQ-4L",
    "outputId": "a0a85740-7f10-48dd-e780-0436062f21f9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Figure: Tract-Level Sampling Bias\n",
      "pearson r: 0.122\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGqCAYAAAAWZQ0OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxU9b3/8dcnyWSbBJIQtoQdkUWDGyKgrUt73a1LrVbr3p9aa+u9rbcu11at9tbautVeFeh1AddyUat1KaIi1gUUXIIsyqIsYU8ISSYJmSTf3x9zAkMIkJnMZJm8n49HHpk5c2bmO5nJvM93PeacQ0RERBJHUkcXQERERGJL4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCSalowvQEfLz892QIUM6uhgiItJJLFy4cKtzrndHlyNWumW4DxkyhAULFnR0MUREpJMws9UdXYZYUrO8iIhIglG4i4iIJBiFu4iISIJRuIuIiCQYhbuIiEiCUbiLiIgkGIW7iIhIglG4i4iIJBiFu4iISIJRuIuIiCQYhbuIiEiCUbiLiIgkGIW7iIhIglG4i0jMBQIBSkpKCAQCHV0UkW6pW57yVUTip7h4EVOfnkmjL5OkYDVXX3QuRUVFHV0skW5FNXcRiZlAIMDUp2fS54iTGHHs2fQ54iSmPDVTNXiRdqZwF5GYKS8vp9GXSXZuPgDZufk0+jIpLy/v4JKJdC8KdxGJmZycHJKC1VRu2wpA5batJAWrycnJ6eCSiXQv6nMXkZjx+/1cfdG5THlqJhvD+tz9fn9HF02kW1G4i0hMFRUVcfetwygvLycnJ0fBLtIBFO4iEnN+v1+hLtKB1OcuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4S4iIpJg4hruZvaYmW02sy/Ctv3JzJaZWbGZvWhmOWG33WxmK8zsSzM7KWz7yd62FWZ2U9j2oWY239v+NzNLjefrERER6QriXXN/Aji52bbZwMHOubHAV8DNAGY2BvghcJB3n4fNLNnMkoGHgFOAMcAF3r4AdwP3O+cOALYBP47vyxEREen84hruzrl3gbJm295wztV7V+cBA7zLZwLPOed2OOe+BlYA472fFc65Vc65OuA54EwzM+AEYKZ3/2nAWfF8PSIiIl1BR/e5XwG87l0uBNaG3bbO27a37b2A8rADhabtLTKzq8xsgZkt2LJlS4yKLyIi0vl0WLib2S1APfB0ezyfc26qc26cc25c79692+MpRUREOkRKRzypmV0GnA58xznnvM0lwMCw3QZ429jL9lIgx8xSvNp7+P4iIiLdVrvX3M3sZOAG4HvOueqwm14GfmhmaWY2FBgBfAR8DIzwRsanEhp097J3UDAHONe7/6XAS+31OkRERDqreE+Fexb4EBhpZuvM7MfA/wDZwGwz+8zMJgM45xYDM4AlwD+Ba51zDV6t/GfALGApMMPbF+BG4JdmtoJQH/yj8Xw9IiIiXYHtahXvPsaNG+cWLFjQ0cUQEZFOwswWOufGdXQ5YqWjR8uLiIhIjCncRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBKMwl1ERCTBKNxFREQSjMJdREQkwSjcRUREEozCXUREJMEo3EVERBJMXMPdzB4zs81m9kXYtjwzm21my73fud52M7MHzWyFmRWb2eFh97nU23+5mV0atv0IM1vk3edBM7N4vh4REZGuIN419yeAk5ttuwl4yzk3AnjLuw5wCjDC+7kKeARCBwPAbcBRwHjgtqYDAm+fK8Pu1/y5REREup24hrtz7l2grNnmM4Fp3uVpwFlh26e7kHlAjpn1B04CZjvnypxz24DZwMnebT2cc/Occw6YHvZYIiIi3VZH9Ln3dc5t8C5vBPp6lwuBtWH7rfO27Wv7uha2t8jMrjKzBWa2YMuWLW17BSIiIp1Yhw6o82rcrp2ea6pzbpxzblzv3r3b4ylFREQ6REeE+yavSR3v92ZvewkwMGy/Ad62fW0f0MJ2ERGRbq0jwv1loGnE+6XAS2HbL/FGzU8AtnvN97OAE80s1xtIdyIwy7utwswmeKPkLwl7LBERkW4rJZ4PbmbPAscB+Wa2jtCo9z8AM8zsx8Bq4Dxv99eAU4EVQDVwOYBzrszM7gQ+9va7wznXNEjvp4RG5GcAr3s/IiIi3ZqFur27l3HjxrkFCxZ0dDFERKSTMLOFzrlxHV2OWNEKdSIiIglG4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4S4iIpJgFO4iIiIJRuEuIiKSYBTuIiIiCUbhLiIikmAU7iIiIglG4R6FQCBASUkJgUCgo4siIiKyh5TW7mhmh7eweTuw2jlXH7sidW7FxYuY+vRMGn2ZJAWrufqicykqKuroYomIiOzU6nAHHgYOB4oBAw4GFgM9zewa59wbcShfpxIIBJj69Ez6HHES2bn5VG7bypSnZnL3rcPw+/0dXTwREREgsmb59cBhzrlxzrkjgMOAVcC/AX+MR+E6m/Lychp9mWTn5gOQnZtPoy+T8vLyDi6ZiIjILpGE+4HOucVNV5xzS4BRzrlVsS9W55STk0NSsJrKbVsBqNy2laRgNTk5OR1cMhERkV0iaZZfbGaPAM95188HlphZGhCMeck6Ib/fz9UXncuUp2ayMazPXU3yIiLSmZhzrnU7mmUAPwWO8Ta9T6gfvhbIdM5VxaWEcTBu3Di3YMGCqO8fCAQoLy8nJydHwS4ikgDMbKFzblxHlyNWWl1zd87VAPd6P811mWCPBb/fr1AXEZFOK5KpcEcDtwODw+/nnBsW+2KJiIhItCLpc38U+AWwEGiIT3FERESkrSIJ9+3OudfjVhIRERGJiUjCfY6Z/Ql4AdjRtNE590nMSyUiIiJRiyTcj/J+h48mdMAJsSuOiIiItFUko+WPj2dBREREJDb2G+5mdpFz7ikz+2VLtzvn7ovmic3sF8D/I1T7XwRcDvQntEhOL0ID9y52ztV5C+VMB44ASoHznXPfeI9zM/BjQoP8rnPOzYqmPCIiIomiNcvPNk3ozt7LT8TMrBC4DhjnnDsYSAZ+CNwN3O+cOwDYRii08X5v87bf7+2HmY3x7ncQcDLwsJklR1MmERGRRLHfmrtzbor3+7dxeO4MMwsCmcAGQv33F3q3TyM0r/4R4EzvMsBM4H/MzLztzznndgBfm9kKYDzwYYzLKiIi0mW0pln+wX3d7py7LtIndc6VmNk9wBqgBniDUDN8edi54dcBhd7lQmCtd996M9tOqOm+EJgX9tDh92n+Oq4CrgIYNGhQpEUWERHpMlrTLL/Q+0kndD735d7PoUBqNE9qZrmEat1DgQJCTf8nR/NYreWcm+qdrnZc79694/lUIiIiHao1zfLTAMzsGuCYppq1mU0G/hXl834X+No5t8V7rBeAo4EcM0vxnmMAUOLtXwIMBNaZWQrQk9DAuqbtTcLvIyIi0i1Fcj73XKBH2PUsb1s01gATzCzT6zv/DrAEmAOc6+1zKfCSd/ll7zre7W+70OnsXgZ+aGZpZjYUGAF8FGWZREREEkIki9j8AfjUzOYABnybXYPcIuKcm29mM4FPgHrgU2Aq8CrwnJn9ztv2qHeXR4EnvQFzZYRGyOOcW2xmMwgdGNQD1zrntO69iIh0a60+nzuAmfVj10p1851zG+NSqjhr6/ncRUQksSTa+dz32yxvZqO834cTGvy21vsp8LaJiIhIJ9KaZvlfEppCdm8Lt2lteRERkU6mNaPlr/J+a215ERGRLqDVo+XN7Admlu1d/rWZvWBmh8WvaCIiIhKNSKbC/cY5V2lmxxCap/4oMDk+xRIREZFoRRLuTVPMTgOmOudeJcoV6kRERCR+Ign3EjObApwPvOadhjWS+4uIiEg7iCSczwNmASc558qBPOBXcSmViIiIRK3V4e6cqwY2A8d4m+oJnUBGREREOpFIRsvfBtwI3Oxt8gFPxaNQIiIiEr1ImuXPBr4HBACcc+uB7HgUSkRERKIXSbjXeWdicwBm5o9PkURERKQtIgn3Gd5o+RwzuxJ4E/hrfIolIiIi0WrVKV+9c67/DRgFVAAjgVudc7PjWDYRERGJQqvC3TnnzOw151wRoEAXERHpxCJplv/EzI6MW0lEREQkJlpVc/ccBfzIzFYTGjFvhCr1Y+NSMhEREYlKJOF+0r5uNLNc59y2NpZHRERE2qjV4e6cW72fXd4CDm9bcURERKStYnniF4vhY4mIiEiUYhnuLoaPJSIiIlHSKVtFREQSjJrlRUREEkwkZ4UbbmZp3uXjzOw6M8sJ2+U7MS+diIiIRCySmvvzQIOZHQBMBQYCzzTd6Jwri3HZREREJAqRhHujc66e0Klf/+Kc+xXQPz7FEhERkWhFEu5BM7sAuBR4xdvmi32RREREpC0iCffLgYnAfzvnvjazocCT8SmWiIiIRCuSFeqWmNmNwCDv+tfA3fEqmIiIiEQnktHyZwCfAf/0rh9qZi/Hq2AiIiISnUia5W8HxgPlAM65z4BhcSiTiIiItEFEA+qcc9ubbWuMZWFERESk7SI55etiM7sQSDazEcB1wAfxKZaIiIhEK5Ka+8+Bg4AdhBav2Q78RzwKJSIiItGLZLR8NXCL9yMiIiKdVCSj5WeHryVvZrlmNis+xRIREZFoRdIsn++cK2+64pzbBvSJfZFERESkLSJaW97MBjVdMbPBgIt9kURERKQtIhktfwvwnpnNJXTu9m8BV8WlVCIiIhK1SAbU/dPMDgcmeJv+wzm3NT7FEhERkWjtt1nezEZ5vw8ntK78eu9nkLdNREREOpHW1Nx/Saj5/d4WbnPACTEtkYiIiLTJfsPdOXeV9/v4+BdHRERE2iqSee7FZnazmQ2PZ4FERESkbSKZCncG0ADMMLOPzew/w6fGiYiISOfQ6nB3zq12zv3ROXcEcCEwFvg62ic2sxwzm2lmy8xsqZlNNLM8byW85d7vXG9fM7MHzWyF14JweNjjXOrtv9zMLo22PCIiIokikpo7ZjbYzG4AngNGATe04bn/DPzTOTcKOARYCtwEvOWcGwG85V0HOAUY4f1cBTzilScPuA04itC55m9rOiAQERHprlo9z93M5gM+4P+AHzjnVkX7pGbWE/g2cBmAc64OqDOzM4HjvN2mAe8ANwJnAtOdcw6Y59X6+3v7znbOlXmPOxs4GXg22rKJiIh0dZGsUHeJc+7LGD3vUGAL8LiZHQIsBP4d6Ouc2+DtsxHo610uBNaG3X+dt21v2/dgZlfhrag3aJCGCoiISOKKpFm+3MweNbPXAcxsjJn9OMrnTQEOBx5xzh0GBNjVBA+AV0uP2dr1zrmpzrlxzrlxvXv3jtXDioiIdDqRhPsTwCygwLv+FfAfUT7vOmCdc26+d30mobDf5DW34/3e7N1eAgwMu/8Ab9vetouIiHRbkZ7ydQbQCOCcqyc0NS5izrmNwFozG+lt+g6wBHgZaBrxfinwknf5ZeASb9T8BGC713w/CzjRO7d8LnCit01ERKTbiqTPPWBmvfCayptCtg3P/XPgaTNLBVYBlxM62JjhNfevBs7z9n0NOBVYAVR7++KcKzOzO4GPvf3uaBpcJyIi0l1ZqGu7FTuG5pb/BTgY+ALoDZzrnCuOX/HiY9y4cW7BggUdXQwREekkzGyhc25cR5cjViI55esnZnYsMJLQ+dy/dM4F41YyERERicp+w93MztnLTQeaGc65F2JcJhEREWmD1tTcz/B+9wEmAW97148HPgAU7iIiIp1Ia075ejmAmb0BjGlaZMabqvZEXEsnIiIiEYtkKtzAsNXjADYBWupNRESkk4lkKtxbZjaLXeu2nw+8GfsiiYiISFtEMlr+Z2Z2NqETvgBMdc69GJ9iiYiISLQiqbnjhXmLgW5mHzrnJsakVCIiIhK1iM7nvh/pMXwsERERiVIswz1mZ3ATERGR6MUy3EVERKQTiGW4WwwfS0RERKLUqnA3s2Qzm7Of3S6OQXlERESkjVoV7s65BqDRzHruY58vYlYqERERiVokU+GqgEVmNhsING10zl0X81KJiIhI1CIJ9xfQSWJEREQ6vUhWqJtmZhnAIOfcl3Esk4iIiLRBq0fLm9kZwGfAP73rh5rZy/EqmIiIiEQnkqlwtwPjgXIA59xnwLA4lElERETaIJJwDzrntjfb1hjLwoiIiEjbRTKgbrGZXQgkm9kI4Drgg/gUS0RERKIVSc3958BBwA5C53SvAP4jHoUSERGR6EUyWr4auMXM7g5ddZXxK5aIiIhEK5LR8kea2SKgmNBiNp+b2RHxK5qIiIhEI5I+90eBnzrn/gVgZscAjwNj41EwERERiU4kfe4NTcEO4Jx7D6iPfZFERESkLfZbczezw72Lc81sCqHBdA44H3gnfkUTERGRaLSmWf7eZtdvC7vsYlgWERERiYH9hrtz7vj2KIiIiIjERqsH1JlZDnAJMCT8fjrlq4iISOcSyWj514B5wCK07KyIiEinFUm4pzvnfhm3koiIiEhMRDIV7kkzu9LM+ptZXtNP3EomIiIiUYmk5l4H/Am4hV2j5B067auIiEinEkm4Xw8c4JzbGq/CiIiISNtF0iy/AqiOV0FEREQkNiKpuQeAz8xsDqHTvgKaCiciItLZRBLuf/d+REREpBOL5Hzu0+JZEBEREYmNSFao+5oW1pJ3zmm0vIiISCcSSbP8uLDL6cAPAM1zFxER6WRaPVreOVca9lPinHsAOC2OZRMREZEoRNIsf3jY1SRCNflIav4iIiLSDiIJ53vZ1edeD3xDqGleRDpQIBCgvLycnJwc/H5/RxdHRDqBSML9FOD77H7K1x8Cd8S4TCLSSsXFi5j69EwafZkkBau5+qJzKSoq6uhiiUgHi2SFur8DZwBBoMr7CcSjUCKyf4FAgKlPz6TPEScx4tiz6XPESUx5aiaBgP4tRbq7SGruA5xzJ8fyyc0sGVgAlDjnTjezocBzQC9gIXCxc67OzNKA6cARQClwvnPuG+8xbgZ+DDQA1znnZsWyjCKdVXl5OY2+TLJz8wHIzs1noy+T8vJyNc+LdHOR1Nw/MLNYt/f9O7A07PrdwP3OuQOAbYRCG+/3Nm/7/d5+mNkYQl0DBwEnAw97BwwiCS8nJ4ekYDWV20LncqrctpWkYDU5OTkdXDIR6WiRhPsxwEIz+9LMis1skZkVR/vEZjaA0FS6//WuG3ACMNPbZRpwlnf5TO863u3f8fY/E3jOObfDOfc1oZPbjI+2TCJdid/v5+qLzmXzwlksn/simxfO4uqLzlWtXUQiHlAXSw8ANwDZ3vVeQLlzrt67vg4o9C4XAmsBnHP1Zrbd278QmBf2mOH32Y2ZXQVcBTBo0KDYvQqRDlRUVMTdtw7TaHkR2U0ka8uvjtWTmtnpwGbn3EIzOy5Wj7svzrmpwFSAcePG7bGMrkhX5ff7FeoispuOWoTmaOB7ZnYqoaVsewB/BnLMLMWrvQ8ASrz9S4CBwDozSwF6EhpY17S9Sfh9REREuqVI+txjxjl3s3NugHNuCKEBcW87534EzAHO9Xa7FHjJu/yydx3v9redc87b/kMzS/NG2o8APmqnlyEiItIpdbblY28EnjOz3wGfAo962x8FnjSzFUAZoQMCnHOLzWwGsITQqnnXOuca2r/YIiIinYeFKsDdy7hx49yCBQs6uhgiIp1ed1ne2MwWOufG7X/PrqGz1dxFRKST0PLGXVeH9LmLiEjnpuWNuzaFu3RpgUCAkpISfeGIxFhLyxs3essbS+enZnnpstRkKBI/4csbZ+fma3njLkbhLl1SeJNh0xfPlKdmcvetwxJ60E+i6S6DtbqipuWNpzw1k41hB9B6n7oGhbt0STojWtc3b958Jk9/Dl9WHmlWr5aXTkjLG3dd6nOXLklnROva5n/0ET+56bd8ud34sqSUxt4HaLBWJ+X3+yksLFSwdzGquUuXpCbDrisQCDB5+gyyx55IwUFHUVtRxhefzmZkYS+1vIjEiMJduiw1GXZN5eXlpPhzSCeTHTXVpPfII0gSwaoytbyIxIjCXbo0nRGt68nJySGNIKMH92Xp6pWU1lRTU7KMa667Te+lSIyoz11E2lVTl0rj6k8YwmYOaFzH1D/ezlFHHdXRRRNJGKq5i0i7U5eKSHwp3EWkQ6hLRSR+1CwvIiKSYBTuIiIdROdGkHhRs7yISAfQuREknlRzFxFpZzqdqsSbwl1EpJ3pdKoSbwp3EZF2pnMjSLypz70NdLpKEYmGzo0g8aZwj5IGw4hIW2ghH4knhXsUwgfDZOfmU7ltK1Oemsndtw7TP6h0CLUidU1ayEfiReEehZYGw2z0BsPoH1Xam1qRRKQ5DaiLggbDSGehKVUi0hLV3KOgwTDSWagVaU/qohBRuEdNg2GkMwhvRWoa/9GdW5HURSESonBvAw2GkY6mVqRdNNBVZBeFu0gXp1akEHVRiOyicBdJAGpFUheFSDiFu4gkBHVRiOyicBeRhKEuCpEQhbskFE2DEnVRiCjcJYFoGpS0hQ4MJZEo3CUhaBqUtIUODCXRaPlZSQgtTYNq9KZBieyLlvCVRKRwl4Sg9f4lWjowlESkZnlJCJoGJdHS/HhJROac6+gytLtx48a5BQsWdHQxJA40KEqisWjRIqY8pT737szMFjrnxnV0OWJF4S4igg4Mu7tEC3c1y4tESWGQWDQ/XhKJwl0kCpo6JSKdmUbLS7cWCAQoKSmJaNqTpk6JSGenmrt0K+FN6StXroqq9q1Ti0ok1H0jHUHhLt1GeFN6Q6CcyqpKRv3bhRGvaNc0dap0w1qSfak0BOs0dUpapO4b6SgKd+kWmi9PW7JyKfNffIZDs3oAkdW+/X4/x48fy51/uReXkYvVbOM3P/9xzGtlqvGFdNW/g5ZElo7UIeFuZgOB6UBfwAFTnXN/NrM84G/AEOAb4Dzn3DYzM+DPwKlANXCZc+4T77EuBX7tPfTvnHPT2vO1SNfQvCk9v2AwLljD1o0l9B84NKKFSwKBAHM+KuaMK2/Al5FFsKaKOR/N4eSTT47Zl7ZqfCFd+e+g7hvpSB01oK4euN45NwaYAFxrZmOAm4C3nHMjgLe86wCnACO8n6uARwC8g4HbgKOA8cBtZpbbni9Euobmy9PW1VYzekghFYvnsnzui2xeOKvVK9qtX7+eipog/uwcsrN7kNenIKbLlWrAXkhX/zt01JLI0QwSlcTTITV359wGYIN3udLMlgKFwJnAcd5u04B3gBu97dNdaMWdeWaWY2b9vX1nO+fKAMxsNnAy8Gy7vRjpElpanvaWX/yUYcOGRdTkW1y8iIenPUvxlytZvv4RjjnxdLJzesX0S1s1vlBALV26lB34uuzfoSOWRO7KLR0SWx3e525mQ4DDgPlAXy/4ATYSaraHUPCvDbvbOm/b3ra39DxXEar1M2jQoNgUXrqUoqIi7r51zzBv7ZdtU02y4KjTOeXwBt7/cD5vPPdXJhxyED+74sKYfWl397XOmwJqh0vh448+JGNwEUNHjO6Sf4e9febiQX38Eq5Dw93MsoDngf9wzlWEutZDnHPOzGK2Nq5zbiowFULLz8bqcaVracsqZOE16mzg9FNPZrGvgv/86cWMGDEipmXsrifBaR5Qvj5DmPPsQ2wtOpTstOSYHkS1l/Za+U4tPhKuw8LdzHyEgv1p59wL3uZNZtbfObfBa3bf7G0vAQaG3X2At62EXc34TdvfiWe5pftqXqOuraogO9UoKCiI+XO1Z42vM2keUP4eOZCUSrCxEUvq8IbGTq27t/jI7jpkQJ03+v1RYKlz7r6wm14GLvUuXwq8FLb9EguZAGz3mu9nASeaWa43kO5Eb5tITDVNx7rk+6exeeGsiAfhRcPv91NYWNhtgh12D6gdNdW898Yr5BYdx5Hfu4yCo07vUgPq2ltTi097fT6lc+uoQ+GjgYuBRWb2mbftv4A/ADPM7MfAauA877bXCE2DW0FoKtzlAM65MjO7E/jY2++OpsF1IrHSfJDSxeecyoABA7pVjbq9hHdJVNQEqaks57izj8LnS8WnZub96q4tPrKnjhot/x5ge7n5Oy3s74Br9/JYjwGPxa50Iru0NEjpyRde4+5bb9AXZ5w0BdT69eu5d/LjZPiSgfabStbeYr1IT1v7+LvqokGyO3ViieyDBinFx/4CxO/3M2LECK697IKEHljY2aaudbbySPQU7iL7oEFKsRdJgCRyM3Nnm7rW2cojbaNwF9mH7jwtLR6iCZD2mkrW3jpbq1BnK4+0jcJdZD8SufbY3hIxQKLto+5srUKtLY/65LsGhbtIKyRq7bG9dbZAa6u29FF3tlah1pRHffJdh4UGoncv48aNcwsWLOjoYoh0S4sWLWLKU9EFRGeqNQYCAW6844+7dTFsXjgr4pkUnek17as8sXq9nZWZLXTOjevocsSKau4dpLP9Q4u0l2i7OfZWa+yo/6VYdTF0tlahvZUnEbtUEpnCvQOoaUu6u0gDbW8D8S46u5qnXnytQ/6XEq2LYX+62+vt6hTu7UzTTToXtaB0DS3VGtfiY/L05xh23A865H+ps/WZx5lgJhMAACAASURBVFt3e71dncK9nalpq/NoqQUl0vO7d2Vd6cCmpVpjfaAcX1ZeVP9LsXrt3W0mRXd7vV2Zwj0KbfliUNNWx2j+nrXUgvLf9z9MTm4uSRk9E767pKt1DbVUa7zmkvOY/vyrLf4v7et/NNavvSP7zDviAK2zjRGQlincI9TWLwY1bbW/lt6zvLy83VpQUtMzWfpNCScdcxb9Bw5N6O6SeHYNxTNsWqo1ZmRk7PG/tHLlqr22yPh8vpi89s7Q6tHVDtCkfSncIxCrL8Wu1rTVGb7IorW39+zW66/drQVl6/rVmC+D/H6FAKRn9aCyzrF+/XpGjBjRwa8ituLVNdQeYdO81tj8fwnYY7pWeItMTfkWqnc0MLwNr70zhKrG7sj+KNwjEMsvxa7StNUZvsjaYm/vWTAY3K0FpSFQzphBvamtqqCsroEP5s1n+5LPuOfhID++4GwKCwu75MFNS+LRNdSRYRP+v1RSUrLPFpmyzev5x1//yIgNa+nVf2DEr709Wz32dVC9v++irnxALrGhcI9Ad+svT4Tawb7es8LCwt1qfatWreKhJ55l3hcrycjO4cQfXkltsJGrbridI8dPJM3qu9zBTUui7RqKNmyabm+PoGn+fjdvkcnrU8CoMUWsm/cKZXn9Iu4Wa69Wj+PHj2XOR8V7Paje1+e6qx+QR2P+/Pn84he/ICkpiSOPPJL777+fnj17cthhhwHwwgsvkJeXx9NPP81DDz1EXl4ezzzzDGbGWWedRTAYBDjAzLKdc5Vm9g6h05I74A7n3NtmdguhU48/5pz7NYCZPQf0A9KADOfcoWZWADwFpAO3OufeNLObgFOATOD3zrkX97LfHtui/Zso3CPQ3frLE2Fk//7es/BaX1FREdf/JJ07//I4Y757HkkpKcye+wEZhaMoOOw4kpOTu9zBzd5E2jXUPDAu+f5pu7Vm7C1s1q5dx10PTm23oGn+foe3yPi8cuVnp3Pr9dcTDAY7xXrwzQ+iSzes5c6/3MsZV95AXp+CFg+q9/a5Brr8AXlL9tcSMXjwYN5++23S09P50Y9+xKJFiygqKuKdd97ZuU8wGGTy5Mm8++67PP/880yZMoWf//znPPXUU/Tv3x8zKwcuA/7i3eU7zrn6sKf5X+AD4DtNG5xzPwQws7OBI7zNNwG/AT4HXgHeBO51zv3BzLK86y/uZb+WtkVF4R6hrtZf3haJ0lIRyXtWUFBAjwwfdbXVkJJKbU01PhrJ8GeTlpHZ5Q5uwjX/gmxt11Dz8Fn11dIWWzMuPudUHpn2HL6sPNKsnku+fxrTn3+13YOm+fu9atWqPUKwT58+UT12PA7wmx9EJ/tScRm5+DKygL0fVLf0uW7eLdEVD8jD1dXV8c47c/nfp2YQ2FFPTUUpRQcOxczYtGkTmzZtorq6mnnz5u28j8/nIzk5maVLl/Ktb32Lo48+mrvuuovly5dTVFRESkoK3/3ud7nyyiv51a9+Rf/+/Zvu6oAG73Ij8KaZbQR+6pwrc85tMrPReynq2cAD3uUi4N+dc87MKs2sh3OuwrstA/hib/vt574RUbhHoav0l7dVV2yp2NsRfmvfs/DXvAMflcUfMOn0C0jLyOyyBzfQtrET4eETDNaxbM2mPVozLv5+DU++8Bq+7F7UB8q58pLzKCwsjDhoYtVX3LxFJpYH5EVFRdx6fV/WrFnDoEGDoj5QaNL8ILohWIfVbCNYUwXZPfb5uWv+ue4KB+S1tbU7g7npZ+PGjXts27RpE4FAgNS0dHrk96dnfl/S/Vn8a94Czjv7DI488kj69u1L3759dz52cXExW7ZsYcyYMSxfvpzc3Fx+8pOf8I9//IP8/Hx69OgBQM+ePXd2GQFUVVUB9Aae8Tad65wrM7MLgV8Dv9zb6zEzH1DknPvE25Tsdp20ZTuQA1SY2cOEDgKu38d+Ld43wj8xoHCX/ehKLRWx6msMf83r1v0bT77wGsvnrm/3g5tYBd3+xk7s73nCA6Ol1oy1LoXJ02cw9Nvn7Hz86c+/useMhP0FTTz7ipteV9MXeiz7xy8+51QGDBgQ9fvU0kH0b37+Y+Z8NIfSCA+q23pAHu1nLhAItBjOLf3U1tbSp0+fncHc9DN8+HAmTZq027aamhrumvIMI449e+dzLZ/7IhdffBGFhYW7laGsrIyf/exnzJgxA4C8vDwAzjrrLD799FPOPPNMKipCOVlRUbHzc+ic44orrgAocc6Ve9vKvId9kVBT/b4cB7wTdr0x7HIPoOkxf2pmNxNq2n9mL/u1eN9oKNxlv7pCS0WsB/81vebCwkIOPvjgVo9ijpVYBt2+xk60NCe8+fP4/X4u+f5pPDJ9BqRlUVn80W6tGcGqMnzZvfY7I2FfQRPvwZut+Xu25n1tXs7ln8/n6ht/y5FHTSKNYKvep5aep6WD6JNPPjmqz1ksTsxjdQEuOvsU8vPzWxXYDQ0Ne4R13759GTVqFMcee+xu23JycjCzVpUpEAi06gCxvr6eiy66iHvuuYd+/foRCARIT08nOTmZ999/n6KiIg488EC++OILGhoaePPNN5kwYQIAt956K0cffTT/93//V9n0eGHN4UcDK/dTzLOB58L/lGY2ESgGejjnKswszTm3A6hhV028pf322NaqP1QLFO6SEOI5+C/84KY9RiLHOuj21lTb2gVdiosXMf35V0nJ7EmwqoxfXfkjPiz+amdrxjWX/pAnX3itVTMS9lb+eL5/rfl7tvZ9DS/njppqvvh0IdljT2TgUcdAfd1+36d9PU/456ytB5AtHZA75ygvL28xnEtKSnjnXx8QtBSqK7dTWbaFqQ/8gYEDBtCvX7/dwnns2LF7hHh2dnarAzvS17GvA8SmUfKlpaWsWbOG6upq3nvvPcaOHcvy5cs5+OCDOfDAAxk5ciTHHnss1dXVTJw4kdzcXGpra3nhhRf46KOPmDRpEsBIM5sLjAECZrYWqAUuM7P+wGrgK8BvZrnA74BV3s/PzGwG0JdQv/p4Qn34W8xsC7DY+/ukAq+b2Yfe7UcASYCZ2b8IjdBv2nZbW/52CndJCO3R1xiP2mVLX+LNA6ShoYEdLiXqoNvbF2QwGNxjTnhFTXC3hXtaes0fLpzFrddfu9to88zMzFbNSNjb36C6uprGmu0xef+a/03Xr19PZZ1jUFaov7WlOeGtfV/DP2cNDQ0ESSIrI5P09HR8vh4tHpA0lSeSg6nWHkA2NjZSVlbWqtr15s2bSU9P32sNe2NVAwcecxrZuflk5fRizUezuOWaPZu/2yqSA5dAIEBeXt4en7cmvXv35sknn6Rfv35cddVV3HTTTVxzzTW89957O/d5//33OeWUUygqKqKgoIDjjz+eO++8k4MOOoj09HTuu+++phr+FkJT1UoIjYyfBvwdeBYoAOYDl3u3Hexd/xI4F7gZOAZ4DJgNnEdo1PsTwBagzDl3jpkdDLwHLAUeAiYB64HrgL8Cc4GTgceByWZ2G3CNc664pal4+6Jwl4TQHoP/Yl273NuXeFOALP98Pl98upAgSdSULKOk5KSov2hbaqoNb/KsLC/lvTdeoaaynHsnP861l11AUVHRPhcBCi9LLJqCy7dto2L2M+T0HRj1+9fSfPE33lvAZ58v5utgD46eeBQZvuTdDhz29hrXr18PQHV1NZmZmRQUFOw+4NKlUFOyjCOPPhafL7XFA5Lw8lSXbaSW1N1Wx1uLj6VLlzJo0CCCwSA+n4/JT84ga8RRONfI5rWr+OXNt3LC0ePZtm3bHoG9detWsrKyWgzs4cOH73a9T58+ZGRktPh3CwQCrL3jj/TqP5Ds3HzKNq+npnwLPp9vr/tH07IQyYFLS/uGf+aa315VVdXiKPn6+nouuOACpkyZwrnnnsusWbMYO3bszgOAJ554gsrKSggNqDsQeBA4llCIXwQEgZmE5p43OQ+4nVAtHXafJpdFqFn9NuAGYAK7svb3wCfAxcDDQE/gjrD7BgmN2L8TuIBQH/zDwJm0MBVvXxTukjBiPfiv+RdYS60DjTXbqa6uJhAIRDwAqeVlcfsSDAb5wakn8O+33U322BPJysjkyKOPZfrzr3LQQQfttc+6KYyaQqi55jXopqC6f8oTzPtsCb2KjuOUsyeR4ds1nz/8NaemZ7J1/WoaAuUt1qgjHZvR0t+g5MN/8ItLz9nra9jf4z30+DNkj5xAj7w+lG5ax+0PTOHMK2/kxNGTeO+NV3h9+sdMOHg41152wc7Hb+l9LVv/Df/1+/so/vJrysvLyO3dj8MPHMTN11292+espOQkpj//KsvXL9vjgKTp9eUc9G2wJNi4hg+em0pSVi9wjtXLl7B84bvMmvkk28pKSU5KYkdNgB07asnM7kl2Tj7Zub1wwRpW9c9n+PDhjB49eo/ATk1N3fl80X72ww9allTWsmzJIkaOGMkd9z60RwBH2zUVSQtJawaB7jb24bN5fP7FMwwePJjly5fz1VdfceaZZ/LSSy8xduxYhg4dSs+ePUlPT6empobk5GS+9a1vkZWVxXvvvcchhxwCoSb4DOAEoJBQoF8A3Ego2I8AfkAocD8CFgM7CE1tOxhIBi4hFNx9CYV4LqEDhWQz+xGQDawAbiG0qM03YVPshhBqKVgK4Jxba2YnAMeZ2RzvcZtG0u+Xwl0SSqwG/+3tCyx8LndV2SZcYwP3T3uhTdPLIFSLW1JZyw2330VmXj+qyzZywMjRjD3uGK/JN5Xl65e12EpQXLyIux6cwtK1W3DBGkYPKeSWX/x0Z1n29aVfXVPD9u3bIS2L1BaarQsLC7n6onP57/sfZuk3JZgvgzGDerNq1ap9vtZ9PWfTgcjGjRvZ4VJITc+kfOsmkpNTqGkI9du2ZuBi89s//PBD5n2+GN+6CjZ9/SW5/QawubyarRvWMvLwSZx5yUCWvDmD639y+W7nC2hp4ZsGB+tTB9KQE6TfkedSu3Ut37gUHnriWe69IxQwPp+PxsZGzj/9u3z99deUlZXxzDPPUFFRQVlZGWvWrOGLJcuorvkjDXV1pKSlk5aWytwtJWRm57C9ooKhBx1B+fZK+hz9A1KSjCOKRvHey0/yvatu2rmAzeaFs7j71hv2+bmOxViQpil+N9x+Fyf+6KdkZPWkIVi320FnU9dC7tjj8WVkEaypanXXVCQtX833LdtYwsynn+CjOa8zfvx4LrvsMqY+cDeDRv2Thvp6XGMjx51xPo899hjPPvssfr+fBx54gGXLlvHnP/+ZQCCAc47s7GyOPvpoXnnlFebNm0d9fT033ngjmZmZvP/++35gDqHa+zfAP4ENhPrJjV3z0/OA/oRq+g2EFpyZQWiOeg2hZnc/oYOBDGAQsJnQAjXvAIcQOpDYwq6wTgUuJDRX/m/AP8zsPe8xV3n3/ROhg4gcM5vhnCve199b4S7dQqT9fC3VGsLnctdu30pwRy1jz7giqv735rXFss3rWbZkEadfcT29+g+kdMNaXnnsXkZ/u4rsfcx1DgQCPDztWUrShzDklLNxwVpWz3+Z/3nsGe773a/3ORp+3rz5/OSm35I5+liqGxZQ7/OzoHgpE8Y27PZcw4YNIyc3l5OOOYv8foXUVlXs87WGB01jzXbOPeUEJk6ciN/v3+1ApGFHgIqt68kpXoIvO4+SLxeR4UvizmA1xxw2hk++/Gavp99taTrazNffpsfIo9nW4KP3kIls+2w2KX0PYP7c2QwZfSh1tdX0yPBRUFCwR5mbauMbN25k7dq1/GnKk2zbvJJt3yxh81efUl8bYE1tFcuSG3jzpRmUlpYSCATo3bs3vXv3xudLpSJQQ2ZOL/zpqZx+4gmcc845/PruByn4zmXkDxlJ5daNbHzrMV58/CFKS0uZMuN1Cg47jrffeJ2CCadR+s0y+g4dyagDR7Puw1cozendqu6JaMeCtPQ/EQwGqQnCB2/PoiE5jeSGHfjryrnxjj+SkdOb6rKNbKqooZ6vcEk+rDFIr8raVnVN7WtgZ0lJyW7laL6vLz2ds86/iMt+8D0uueL/sSEAvvRMio47g9em3oVrbOD//rqcWT2yCQQCVFVVMWfOHIYOHcr48eNZtmwZa9eupaSkhLfeeouDDjqItWvXUlpayn333cddd90Fodp504o2Q4C7CQ2iG06oVn4oMIpQ03kVofzMBg4D7gHmAVcQ6jOvA+q9H0fogOAd7z5fEmq6HwxsN7MB3nMZcDWhZvxcQv37k4FDnHPzzey7wE8J9e3fSaipfq8U7pLwIq3VNNUammqUGf5sdjSby71h7dfMeuohxqVnArv302ZmZu73IKKptvjQE89S05hCY00FI0eMpFf/gQD06j+QkSNG7vdLvry8nJrGFFKze5GWkQkZmaRk96KmAVasWMH9U55gwITTd54kpelLH2Dy9OdI7TucXoNHk5LmZ/UH/yA7LYV1ga/45dWX7jbALymjJ/0HDiUYrKOuvoGyypoWz5gXHjTV3gl4rv/9g0w45C2u+OFZPD7jpZ0HIsHqShY8dy/bffk0VtcT9OfjGmp565PlvDrrTcacdDHHHL57N0FTjb6p+T2/YDB1tdWh1pTsXhx66CG88cEn1FVuZ+v6b0j392DFljXM/ON/0qdPL46bcARPPfUUa9euZd26dZSWlrJ582bKysrYvHkztbW19O7dmx31DdT7sqmuqiRjyCGk9R5GUnIyDas/4ZZbbuaMM84gLy+PL75YzMPTQucjyBucwzEnnk52Ti/WLZzF9w89lMPGT6KUBratWY41Bhk1pojU1FRGjx5Nmv0jtGBNXTXbN5dgjUGCNVX0zuuxcwCZz+cjGAzus9tn/fr1VNQEGdDss7ivwN3bcsJ1dXV8vqgYO+R7pOf2p7Z0HcvmPsWhp15Iv0HD2bhmJa/fcyuHDD+WvIIhbN9cwrL5f99r33xLn/nwcTHHjx/LHfc+tMf/ZtP0y7889gyNqZlk++CK889i+vOv4u8zkAOOPpX65x7j5UfupLG+nszsHlhDEOcc9fX1+P1+8vLymDFjBjfffDOlpaWkpKTg8/nIyspi8eLFBINBzIyMjAw+/vjjpmKuI7RwTBqhsF1EaKBbf+AaYA1wAPApoZp5EFhLaD34BYRq4KVAtXefXEK1+xTvsQu9x+1FKPTXEGqiP9S7vsF7/iRC89x3Ts8DfIRq+5XePvukcJeEFk2tJicnh7L13/Dp9EdIye5FfWUpvZOr6Ttk5M5mwvx+hZgvg63rV1M4fPTOftp7Jz++s7bZfP315pwD1+hwDfX4UpJpaNixW60m/Et+XwvMZCTVU1dZyo6aalywlvrKUqrrG/jDX/5K8erNrNj2dyYefyIFQ0fudjKXipogG1YuY1utI6mxnp59BzAyo5q7b71ht1XXfD4f1WUb+XLRp3z2xRJKviqmrqKMiq0b+M31P+Ooo47auW/TgVF6Vg/emT2HjPyB1BccSPbICTwy7TmCKRk7D0RqAxWYL5PsrDTK1pYw+LjzKV2/hkx/FoE500nvN4wFxUv5t2Mn0Rh24LRw4ULeX/gZ6Wu2EawoZcjwA9i6cgWB8nkkL/iIzatWUrWtlIbq7VQ1hCpOn65LJrdXL5YWf8aQgQWUVdZgaaFBfEOHHcABReP46eUXMmnSJMyMRYsW8Zu77uX9xauxwtGkZeeRnZVJSq6f1+fO47TTTqOmpoapT8+kx0HHkp86An9eH+bPnc2p519Goy8UtPnZ6YwYe+DOJuxtxbtqqE1B1z/T+PKtxxg1pohtxSU7l8edN28+k6c/B2nZWLCaay45j/Hjx+/2/hcXL+Lhac9S/OVKlq9/hGNOPJ209EwqN69j69ate+0W2dtywrUVpdTWN5DlGqiv3ELDjhosM5cUX6hf35/bm569+lDxxRxq1uRjddWMHDGS7du3t2qt/vDxCj6fjzvufajF/82VK1dx3+THWbxyNY0kc9DQAioqKthStp0dNdXk9ing8NMugfwhVC+Zy8jDJzH/H9PxmaO2tpaqqio2btzIGWecweWXX87bb7+987Pcs2dPSktLGTRoEOvXr+ewww7jkEMO4W9/+5sjFM7bCPWZbyI0h70AmAgsJFRjH+7tl8auWvnV7Gq+/5DQKPggoWZ6R6gZfi6hkfUF3r6bvMc4m9DBQAahfvyehLoG/gcYCXzt/flSgXu950vb6x/Zo3CXhBbtCPek5GR8/YaTmt0L8+eRFlhJUl1gZ/jWVlUwZlBvKr+cx/J1y2gIlJOUnEzBUafv9oV56KFHQF0V11z6w50h2NTn/NDjz9D7sO+S7EulIVjHindmsn7+KzSkZFIfKOeaS87D7/fvtkxmc36/n2svu4C7HpzCktc/xgVrGFHYh+raHVRmjaC+X29KM3N455//4KSzMnc2g27YsIEVK1cy+NgfUF5VTbAmwMaPX+LByffsFuxN4bG1qo75k/+ApWXhH3ooSWlZFJc6rvzV7Tx4x42MGDGCnJycnc2pSxa8x8pP3ye1Zz61JcsIjj8KX1YeDRWl1DWW8s3nH7Ly47epraokMG82Gb0HUldbS13ZelIbcqkPbGf1v14kWFPJ5g9fZOs3S3n+0T9TXV1DoKoCS0omvWc+aVk9WbP4Y/J7ZjFp/BGsWreRpKRkUnIL6XHU90nukU/N8g/Jys7hR5dcxvYt63n5ySkMO+48yiqrGHzscBpWf8LQb5/Bsy+/waGHHrqzufrfr7wU/vdJ1qTm0nPQSII1AdYsK2dxdRI33H4XF51zGo2+zNCB3perMV869ZbC10s+IVi+GWDn+vrhNdOWFq5pqqE3heP8jz7iJzf9FvKGUF62jtzCYVx1w+1Mufu2nYsqNfV/Fxx1Oqcc3sD7H87npf+9F5yRlJnNJb/4zR7jL5r+J3bgI+iMdevW8fny1aT1P5CCw46joqKCmnc/YGC2n7LtlRiO2vLNrF27huycXpSWfE2aBfn2yd/b2Se/8t0Xuefhx0j257SqZazp9S9dupQd+Pb431yxYgV/+p+/sqpsB8NOvxbzpbNxxac8MeNF3n39H1x060OhmSRzXyY5LZPCgYPY8M1XrFm1gtTUVI488kjuuecejj/+eNLT05k2bRpJSUk7B88ddthhJCcnk5OTQ2lpKeXl5UybNg1CteWHCPVv1wLvepdLCAXxYYTC/UvgcHblp59Qc3m1d/10dp1RbjuhmnsFoVp70wFACaGm+97sqqEnEWoluBNocM5N9AbSpZjZeGCbc+4IAO+sdfukcJcuo6UR4ZEsndra+dPl5eXk9B3IoZNOpra2lvT0dL754FXOOe4I/j571s4mxZuvu5phw0JfztXV1dw/7YXd1l93PQtZunoDSZk9uPrG3zL1j7eTnp7B1KdnUlnn+PD9D8lesYHM/AKsrpr+memceuxEXvjn2/iy8rhv8uM0NjaS3qMXGcmOn11x4c7na7662V8fuHu3qVuX3vB7hhw9juzaGtas+op1a9aw+r2/c9ZJx3HHvQ9RUROkqs5RXbKWpNRMXF0Nhxx6KAcccMDOv3N1dTW/e2Aym/zDsT6DyRgeZMfGVSTX15B7+EnsqK7CKtdz7X/9N5OO+TZpVs/F55zKSceM4+a7/4fk/oeTnJFBQd9CZs+cRn9/MmNHDWP+G0/y9Zp1JPtzoDFIfaCCqm8WsfXT2SSnZ7EtPRMfDWxf+SmWnEzBoP70HDmGkceehT+vDx8t+pLgppVkDy5iRwOs+teLjC4aRXa/3lx3+un86g8PY33Gkj7oICwpmdqvP6O2fDPOjE8+fJe0gpHkDBlD5TerqGzwkZmSTrIvlUZfJk8//TT3PfY3yoPJNAS2kZ20g81V8/HlDWRHxRZ6DCsiUB8gY8QEZr7+Nq7RUVtVwbixo5n1ykuUfP4upauXkRSs4dd/fJjeeT32ujzt3j67gUCAydNnkHXwd6n29ycPCCyZS8++w/n9A4/Qr19/kv05u02tywZOPOHbPP3ZW/QYNYn+B0/aY/xF03OsXbuOf82dQ82KCoL1jTTUB/GVr6Gquob+AweT06sPy/85jewDDscFaykcfQQfvvAon87OITnNT//8Xix782/0HXwgyQ01uMYGCieesXM2RfPna664eBEPPf4MFbVBFn/xBRmDixg6Ylfr1z2PPM6nqzZRUVlFbsU2cgYMJ8WfwzuvvsEtN93AspULqCrbzpBRB3P8Bdfytz/8kpLlS0hJ8XHIIYcwd+5cfvnLX9LY2MiaNWtIS0vDOUcgEODSSy+lZ8+eTJgwgdzcXFatWsWqVaua1qdvBP6TUM28gVCAf04ooJ8BziJ0etfRhEbC9yHUT+8Ihb+fUK38HULN7NmEatipQD6h1e7WEgrz/oRq9j7veZ2373OE+ukDXl/8ZOARQn3yi7xtjbRivXmFu3QJLY0I/8Fp393nOa8huvnvTQcEtVUVux0QTJw4kYkTJ+7xhdx0kNEQKKdk5VLSe+RRE6iirORrxp59DVn5BaxfPJ+/PPYMmRkZFE48g7z6ILP/NR/8hQw54iQC2zazdPaj/D0tiWHH/YDU9Eye+ct/s606yOADcmkIbONXt/6OIcOGtzjIzO/37+z/Xr58OQ07AtRWbScrrw+DhwyjcYlx7RU/4vG/vUT2yAlkJqdQM+ddemT0ZMiBo6irKqfknY/5/PNinnzhVZau3UJN+VY2by3l8ItPILdPf0p3wDfL5oMvne3vzqS+rpZg+Qay0lJYvHQppRvWMP2vD5OcnERVRQUN9W+R7EsjNasn1lhPn7FjaGho4NhJ46moT6b/hNPp2buAjJ75LH1lMsP65LJ2Ww3bt5XSIyOFkYMH8JPLLmTUqFHcP+2F/9/eeYfHVZz7/zPb+0pa9WZblmxZ7p1mMKaEnhBIgNxLQpKbkITLTb3c5CY3IaT9UkgFEgghofcSug02timuuMuSLVmS1etqezu7O78/diWEG8Yx73dm+gAAIABJREFUGMR8nmefPWfOnNnZec6c75R35qVm8VLCQR+6HQ0k0TG7diIbtuykYtJkTrn8yyRiER566k6CAR/W2mJSwSFS8Qjx7j0UVk3n1dfX0TUYQEvEiAaG0SGJ+gcxRwOktAShoW5++/IKUnUX4PJUkPD10/PiX3DWLCClM2Ermkxw3ybyF57FzsYWqm0GLlu2kMdfeBRMDpJd9Zx39fV0ekMYXQX0NL5GzZLTufeJ5w+ydj/c/vRGo5H29nYwOzBIE8mkhj2vBF9akor4aIlEqDnjUsomTxs1uqw5tZu8wlJ8/T3ozVbsnpKD7C9GRn+6u7v56wOPY50wC293OzqLnXDzm+TVncS2hn2cYjFRU5ZPk0yRXzYBs9nC3BlTWXHPH1nyyS9g1sP6V5bT6/Xizh3i42efxtNrtzDU28XWDa8jTTYC++tZt24dZ5999kF1KxwO8/Pf38b+QBqD0wMWJ8vv/gOnL12GLhFGp9czcckn6HY0kOrtZ9+mldQ6cuje+gpDfV08/PDDAPz3f/83N7z6Io/94nqGentx5+RisZgZHh6moqKCOXPm8JnPfIaVK1dmPOWlM1u2t7e3M3v2bFatWkU8HiedTlNXV8fNN9/MOeecM0RGEzWgicye8puz518iMyzfTGYv+XlketpfJdNTX0rGSG5km77HyfT0R8J6yAh8KRlB12XT7iNjST8dWExm45tTyTQuHs6mtUxKuU0IMWtM2HWHfYFlUeKuOC68l/utH8oivOX1x/nx7//Kx7/y3cP6vB7hSOvfD5Xvd2oQjIj5WAvfXfX19PX3s27XfUgtStg3RG5ZDY78UuLRCBarjTQ2oikI+oZ4dfkzCLMN774d7NMncbrzqJgwkaTOjMlio7N5N0NeH7mnfIrcSZOIeHvY8vgtTDvvs1RW1Rzx/w4NeUmE/ex64V70NhduY5pZUyYyMDjIG1t3Ipp66NvfRDyhEVz7MMN7KkmGfJiSQa795v+gyy3DYLYx0NNFxDfEq7+/DpPRRCIaRAsHCLfXY3B6MDk96EwWQt5+EtEKkuipOO1SaiZPZv2al7BOmMmESZOp9NiR+7dw0/e+RTKZZPXqNTy5ejNhYy5p4UALBHAY9ZSVljDvwrMI+oZYt/pl2n0xnlu9DqfTSdQ3QMPm19iyaT3dg34C++sZaq3HZDJx8Wevw2y10b5nJ2/uaiQtDPjffBZb8SRIJrA73QTbdpCMDJPQ0phMZppefgBPYQne5h1MO/kk9q19ksG+HoY1PSa9FZvOgNFdhM5diIYRW+k0kmEvmG3s37mOHkcuLb42AoM9eEonEBjoprq6mslTp9G1aSfuwjLCLTYSWpI+f5jm5uaRtdT09/fzuzvupnDuWVhsDvY37uAL3/wBNdOm07ZvL5MnVrGvrYVpZ15Gx+49DO3bSWjvm5hKSkhGAkgE/Z1tuD2Fo0aXvfYcor4BppR6GAz73mZ/Yc3Jo6Ojk1/88Q4CUY03dzXhmL6U4rK5mMxm+lMRkj1NdEcH2NT6GhaTAZNBj0xpzJu9AGIB9BY7peXlrHj8ATwLLsQ02EM6x8wPb/4LmKz4X11H7QVfwFlQirDn8dgLq0ZXSIylu7ubhrYuJp5/LRZXHnkBLy3P3spnLzqD4uJifnf3E+QVlrJorp616zezr6GD1udvZ96UCTx684a3Nd4vueQSbrrpJm666SYCPi9WqxW73c7AwACrVq0a3XiotbWVe++9l6997WvEYjGklDQ0NLBgwQLC4TCNjY0jjmN0ZHaGW0Jm+P1SMiL+Apn58DlAmMza99vIGM79msyyuXPI9LKLyRja5ZHxJvczMhb3UTJz6XPIWOI/RmZJ3C4yHud+BQyS2RinHrhSStk+tuyyS99OPcxr8iCUuCv+Zf6VNbZH0yg4lEU4Jjspk/MdfV6PMNbye+T8SPk+UoPgwPtOnjWFX//1fpyzzqWgxMzEohyGd6yiqW0/HdvWohcwtbKY0MAwsXicNS/8k4KFFxFqayEajdC74yXE5DqG+/cR6NHT2NFP2mAmGPBj7G2lTYuSCAdI6K3EkvKg/xsOh0fdj0LGCn7GaRew7vXVxPoD+AbbmF90Ov/7o5+w3xcn5n0TmUqS1uLIZIxwyxb0VgdCprHmFpEa3o3J4cJVvQB9KEhaClJD+3HklWAWKSyeMpLCSFH1TAZDMby7XqVjOEQsGKW43MiuniCOGWcyuOl5Eh072eUbZPbc+Zx35ReJhAN4/RFSWpzwir+ht+dCeJh50yaSSOvRtATPP/YAjlnnkNZi9Ekr37jxV1RWVLL+5dswT12CrsCDp3Ayw9uWo08MsHlXE83121i7/BmMhZMR+jDO4slEuxpxVUxBH+5C56pg8sVfAZ2J9r27GNj0DJUVufzfTTcwZcoU/vi3e5ly1qfZdfefkVqceMhPKh4mFRhAVzyFdDxApOVNZCKGtDhIGSyEzR62hV18ctopuKpjPPvXX+OZuptoyMdQVxu+vi4effo5Ii1b+Pf6Hdz4jWupqZnCb/9yF2/s2Eti137sVhudTbtwzziD+oid0lOuZGD/FhaecykbXnwUnTQgE2mmnHI+Um+k+eUHePCOP2Cwuch1WJhWmsPnPnUl/3j4KazufJJaN46+7bQ1bRod4friVZdyz+PPUTj/YxQi2Nm8n/7ONpwTXaTTaUwmC063i6S/j04KcOYVMv/089izeycvP7id+VMrqKsswNffQyQl8PV7CXV10NwUBXMh+UXFQDed7S2Ux0OcevJiBrevpqGhgWnTph1UH4XRijBasscW9BYnxcXFlJaWjk6hFRYWseykuVRp+/nOVz9PdXX1Qel4vV5efvllOjs7KS4uHg1/4YUXRr3A3Xrrrej1ei655BKefvppnnjiCaSUXHHFFdx///0sW7bsrXwJEQT6pJTnCSFygD9IKbcKIV4ks4/8YDbsCSHEJcBPpZTN2Xt/D/xeSnnPmPS+ATwmpfxJ9nx1Nn6zEGIt8O9ZZzGXAiEpZfkhX4DHiBJ3xdt4tz3wf8Wd6NE2Cg5lEU4ijD4RPCqf14f6rRFDpyNZ0R/KkceB+4N7+7v5+S03YZ8wk9Lpi4lHI7T37GNi6QTOXDybvzz4NHGpp+HlPlx5haBFCSbSGCtmY0xF8fqHSdoL8Wp6HMV17NqyCkftJIymXCxVcxlqWI/ptMsxWJ2Y0Nj8xqvo4wG8vR307HiNr39jM6vf2EQyLUmEfAgpicVjpFNJTDYXFrcHUhorVr5CXJhwFlcTD4dxzT0fYXZgNJvxb36WgspqXA4z1qIqOhq2YnEX4ll4EV3160lGw0SD/UQDg1jmnIvO6aaiYhKRPa+TiBsR9jxM05eSHuygq2k7rrrTSHsHSMQiaMKAzppP/ZCGpXwBgd2vYSytwDPrY2ihIdIDLQh/FxNOuoCNzz+CeX8XlrJanBNmEuxtZdPGl6mctYwpc+toGkowqKVwFE8grTNgrZxJqm8P3j0b2Nq8B9f8SyidfRox3wA9rz1CQY6dmbkJrrrmm/z0r49idXnQUml07kKSiTi797Xy8MvrST71Ig2tXZRZa3EVT6T79QfBZCMVGsZi0hNpXEvcnou5agEGlwd0BuLNGzDklmHIK2PFqjXEUjAYg3t/fxMmq5N0UkNvMJBXNZM5n/gKwmThx7//KwvnzqBw7llEN+7Efeql+DqasE40gMWNzu4kkDRgN1gorapl+owZJJNJZp//GeJailff2IDZU4p7wcUYnXmE2raT0AI8+NSLVC391Nt2+fvFNVeO9lxHjEqDviHeWPkiiZTEv2MlsrueaDiE02IkKdKcduV17PfGsOcV0rLtJc7/9DU0v/YM373+88RiMX7/17vp2LUBt96F3mRGpzOi+XspqjsJX88j5NgtnLJwDkMdzWzauA6EwMwzb6vXpaWl1FUW0Nm8FZPTQyI4RF1lwUFb+46MmN1w3RcP+U44Hl7gxgr7IZgKTMkK+0ygnMyStQPDzhZCfBGQBwj7uWQs5q8Y+xo6hBe4WWSG2C88UmaOBSXuilGOpQd+rO5ED2wUDPV0cPOf7+JXN37vbdba8JZF+E9+exu7/vk6UosxY3I5V179paPyeX2oBsif73kEg8190Fr2Q/U2RgyAoilBOh5GS+tG9wc3Wh3oHAXIWBBfXweOvOLMfHuwna3hYS76wjd55YV/Ypz+MYRMUpyfR/1zdxH19WOxu7DZ7dAfI7h/N95IgLSWxLtzDXqjEYNOEOnax77mTchkHCElw01baV7hwuV0Mr12Cq9u2ELB7DNxFpbRvmM9upwS6ibVEDS4kd0NeMom0rb+RWxON/H+LrDlYSmrwTH1ZLSgF6PFisnhwlM1ncGmbWjWAIG2esLGJjSdESkM5BZXEKlfiT6nmITUkfT7CIZ2kupsJ9zfgWXiLLSQF4SOZGCI4fX/BL0B27QlGNyF6AwWgjuWkzTnkjbZMRhM6GwujHo9WthLqL2e1avXEBoOoAtFMebEiTXtRApBMCHp6usnEJ5EKOCDnBxSehNaYJDIUA82s4OJkybT0d2L3mjAN9BDTn4xdnce1Tlpfvid6ykpKaFu+Wo69m6meziMv2Mv8YRGtOIk3vSayLOV4gs0kmxrIZlfizmqEd+zjpy8fEomTaWsrIJ1GzYgjCaIhTAVVZFy5RMf6iI01IMvoZFbWYu0e8idcS7BjU/gmrYAX/2rVC06m1RSw+7KJWVyEogmKXe4KZo0lWjYRxpBMhZCp8WJhCAa15C+flJaApfFSDSWpKu1CVd+CZGgD7PdTe2M2UgEAUMCvb+RaEq8zflPNCWw2Wxv24MgFfax/NF7iOgdCGsBtsIJ1BZbuPF/bkLTNO55dg0102fSvuYNUmlJTEsx1NM+uumP3W7nu9dfS/+Qn95wB72DPqKBYdwzTseWW0hexWRCO1fSYQ2zbetmTv/0taNGcmMbzXa7ne/917WZ/R1iAazmJNdde+0hVxAcqYPx6KOPsmnTJm644QYAfvGLX3DdddfhcDiYNGkSP/7xj9Hr9XzpS19iyZIl5Obm8sADD9Dd3c0vf/lLTjnlFJ588kmuuOIKvvrVr/Kzn/0MMkPq3xJClEopbyKz/A0hxD94q5f+trBsdm4DNmYt2NdIKX8E/ImM0dsrQog9UspryQy930NmiH/E49uvySy7Wy6E8Espj7gxzbtBibsCOPZdrt6NO9Fb//Eg3/6K5W29CWduPt2te9iwZiV+7wDf/uHP+dzlFx00X5dZE54mrcUhGUemJVOmTHlHn9fhcJiGhgbCSQEGE+FwmEQyhT+awLdvM9v3tmLOLWZ4fwOR4X40TcNhequ3EQ6HuenXf6DVp2HJLSLuD+Btq2fCorPIKyqnb38zgZ59xFweup++M9OD9nZBXS3bd9ZT5YvTv7+JUEcPscEu9kf9xLzdDO7ZTDqlZf6Y0CMMJgyuAvRWByZPOWlfL8b8SmzJJHkzz6SocjI9+5uJNW/Ak+NgRk0Vl194Nh13PUXVOZ9hqLMF6WpDlz8JTeqIDXYS3r+Xru1rKTjpEwiDAVdeJT1bVqEz2wjsXInOaCaWiJIODdG04n6E2Urvnq0Ik5WUMDK8fSVSbyZstqBFI+iEkVTETyLkxVgwEYTANHEuiYFW9DY3WmCQVDyMvXoh6aSGMa8czduJuXQqOouDVMRHKuLHYHUT7d9PKtBHtGUr1kkLIK+I3IJKgrtfQxMmYk0Zw7147z4iuWWseG0T5px8BjY/Q6RpAzqrk3QiQjw6RGuLA4PFjrOgjGh4mO6uPURbt2NYPJ/bH30RnRbhwqUn88izL9G4ozljA1FWi3nCbBIBL0ORfoTZgW/fdizFVSSHujC7PFhdecR9fRgK3OQYNMLxANKWQ6R1K9H2ekxWB6n27UQiUZL9LZjKpuGYMIv4vg24y2sZ2r6GzY/dhtlThhYYwJEO4po+gZSWwGYyYrTbGO7XMOSUENy5CmdRJaGuvRgL8hncvpIpFYXc8fAz+DfuIh0LYxUa7pIqIr4BTDYXybAPl0mHQS+P6PzHbrdz0bJTeOqlNRSceQ0mm5OJc0+lZ819eDweSktLefjZl4mFAlTmmFjx9J1osRCxniZ+9I0vjdar0tJSJpSXMq16IWuefwrqFtK3dws96/zInr3c9vPvY7fbud1kYlLNNODQ02UzZ87k5psOL+BHs4X0VVddxVVXXfW2sC1bthwU7+qrr+bqq68ePXe73SQSiYPiff/73+cHP/hBvZRywYHXpJTXHClMSnnQmnMp5dRDhHWS2bd+bNjHDsrMcUKJuwI49vXghzM+O9CdaCSRYv2uffzkT3/HZTXy2csuRKdFGOrpYMOaldhqTyPV30FbKjG6o9l/fuEzowJ7290P0uesofbyTyO1GF3ZJT4/+d43D5u3kZGInr4hXtuwgeIhgW+wl2hvG/GQH10qjtkeJCcOXW37cZRV0dw7zKz5i0YbNqtWrWLt5u3kzL2AoUAALZQkMOTjbzdeTzKpEY/F0RIxhN6ATMRIJ+MIoWPzYCd6s53+zv1o6NDlT0Ln8JDWG8itXozNZsWkk7RtXoWpoBK9yYa1ZjGpwAAGm5vw7rVow11ULTqHnj1b2bN3IymDhaIFF1FRU0vXzpXc/ehT9LTtJfjGK0TCYeLhEOlIiHBuFdH+nQQG2rEUVRHt34+rqo7+oBG9LYdUoI/InjcwecqRYS85xZV4O1vRgkHMJTW4T/40aS1GtHkTsb2vY5u8kJTchWPehZBOk3YXEd69GkvxZIyOfEglwGBFZ3Git+dgyCkl1tOITGnoTFYSPU1o3i7isb0InY5Yxy4iHfXo7TlILYbe6gCbG0NSD0Yz8Z5mdPkTIR7BXncmsc6dxGWa5HAPBqsD66S5mHIKSAz3Eq5/hcFhP8WTZxLYvhzMDgLN27AbJF2mCiZNmEs6EeXXf72dyTVTSUR8JNMCnZZAC/kxuPKJDrSSTkQoWfpZdEKS6GsjZ/YyivLzWDhtEr2bXqDAacbfsR1NbyHh6yetxZk4cx7acDeRwR50BRXEepsRRjOpkA9kCp1ej3P+xZg9pSSGexA7n+Xy88/iudVryTcleX3F33GVT2Gwsw17XhHBjgYqq+uQUS+nza3jN3c+wKSLvoY9t5Dh7jb2PftnIgMdbHv0j+gMJuZWl/GNH2Z6rrfc9QDrt9fjrjudpZcuPmhXv9raWopKyiivKMfmcJFMxBFG69vq8C13PcDm7fVMnH0qc2bPJsdh4ZWNr3DeeeeNCu5IXZ+Y72TP3ldZUDMVhxG++o0fsXjxYsLhMGaeecelp8fLB4Ti8ChxVwD/mj/0d3InanG4eGP9BqzOHOrO/jSJWIR7Hn+Oy88/k1vvup/h/gDGkgGEEHimLICoH+fUxaMvp8NtsdrdP8DX//cm4hhJhr2cMnsaS5cuHR36v/XvD6AvqaVlZxMls8+kfctapExhKJ1GyUlzCA32Ea9fwf6tazDmFBHqaaM14qN54yoM6QQv//Nh2js6iUajDLXUI8w29LYcdOkElhwPyWAAndOBvbQOU34l6agfQkPkTjsZo9lMcPerBAe60KMjFfajy7WiM9kwVcykwGkg5BvElFOEa875xNt3YnDmkw77SXi7SAUHiKTTdDZuJR4YQmfLQW+2YXXlMxRJIaWRjTuacE6eRf/m59E7C0mnkqTbd9GzfzvacDeuORegM5lxFlbQuXUFupJaZCoBVje2SfMxF1ai1+mJtW3FWDARg8VFKuQlFQuTCg5iyCkGq4todwN6kwW9zY3ObCfR34bRUw5GKzqjGWG0YCqpId6xG5nSMJXUYMgtJtzwKqmwD709B9vkRUTbd+CYPBet5U2EZwLSbEdnMBPprMeUW0zI5yUaDpOMhrDac7BPWwJCh9bfQmqojeLpJxEIR3FWzSHRtw9XfhGmybMpdlnwzDiJQH4xw70dBNJJwvocuoIpHnn6BQrzPVhKpjD5pPNY1x4hVv8qUm/E98bD6G0udIEe3BY90b5WkmEvqZRGPCnp7+/HX1nC9n1dDAxGcZZWUzp1NgajhfoX/s6cOXNo2i2wTTuD7u4udJqGb92jOHKL8W5+Dnt+CfNPOoVAIMiw1Uzv3hzuevAJ/vs/v4ROp+OWnDzK5p7B2vWb6fdHsHuKKV+4jOhgB08uX03K5MJdmHFxmls6kbjRzrmfupqiyiqCQ73EWrZQVZUR7+987Qv89LZ7mX7OeRizu8mNbZyXlpYyY2IRne27SRww1z1Sh99K48LRNIYO0es+3MY78P64XlYcHUrcFcCxVcoDjeXGxh2bXjAh8e/exrlXfgmz1YbZaqPeG+BvDz1F2uIk5mui0JRgwOhGajFEIoIrr5D+3Rrd3d2UlpYeZFAXH+5lZ2sDcsIiwsFhIj0trHhjO+L3d+JxGJg/vZbXtuwkGo2iaSl0BgOpeAQhJenmLQy9fGdmDlVKDE4P6VQSndVFLOTD4siD6DCFk2pI5ZTj0/QYa5YgdfqM2DSupnLBmbRsXk0KHZbKmejMNtKOPBJJDd9AD9aSGqTZjlEH+opZiEQQY/FU/I1vkE5pdHf1QTKGjIUxuQvQTV5IuOFVtKFOdDpBOi3JWXgJkdatOBZdTioagHSS7vXP4K6Zh3/vThJGG0ZjHul0mvRwN7Ypp6I3Wwk1b8BcVI2xZAqpiI/hrhaiAS/pvuU4Z58N8QiGkqnIqB9jYQkRgwVtuBejI59UYIBE3z6sVfNJ+nrRuwpJDnVgdLjReprAaCEd8ZHo3ovObEPoDViqF4FMI1MJdAYzofqV6O15pBNRpMGAdfJCdGY74Ya10NtELOhDr7dgcBWQlH5SQS/ezc+ii/mxFU0gbs8h2vom2lAnpJOQ1khLwFWIOdyMQQshbDbS4QClHjdTpkxl7XN/IxJPEBkeQJjt6IQOQ24JKS1Bd3crlnA3Xn8Qc04B5pIpmEtq0FXOJLTlWXKsem6/+Wfc8vcH2NDSCJYc7O48ysrLWfP6a9gtNuzVC3FOnMVA8wZqTruIPE8hwaaNhPxRCifYOeXCcwmFQrwebqO6tACnzcS2XXuI+QbwR9KZaQqpMWnJx7nn8ef44bevw2ESmC1W5syezfMr15IM+9H8/Zx68mJ6NgdI71uHv78Ld2EZw91tyHiECVOmYXfm4Ckooqmv5W3i7TQJYqEAxkM0zt9prht4xzTG1u2j3WL2vVgaqzg6lLi/B7yXa77fS95NpRwZ8s70mn2H3Pd6JL3u7m5+c5uGM8cDwFBPB9t2bMcx7XT8SR3x4tksf/hvFJRW4C4opbxiAk/efSupRJxf3XonV158LssWz6HxrnvZvOohtIgfi07S5w3CwBCJoS5IpzK9Up2B7iFB1/42hMGUmcMunYrRU046FibRuxfHhLmUnnQBsb5WBre/gm3SXIyeisya7R0vEeysp7yqhri7kuF9r5HKr0J2bCfmH0KmNIy5ZXS0NBHzDaAzO5DpNJq3G5mMk/T1YvKUkxjuJtS0hdLSMkJ6A85Js+nb+Bw6o5XAxidxFE8g2N6IMBgJN29C78xHZ7KATOJefCXBnS+TTiYw2HMxF00kNdxNIhYlGhgisv450rEA1oJKop27sVUvRphsxLsbSUX8yGQCKfSk41HMJVOIRoNImUbY3JgrZxHc8hwiLwxaHP/eTcS692ItryMx3EMqESG0YwWJvhZ0JiuO6UsJbn6aVNhHtHEN0mQjHfZjEho6nZnceeczuGcT6ZCXUOPrGGw5mDyVSJnG4Cok0b+PRG8z6XgYkYpTWjuHPT37cM5chqloMvH+NhJdjaSCQ1TOPomEvZh481Ys5XWkA4NgMKJ5B9GRoqSyGsonsn/jchxmHYZklDnnf4o9u+upPutKGl96CMP0M9F5KhA6GH79EQw5xSQ6d2OatpjtDU3IpIbJoMdqd6IzmEjYHEyvm8TMmTP53U//jxv/9HfsFXXs3rKB4O4WfLtew1o5iXRgiJ7NyzGkEww0bGD+9Gp+9r/f5qabb6V81pTRvRaM86fz7a98ntLSUl588UV++NvbGUroMUqNsy74BCUTp9C0P+O0ZKTxG8eIvXcr887+OHULTiMWCuAwCf73umv49V/voseaiwwNMW/qBNLJJMAhxXts43zEG9+B9fGd5rqPV69bDbufeISUR+37fdywYMECuXnz5vck7ePhV/mDxKEaKuFwmG//8GdoxTNp6/eRiMcJ7ljB7b98a97twCH6devW8ehzK9HbcxjubuP1bQ0kKhaQTmokg14ie17HTpR4PEE0FkcYjJBMILU4JrOZqkkTcTgcDAwHMNqcyJRGe1cv+uIppMLDuE+5Ep3Zjs5kZXjN3aRTKawV00GmsU9fSrx7D0hJeNdKZMSH0enBKBMktQTGCXMx5pUhkwmC25aTTiXJLasiEfASHOoFqwu9IxcZ9mXmzSN+ksEB9E4PqUgQncmC3paDTETQ2XMx55eTGOokHRqktLySoc42MFnRUmAvrQGTjbTeRLy/FSJekpqWEV+DGeeMpehdhfhefwjbhNnEvZ04Z51DMjBEWosS278D5/QziPfvJ9nTgBaPk7PkanQmM6mkhm/tPaTjYexTT83MUbsKiXfuwjpxLpGmDThmn4vm7yPp74OUBlqMVDSAzuzAMf0MSGrEA4Novj5yT70SocUYXH4rrvwi7HnF6MtmEBjswWIyMLDhaexFlSSlHoPFTszfj7GoJjO/bLSQ6G1Cb89DZ3UitRimVIzSkiK6e3ow1ZyGMbcEmU4RaXyNWFcD5ed8gbSrmMj+enxbX6Bw2iISvl5c+SUQ6CEcS6BzeIj1tfDtz3+apUuX8ue7H2Jraz+uSbOIDHQRtBXjC4axVc7Av/UFEIJUxw48E6cR9fZiExqJvCpck+eRCnkx9TdwyozJ/PanPwDgf276FYXzP4bJYmN/43Ye/9sfmPuZ/8EMavkqAAAZZElEQVTozGNfYz2Da+/n/DNP4b+/9h/MnDmTnTt3cvt9h6/vra2t3PCjnzNpyccpmTjlIP/sI3Wls7OTe594/qB0+vv7R/cw6OvrO+JvjdTNN95Yx2PPrzzqvd6Pps5/FBBCvHkog7oPK0rcjyPhcHj05TAybz22In9YGKncHR2d3PfkwS+cl156mW//6nbixTMw51dSXFTE8JvPM8Ec4fwzT+WJ514iHNeI+Aaxm3S0dfYQjUaJhUOkkglisRgpKTC4CzHmFCNMVhK9zRj0Oqx2O3FhxVpWixYawqyXFBYUMntiIel0Cs+cc2hs72Ogt5uGFQ+gs7lB6Mk59UqEwYjO4mD4lbtAb8JaOZN43z7stUtIhYdJh4eQ3k5MeoE1McySiz7NKy+9wEBHG8LmJhXyojNZ0dvcmDxlRFq34ZhxFrGevZBO4Vp0KTq9MWMQ1rQBW+2p6MwOglueQaaSGJz55Cz6BGaTCf++bST6mpl1+nlYSLLxiTtJmawYC6tIJ2KYi6shHkLnKgSdkVjzGyTDASzFVWhhPzLiR1icpBNh0skUhrwy0mEvlsqZmHNL0KJBwrvXIFNJcpd9EYO7CG2wA/+6R0h4O7AUT0GSRphtGHNKMRVNYviVf6A32xEWBwarA1NOAXq7B623CYvDianmFPD3Eg36CbVuRRhMiEQEzdfLrIs+j9/nI+EuB5ONuqlTaFv1AInuvaR0OjRPDaSSOGpPxde6k5SWINiyBUv5dEy5xVitNsJbnuHsUxfiD4ZoN5ThDSdIp9MkGldjFCmc1QuJpA1ILU6is57i8orMXp7xEHpnPuVLLiepJQn17GOmI8TNN33/rTo39yzWrVqOLJ1B+749RMNhgo2vYTSbmbH0k+RNnknr87fz4+s/x8PPrWT3/l50KY3pkyve5lhlrFiP+i4vm4fUGZHJGI6hBn7zvevftszsncTwnRoAR5vO0cQZL++gE8F4E3c1LH+MHKqSHa3F+XvdMv5X0h873L7+jbXMPu1jeArK6W7Zw5f/8xucNG8WK1atpbt/kNj2jaS0BLvCw8ikxi6rjdUrnsNTUY3ZkYM3Ihje24R94mwKZ9dgNhnJCeyjoKiEFeu3E0sksRROROh0kAijhbwkjE6stWdgLq+D3mZinfWYK2bimVHDy08+gMvYjLN0MrF+L/ZJszEWTyG8Zx3x7saM0KdTCL0JmU6STkTQmWyEtr9IOhGDVAKd3kjabMGiF5hsTvxDg9jrzkTzdZNzypWkQkPIRIRYV0NmjbbZhqW8jui+zaQjPlKJGMb8Sgx9+zDmlZGO+NHb3CR7mpDJGKGNjxIxWknHI5TULUIKPfVtXaTdJeQuuJhkPEIy6CVY/wrWypmYTDZ0ZhtWTzm5C2fSvfZR3C4ncZsVW8180tEQ3pbtSJkmlYhiKaslGfQiQ0MILQpmO6HtL2FwedC8naQiw5jsuYh0xho60d8O8QiRPa9jySlEi4ex2B3Y684gp6CYYMceSPggGEM31MbHzl5GMhZm7SMNVE2cgMNsZHh4kMBwB8H2/YhKE3lVMxFaFE+Om6qaZZx30nSeX7Oezbv2kpvykjd5Mjqzg92dO7ARQx/po9Rsx7NoATd971vEYjFu/ceDDEeNJMPDXHrjd/j1XY9QftqFxBMJenr7GRzYx4zKQr70758iGo3yf7fcRzwUQKQ1Tj99Cf6GN/D5fJSVlfGtaz/3luvUTU+wcOo0hCtNe6wAMXkJNk8xgY5GZk+ZwPnnn8/5559/kAOiEQ7lkjT3bS5be0eN0EY4XnPQRzOU/U5xjnXVi2L8ocT9GBjxs2x05GEWydGW+NFYnB9p2P54iP6R0pdSEgwG6evrO+Snq6uLzVu2EU+lCfuHSWpJmrauJye/CC2ZAqHD7/NhsbtYdO6pvLF6Jc4p00lrMSbOWEy0YTUz58ylZP45bNu9l0Kzgy3P3k3pWdegk0mKcuy8ef/rnL34EnICToLhKOG969CZbSTDQQwWJ5idGPPKMoZlRgvoTaTjYYaCCdrb2rAYSzD0+5BaDG2wA/eEaZScdhFNqx4CnYF0IorO4UHzdqEN7seYU4zBmY+5vA6Du4jghsdIB/sIWu2se3U1qbTERAqDuxBTfgWxRAST2YpmMJPwtyJTGraqeSR69qAN7Ednc5PsHwLILCsLeRFRP5WzFhEYHqagwMO0uYtobthFQOgx5hSjJfait7sxFU5AF4uRigQwuvKRMT+RpvXoTVbspdXYiydhd+cwt6YEi81Jp7WI3IlL2TjcCzkVyLwygjtXkvIPkGNMsHDBTNqSOfR3t6PXIqRjAUqnLWBg7xb0BjP2KYsxWe1IbzvhljDuWUsZbt6Kc+J0km0b0YWLiDa+yeL580CLYLZFET31mLUI9/7p/436om9paeGWux4gz5xmy/ZXMOjDBPpMzJg7H91AMxdffDEXX3wx69at47EXVuGL+ti9YTsep4WC2tnMnVFLjtuNr37tqJiOnfsFeHbV6+zfugKD04PdP0D1jGr+9KufUFhYSDgc5rnV63DVlJJfXEYsFCA4pl4dzoK7paUla0DWepAB2die94GMFdCROehDuWx9N7xfc9D/yqoXxfhCDcu/SzZs3MiXvvMjTMXVWPQwa/4idAPNo8NeRxqCO9KQ2ZF2czsSUkp8Ph99fX20tbXxuz/fiT6vAi0ew9vXhXf/HvJzXQwMDNDX14der6eoqOiQH4PBwAuvb2XamZdicbhYu3E7gfbdiHiY3HnnkQp5mT+tihf/8VuWXnEtr63bzGD7HpLxKBMmVJIcbKelvRv7zLNJCgOFeS5a1jzJ5E9+i7QEd9pP40sPUbPkYsyectpb9tK9/lmQKXInz6OoZha7/3kbtrplGFwFpPy9pL2d1CxcQl93J9GO3YSG+jFPmIMW6COvZj6hfW9SOfMk9q58mNqlnyRiLyQWCjKw/gkqZ5xMZ8NmSj/2ZULhMAajCe+aeygr9GCvXoTJaqNp7x4iXU0YnR7s088k4e1EHw+RattMXEuS1JkQRhPpsBcSMTBaMNhzSKeSGFwFJHqbKC0pZvpl32Cwez/p8DCR5g2U59pp7vHiLJ/K7nWvYCqagHvhJ5A6A6GG14i3bSV3/gVEuhoxC0kiHsWsB2c6wD23/AYhBF++4UbMJVPoaN6N1JuJx+IYdGCND3H3rTfj8Xj4xR9vZ3NDC96+bnLz8lk4q5azT57Hw88sp30ggC8cp7CoiLmLTmX3ru0MdbViNeqpnVaHwwifv+ITVFdXj778D9e4HGl4Njc38/eHnzqoYXtgPKPRyI4dO3jshVWH9GJ3IDt37uSWuzK7AI64tx0b92iHtg/keDSYP2xz0MdaVh91xtuw/LgQdyHEecAfAD1wp5Ty/x0p/rGKezgc5t++/F9sj+djLZ1COhbC3r+TeVMmcOPXv0BZWdlovEO9DJqamrLrSD89uo60ac2TfPNzn+QPd943Kvr+oX46Xv8nX/nsp4/Y0+7r66O/vx+LxUJRURG5ubkMBiKU1MzCkevBmZNPqLuZr159OTNmzKCoqOgdl7aNbXy0NjWw8r5b0NnzKJ53FgtmTaOwsIgNT/4NYzLCzv39mKx25p18Olanmwf+9HOqzvgUg+178fZ1k9ZiFLptJMy5oDeQZ5b4u1ooPvkTlM08GW9XK7ueuZOKhWdjsTnIqaxl48N/IDI8iMlTjkUvSYZ8uJ02hr3DTDv7Cvbv3oLMm0A8GsRutWLoa2BGRR5+nxefLoehuEBPGnMySH7tSex66WGc007NrE+XKcwDjdg8xbgnzuSMkxfS2NjI8vtvIxkOkpRpcvLyMUuNuro62trbmbDwHHoGhkin0sRbN3L+KXN57rWtSIuT/vZmFi27kHBKj2f2WYR79nHGyQtpWvMk37rmMm77x0NYqubx8orldOzeQgqB0BsxG/UUFRdRuOgS4v5BhBbFu2sNc6ZN4VtfuWZ01cGGDRv4890PEdJgT2M95aVleHJc/NcX/200zljf6yN7iY8YanV3d9PY2Mizq95Ab88hFfbxqQvPYtasWQetT3639eBoxe54xv2wieyJRJXVu0eJ+wcMIYQe2EvG5V4nsAm4Skq5+3D3HKu4NzU1ce6/f41w4Uz0Dg8ypRFr2sBpNfk8/Pe/HLESbd26jd/++S42bNuBuWQKE0uLSIR8DOx9kwmlRby5YzdpoSc4PEg44MNgMFBRXkZZWdlhe9pFRUUUFhZitWZ2mjoexjQHtvovP/9M7n7sWcpPvmh0uU//m8v54bevY8eOHaMW8IOdLexo62f2VTeQjEcZ3LeLhpcfZPb0qfT2dDNpUjUlHjdnLJzJr+64D2tZLbpUnOBQPzlTF7NowXw2bn6Tvi0voYV9oDdhzy1kWkU+n77oHH78hzsoWvpZOra9hrFqAfG+Vi44YzGB+rX88oc30NfXx+//ejfrt+/BnlvAGRdcisXuZMez/yAnNxdfLEVr016qq6tpaW15297X+197iisuOhuXy4XH48HtdqNpGl1dXdzz+HMHLfcbsWAOh8M88NSLozuDnXpyZmewkTJvaWnh9vseYzAYY+e2zRTmeTDpJe7cPITRxp6mPdTWzSTHqufy85cd0kXm2J7wsQqyetErFO+MEvcPGFkvOzeO7NErhPgegJTyF4e751jFff369ZxzxZdwL/sPdDY3mq8f78o7mF1dwf9+5+v4/f5D9q57e3sZGhrCYndhdbrQtCQ6vYGK0iKWLTmF2tpann15LeULzqaoshqJxLt91TFZuB6PIbkDxeCdphp8Ph+JRIJLP38dxWd9AXdhGf7+LrpW3MF9t/yakpKStwnTSI/U6MgjONiNTq8np6hidG3urFmz8Pv9wFsGT48//jg/+dPfiCQF/qF+5s2fz+SyooPyMnbJ3Uheq6rePh87ItpHU0ZH05scmWs+1PDzocQZ+JcFW6FQHF+UuH/AEEJcDpwnpfyP7PnVwGIp5X8e7p5jFfdXX32VT1x7AylHEYFtyzOeunQGnE4nc2bNoLKy8pC963Q6zR2PvcjUMy8HIB6NsPvlR/i/6z8/athzPOfJ3oue2tGkOSLA0pqLiA7zf9d/kcsuu+wd04PDz/WOZaTH7PF4MJlM/9Lw7fEuI9U7Vig+3Iw3cf/IWMsLIb4MfBmgsrLymNKYOnUq+Q4roQmzcNUsREpI7HieC5eexJ9+edMRhcSQ+ueoBWsiFhl1pTjC8dyy8b2wzD2aNC+77DKWLFkyuunGga5bj5Te0eS3sLDwiGm+m7we7zJSO3IpFIoPEuNB3LuAijHn5dmwtyGlvAO4AzI992P5ocLCQn7+P9fxvV/eij+lR0b8zKudxDevveaYPKcdi6vDDzpHK8AKhUKheO8YD8PyBjIGdWeREfVNwGeklPWHu+df3aGuv7+fPXv24HK5qK6uPmpBVkO3CoVC8cFEDct/wJBSJoUQ/wksJ7MU7q4jCfvx4Fh7p+OhZ65QKBSKDz4fenEHkFI+Dzx/ovOhUCgUCsUHAd2JzoBCoVAoFIrjixJ3hUKhUCjGGUrcFQqFQqEYZyhxVygUCoVinKHEXaFQKBSKcYYSd4VCoVAoxhlK3BUKhUKhGGcocVcoFAqFYpyhxF2hUCgUinGGEneFQqFQKMYZStwVCoVCoRhnKHFXKBQKhWKc8aF3+XosCCEGgP3/YjL5wOBxyM5HFVV+x44qu2NHld2xM97LboKUsuBEZ+J48ZEU9+OBEGLzePL9+36jyu/YUWV37KiyO3ZU2X24UMPyCoVCoVCMM5S4KxQKhUIxzlDifuzccaIz8CFHld+xo8ru2FFld+yosvsQoebcFQqFQqEYZ6ieu0KhUCgU4wwl7gqFQqFQjDOUuB8DQojzhBB7hBDNQojvnuj8fBAQQlQIIV4RQuwWQtQLIb6eDc8TQrwkhGjKfudmw4UQ4o/ZMtwhhJg3Jq3PZeM3CSE+d6L+0/uNEEIvhNgqhHg2ez5JCLEhW0YPCyFM2XBz9rw5e33imDS+lw3fI4T42In5J+8vQogcIcRjQohGIUSDEOJk9dwdHUKIb2br6y4hxINCCIt67sYJUkr1eRcfQA/sA6oAE7AdqDvR+TrRH6AEmJc9dgJ7gTrgV8B3s+HfBX6ZPb4AeAEQwEnAhmx4HtCS/c7NHuee6P/3PpXht4AHgGez548AV2aP/wJ8NXv8NeAv2eMrgYezx3XZ59EMTMo+p/oT/b/eh3K7G/iP7LEJyFHP3VGVWxnQCljHPG/XqOdufHxUz/3dswhollK2SCkTwEPAx09wnk44UsoeKeWW7HEQaCDz8vg4mZcv2e9PZI8/DtwjM6wHcoQQJcDHgJeklF4p5TDwEnDe+/hXTghCiHLgQuDO7LkAlgGPZaMcWHYjZfoYcFY2/seBh6SUcSllK9BM5nkdtwgh3MDpwN8ApJQJKaUP9dwdLQbAKoQwADagB/XcjQuUuL97yoCOMeed2TBFluxw3VxgA1AkpezJXuoFirLHhyvHj2r5/h64AUhnzz2AT0qZzJ6PLYfRMspe92fjfxTLbhIwAPw9O6VxpxDCjnru3hEpZRfwG6CdjKj7gTdRz924QIm74rgihHAAjwPfkFIGxl6TUkpArb08ACHERUC/lPLNE52XDyEGYB7wZynlXCBMZhh+FPXcHZqsHcLHyTSQSgE7H43Rio8EStzfPV1AxZjz8mzYRx4hhJGMsN8vpXwiG9yXHfYk+92fDT9cOX4Uy/dU4BIhRBuZaZ5lwB/IDBkbsnHGlsNoGWWvu4EhPppl1wl0Sik3ZM8fIyP26rl7Z84GWqWUA1JKDXiCzLOonrtxgBL3d88moCZrUWoiY1jy9AnO0wknO/f2N6BBSvnbMZeeBkYsjz8H/HNM+Gez1ssnAf7sMOpy4FwhRG62Z3FuNmzcIqX8npSyXEo5kczztEpK+W/AK8Dl2WgHlt1ImV6ejS+z4VdmrZonATXAxvfpb5wQpJS9QIcQYmo26CxgN+q5OxragZOEELZs/R0pO/XcjQdOtEXfh/FDxuJ2Lxmr0O+f6Px8ED7AaWSGPncA27KfC8jMya0EmoCXgbxsfAHcmi3DncCCMWl9gYxRTjPw+RP9397nclzKW9byVWReks3Ao4A5G27Jnjdnr1eNuf/72TLdA5x/ov/P+1Rmc4DN2WfvKTLW7uq5O7qy+zHQCOwC7iVj8a6eu3HwUdvPKhQKhUIxzlDD8gqFQqFQjDOUuCsUCoVCMc5Q4q5QKBQKxThDibtCoVAoFOMMJe4KhUKhUIwzlLgrFAqFQjHOUOKuUCj+ZYQQS0dc1R4hzhwhxAVjzi8RymWyQvGeoMRdoVC8X8whs7ERAFLKp6WU/+8E5kehGLcocVcojhNCiM8KIXYIIbYLIe4VQhQIIR4XQmzKfk7NxrtRCHGXEGK1EKJFCPFf2XC7EOK57P27hBBXZMPbhBD52eMFQojV2eMzhBDbsp+tQgjnYfK1VAixNpv2HiHEX4QQuuy1q4QQO7O/98sx94SEEL8TQtQLIVYKIQqy4auFEAuyx/nZ/fAP/L1FQoh12Ty9IYSYmt2q+Sbgimx+rxBCXCOEuCV7z0QhxKps+a0UQlRmw/8hhPhjNp0WIcTlB/6eQqE4GCXuCsVxQAgxHfgBsExKORv4OhnnL7+TUi4ELiPrqz1LLRkf4ouAH2Wd7pwHdEspZ0spZwAvvsPPfge4Tko5B1gCRI8QdxFwPVAHTAY+KYQoBX5JxlHNHGChEGLEd7cd2CylnA6sAX50FMUwQiOwRGa8tP0Q+LmUMpE9flhKOUdK+fAB9/wJuFtKOQu4H/jjmGslZLY3vghQPX2F4igwvHMUhUJxFCwDHpVSDgJIKb1CiLOBuoxPDgBcWZe4AM9JKeNAXAjRT8bf+E7g5mwP+lkp5avv8JuvA78VQtwPPCGl7DxC3I1SyhYAIcSDZMRSA1ZLKQey4fcDp5PZnz0NjAjwfWQ8hh0tbuBuIUQNGX8DxqO452Tgk9nje4Ffjbn2lJQyDewWQhQddKdCoTgI1XNXKN47dMBJ2Z7qHCllmZQylL0WHxMvBRiklHvJuCvdCfxUCPHD7PUkb9VVy8hN2fnq/wCswOtCiNoj5OVAJxLv1qnESPxD5uUAfgK8kh19uPgI8Y6WsWUlDhtLoVCMosRdoTg+rAI+JYTwAAgh8oAVZIbCyYbNOVIC2WHyiJTyPuDXZIQeoA2Ynz2+bEz8yVLKnVLKX5JxRXwkcV+UdVOsA64AXiPj2euM7Ny5HriKzBA8ZN4NI/Pbn8nGPzAvh5v/dvOWP+9rxoQHgUPaBQBvkHF3C/BvwDuNWigUiiOgxF2hOA5IKeuBnwFrhBDbgd8C/wUsyBqJ7Qa+8g7JzAQ2CiG2kZnj/mk2/MfAH4QQm8n08kf4RtYQbgeZIfYXjpD2JuAWoAFoBZ6UGT/m3yXjv3s78KaUcsR3d5hMg2AXmSmHm7LhvwG+KoTYCuQf5rd+BfwiG2fs1N8rZKYpto0YC47heuDz2f9yNRmbBYVCcYwol68KxThHCLEU+I6U8qJ3cU9ISul455gKheKDiOq5KxQKhUIxzlA9d4VinCCEmEnG0nwscSnl4hORH4VCceJQ4q5QKBQKxThDDcsrFAqFQjHOUOKuUCgUCsU4Q4m7QqFQKBTjDCXuCoVCoVCMM/4/TSRw+blTYTEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Figure: Tract-Level Sampling Bias\")\n",
    "ax = plot_sg_vs_census(panel,'census_population','number_devices_residing', alpha = 0.5, text_col='census_tract', num_to_label=9,labelfontsize=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WNj5fehYTIfB"
   },
   "source": [
    "Figure Explanation:\n",
    "* Above is a scatter plot of `census_population` on the x-axis and  `number_devices_residing` (in panel sample) on the y-axis.\n",
    "* There is one point for each census tract in the panel. The line shows the expected relationship for a true random sample. \n",
    "* Points above the line are over-indexed groups in the panel. Points below the line are under-indexed groups in the panel. "
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "Copy of Quantifying Sampling Bias in SafeGraph Patterns-BlogPost-2019-10.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}