replication_report.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [ADA-2020]: P2-milestone\n",
    "### Paper replication exercise\n",
    "\n",
    "**Objective**: Replicate the results of a given paper\n",
    "\n",
    "**Study paper**: Eunjoon Cho, Seth A. Myers and Jure Leskovec, *Friendship and Mobility: User Movement in Location-Based Social Network*, Standford University, 2011\n",
    "\n",
    "**Study Figure**: Figure 2a),  distribution of home distances between friends\n",
    "\n",
    "**Due date**: 6 November 2020, EOD\n",
    "\n",
    "**Author**: Ivan-Daniel Sievering, SCIPER n°270318, ivan.sievering@epfl.ch\n",
    "\n",
    "**Method**: The method that will be used here is the following:\n",
    "\n",
    "0. Apply some preprocessing to clean data and remove useless information;\n",
    "1. Convert long-lat check-ins' positions to x-y position (Mercator);\n",
    "2. Assign to each check-in its square in the map;\n",
    "3. Find for each user the square that has more check-ins;\n",
    "4. Mean all the check-ins in that square to estimate home position;\n",
    "5. Convert homes position from x-y to long-lat;\n",
    "6. Compute the distance (Haversine) between two friends homes;\n",
    "7. Estimate the probability distribution of the distance with a Gaussian Kernel Density Estimator;\n",
    "8. Plot the result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Extraction of the data and preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Librairies\n",
    "import pandas as pd\n",
    "\n",
    "#Constants\n",
    "PATH = \"data/\"\n",
    "EXTENSION = \".txt.gz\"\n",
    "SEPARATOR = \"\\t\"\n",
    "CHECKIN_COLUMNS_NAME = [\"user\", \"check-in time\", \"latitude\", \"longitude\", \"location id\"]\n",
    "EDGE_COLUMNS_NAME = [\"friend1\", \"friend2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Extract dataframe from .txt.gz (tabulated separated without header)\n",
    "loc_brightkite_edges = pd.read_csv(PATH+\"loc-brightkite_edges\"+EXTENSION, sep=SEPARATOR, header = None, names = EDGE_COLUMNS_NAME)\n",
    "loc_brightkite_totalCheckins = pd.read_csv(PATH+\"loc-brightkite_totalCheckins\"+EXTENSION, sep=SEPARATOR, header = None, names = CHECKIN_COLUMNS_NAME)\n",
    "loc_gowalla_edges = pd.read_csv(PATH+\"loc-gowalla_edges\"+EXTENSION, sep=SEPARATOR, header = None, names = EDGE_COLUMNS_NAME)\n",
    "loc_gowalla_totalCheckins = pd.read_csv(PATH+\"loc-gowalla_totalCheckins\"+EXTENSION, sep=SEPARATOR, header = None, names = CHECKIN_COLUMNS_NAME)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>user</th>\n",
       "      <th>check-in time</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>location id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>423638</th>\n",
       "      <td>1017</td>\n",
       "      <td>2008-06-20T19:56:57Z</td>\n",
       "      <td>39.746180</td>\n",
       "      <td>-104.995574</td>\n",
       "      <td>3e2568f9bafa6dd1713c19baa1b6edf63386177d</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1934953</th>\n",
       "      <td>7603</td>\n",
       "      <td>2009-08-15T08:14:02Z</td>\n",
       "      <td>28.302351</td>\n",
       "      <td>-16.512186</td>\n",
       "      <td>f99be920896a11de9df3003048c10834</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1464215</th>\n",
       "      <td>5694</td>\n",
       "      <td>2008-10-12T21:38:41Z</td>\n",
       "      <td>42.135650</td>\n",
       "      <td>-71.970074</td>\n",
       "      <td>ed6e8820a22411dd837157cf69a71e47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>451774</th>\n",
       "      <td>1094</td>\n",
       "      <td>2008-05-22T07:50:45Z</td>\n",
       "      <td>48.900000</td>\n",
       "      <td>2.300000</td>\n",
       "      <td>d993133ea22411dd8903d7290d7db445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2786481</th>\n",
       "      <td>11619</td>\n",
       "      <td>2009-05-05T10:06:19Z</td>\n",
       "      <td>3.171760</td>\n",
       "      <td>101.603820</td>\n",
       "      <td>9daf7e589cdd11dd9e01003048c0801e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1348408</th>\n",
       "      <td>4681</td>\n",
       "      <td>2009-01-31T21:15:24Z</td>\n",
       "      <td>37.158997</td>\n",
       "      <td>-93.272385</td>\n",
       "      <td>834e8402ea4c11dd8b6c003048c10834</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2884043</th>\n",
       "      <td>12376</td>\n",
       "      <td>2008-06-09T18:26:54Z</td>\n",
       "      <td>53.342031</td>\n",
       "      <td>-6.254642</td>\n",
       "      <td>0947a15f032daf2154d4f9512013a0cc98458c33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3801160</th>\n",
       "      <td>23761</td>\n",
       "      <td>2008-10-04T23:31:47Z</td>\n",
       "      <td>35.184655</td>\n",
       "      <td>-106.598454</td>\n",
       "      <td>e252512a91d911ddba81003048c10834</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>220763</th>\n",
       "      <td>300</td>\n",
       "      <td>2008-11-23T17:36:20Z</td>\n",
       "      <td>33.924269</td>\n",
       "      <td>-84.378538</td>\n",
       "      <td>eb5b200ca22411ddb5b6e79ef862e7eb</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3658232</th>\n",
       "      <td>21206</td>\n",
       "      <td>2008-06-26T16:27:42Z</td>\n",
       "      <td>52.215244</td>\n",
       "      <td>5.191089</td>\n",
       "      <td>93928ad0c7605633b9168b7435aff110e8d6df72</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          user         check-in time   latitude   longitude  \\\n",
       "423638    1017  2008-06-20T19:56:57Z  39.746180 -104.995574   \n",
       "1934953   7603  2009-08-15T08:14:02Z  28.302351  -16.512186   \n",
       "1464215   5694  2008-10-12T21:38:41Z  42.135650  -71.970074   \n",
       "451774    1094  2008-05-22T07:50:45Z  48.900000    2.300000   \n",
       "2786481  11619  2009-05-05T10:06:19Z   3.171760  101.603820   \n",
       "1348408   4681  2009-01-31T21:15:24Z  37.158997  -93.272385   \n",
       "2884043  12376  2008-06-09T18:26:54Z  53.342031   -6.254642   \n",
       "3801160  23761  2008-10-04T23:31:47Z  35.184655 -106.598454   \n",
       "220763     300  2008-11-23T17:36:20Z  33.924269  -84.378538   \n",
       "3658232  21206  2008-06-26T16:27:42Z  52.215244    5.191089   \n",
       "\n",
       "                                      location id  \n",
       "423638   3e2568f9bafa6dd1713c19baa1b6edf63386177d  \n",
       "1934953          f99be920896a11de9df3003048c10834  \n",
       "1464215          ed6e8820a22411dd837157cf69a71e47  \n",
       "451774           d993133ea22411dd8903d7290d7db445  \n",
       "2786481          9daf7e589cdd11dd9e01003048c0801e  \n",
       "1348408          834e8402ea4c11dd8b6c003048c10834  \n",
       "2884043  0947a15f032daf2154d4f9512013a0cc98458c33  \n",
       "3801160          e252512a91d911ddba81003048c10834  \n",
       "220763           eb5b200ca22411ddb5b6e79ef862e7eb  \n",
       "3658232  93928ad0c7605633b9168b7435aff110e8d6df72  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Sample dataframe to analyze and verify structure\n",
    "#loc_brightkite_edges.sample(10)\n",
    "loc_brightkite_totalCheckins.sample(10)\n",
    "#loc_gowalla_edges.sample(10)\n",
    "#loc_gowalla_totalCheckins.sample(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Remove from checkin the useless information (location ID and check-in time)\n",
    "TO_REMOVE = [\"check-in time\", \"location id\"]\n",
    "bright_checkins = loc_brightkite_totalCheckins.drop(columns=TO_REMOVE)\n",
    "gowalla_checkins = loc_gowalla_totalCheckins.drop(columns=TO_REMOVE)\n",
    "\n",
    "#Remove from checkin some invalid points (null data and incoherent latitudes)\n",
    "bright_0_lat = bright_checkins[bright_checkins[\"latitude\"] == 0]\n",
    "bright_0_lat_long = bright_0_lat[bright_0_lat[\"longitude\"] == 0]\n",
    "bright_checkins = bright_checkins.drop(bright_0_lat_long.index)\n",
    "bright_big_lat = bright_checkins[bright_checkins[\"latitude\"]>90]\n",
    "bright_checkins = bright_checkins.drop(bright_big_lat.index)\n",
    "bright_small_lat = bright_checkins[bright_checkins[\"latitude\"]<-90]\n",
    "bright_checkins = bright_checkins.drop(bright_small_lat.index)\n",
    "\n",
    "gowalla_0_lat = gowalla_checkins[gowalla_checkins[\"latitude\"] == 0]\n",
    "gowalla_0_lat_long = gowalla_0_lat[gowalla_0_lat[\"longitude\"] == 0]\n",
    "gowalla_checkins = gowalla_checkins.drop(gowalla_0_lat_long.index)\n",
    "gowalla_big_lat = gowalla_checkins[gowalla_checkins[\"latitude\"]>90]\n",
    "gowalla_checkins = gowalla_checkins.drop(gowalla_big_lat.index)\n",
    "gowalla_small_lat = gowalla_checkins[gowalla_checkins[\"latitude\"]<-90]\n",
    "gowalla_checkins = gowalla_checkins.drop(bright_small_lat.index)\n",
    "\n",
    "#Rename friendships data\n",
    "bright_edges = loc_brightkite_edges\n",
    "gowalla_edges = loc_gowalla_edges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Compute home position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a) Determine the 25km square with most check-ins\n",
    "Data are projected to a Mercator plan, where squares are easy to define. The square with most check-ins is extracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Librairies\n",
    "from pyproj import Transformer\n",
    "import numpy as np\n",
    "\n",
    "#Constants\n",
    "SQUARE_SIZE = 25000\n",
    "\n",
    "#Define the transformation from long-lat to x-y\n",
    "trans = Transformer.from_crs(\n",
    "    \"epsg:4326\", #long lat notation\n",
    "    \"epsg:3857\", #pseudo mercator https://epsg.io/3857\n",
    "    always_xy=True,\n",
    ")\n",
    "\n",
    "#Transform data long-lat -> x-y\n",
    "bright_checkins[\"x\"], bright_checkins[\"y\"] = trans.transform(bright_checkins[\"longitude\"].values, \n",
    "                                                             bright_checkins[\"latitude\"].values)\n",
    "gowalla_checkins[\"x\"], gowalla_checkins[\"y\"] = trans.transform(gowalla_checkins[\"longitude\"].values, \n",
    "                                                             gowalla_checkins[\"latitude\"].values)\n",
    "\n",
    "#Some checkins at extreme values cannot be converted due to the projection -> remove them\n",
    "bright_checkins = bright_checkins.replace([np.inf, -np.inf], np.nan).dropna(axis=0)\n",
    "gowalla_checkins = gowalla_checkins.replace([np.inf, -np.inf], np.nan).dropna(axis=0)\n",
    "\n",
    "#find square position to which the point belongs (take the divisor of the euclidian division) \n",
    "#(!i_x/i_y may be negative because the origin is at the center of Earth (Null Island))\n",
    "bright_checkins[\"i_x\"] = (bright_checkins[\"x\"]/SQUARE_SIZE).astype(int)\n",
    "bright_checkins[\"i_y\"] = (bright_checkins[\"y\"]/SQUARE_SIZE).astype(int)\n",
    "gowalla_checkins[\"i_x\"] = (gowalla_checkins[\"x\"]/SQUARE_SIZE).astype(int)\n",
    "gowalla_checkins[\"i_y\"] = (gowalla_checkins[\"y\"]/SQUARE_SIZE).astype(int)\n",
    "\n",
    "#Now we can drop x, y, latitude and longitude\n",
    "TO_REMOVE = [\"latitude\", \"longitude\",\"x\"] #keep one column for the count\n",
    "bright_checkins_light = bright_checkins.drop(columns=TO_REMOVE)\n",
    "gowalla_checkins_light = gowalla_checkins.drop(columns=TO_REMOVE)\n",
    "bright_checkins_light = bright_checkins_light.rename(columns={\"y\":\"count\"})\n",
    "gowalla_checkins_light = gowalla_checkins_light.rename(columns={\"y\":\"count\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For each user find the most used square\n",
    "\n",
    "#Group by user and square position and count the nb of occurence\n",
    "bright_checkins_zone_count = bright_checkins_light.groupby([\"user\",\"i_x\",\"i_y\"]).count().reset_index()\n",
    "#create a dataframe with the max value\n",
    "bright_checkins_zone_max = bright_checkins_zone_count.groupby([\"user\"]).max() \\\n",
    "                           .rename(columns={\"count\":\"count_max\"}).drop(columns=[\"i_x\", \"i_y\"])\n",
    "#merge both dataframes\n",
    "bright_checkins_merge = pd.merge(bright_checkins_zone_max, bright_checkins_zone_count, left_on=\"user\",right_on=\"user\")\n",
    "#keep only the lines where the count is equal to the maximum\n",
    "bright_checkins_home = bright_checkins_merge[(bright_checkins_merge[\"count_max\"] == bright_checkins_merge[\"count\"])]\n",
    "#There may be still some duplicates if tie (keep first)\n",
    "bright_checkins_home = bright_checkins_home.drop_duplicates(subset=[\"user\"],keep=\"first\")\n",
    "#by rehandling the dataframe we have for each user it's prefered square\n",
    "bright_square_home = bright_checkins_home.drop(columns=[\"count_max\",\"count\"]).set_index(\"user\")\n",
    "\n",
    "#Group by user and square position and count the nb of occurence\n",
    "gowalla_checkins_zone_count = gowalla_checkins_light.groupby([\"user\",\"i_x\",\"i_y\"]).count().reset_index()\n",
    "#create a dataframe with the max value\n",
    "gowalla_checkins_zone_max = gowalla_checkins_zone_count.groupby([\"user\"]).max() \\\n",
    "                           .rename(columns={\"count\":\"count_max\"}).drop(columns=[\"i_x\", \"i_y\"])\n",
    "#merge both dataframes\n",
    "gowalla_checkins_merge = pd.merge(gowalla_checkins_zone_max, gowalla_checkins_zone_count, left_on=\"user\",right_on=\"user\")\n",
    "#keep only the nb of checkins of the square is equal to the maximum\n",
    "gowalla_checkins_home = gowalla_checkins_merge[(gowalla_checkins_merge[\"count_max\"] == gowalla_checkins_merge[\"count\"])]\n",
    "#There may be still some duplicates if tie (keep first)\n",
    "gowalla_checkins_home = gowalla_checkins_home.drop_duplicates(subset=[\"user\"],keep=\"first\")\n",
    "#by rehandling the dataframe we have for each user it's prefered square\n",
    "gowalla_square_home = gowalla_checkins_home.drop(columns=[\"count_max\",\"count\"]).set_index(\"user\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b) Determine the position of the home\n",
    "Remove all check-ins that are not in most used square, mean the position of the remainings check-ins of each user and then reproject to the long-lat space (no loss in precision, the transformation is reversible + we could not have used mean consistently in long-lat space due to the nature of the unit) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For each row add the home square of the user\n",
    "bright_checkins_full = pd.merge(bright_checkins, bright_square_home, left_on=\"user\", \n",
    "                                right_index=True, how=\"outer\", suffixes=('', '_h'))\n",
    "gowalla_checkins_full = pd.merge(gowalla_checkins, gowalla_square_home, left_on=\"user\", \n",
    "                                right_index=True, how=\"outer\", suffixes=('', '_h'))\n",
    "\n",
    "#Remove checkins not in home square\n",
    "bright_home_checkins = bright_checkins_full[(bright_checkins_full[\"i_x\"] == bright_checkins_full[\"i_x_h\"]) \n",
    "                                            & (bright_checkins_full[\"i_y\"] == bright_checkins_full[\"i_y_h\"])]\n",
    "gowalla_home_checkins = gowalla_checkins_full[(gowalla_checkins_full[\"i_x\"] == gowalla_checkins_full[\"i_x_h\"]) \n",
    "                                            & (gowalla_checkins_full[\"i_y\"] == gowalla_checkins_full[\"i_y_h\"])]\n",
    "\n",
    "#Remove useless columns\n",
    "TO_REMOVE = [\"i_x\", \"i_y\", \"i_x_h\", \"i_y_h\", \"longitude\", \"latitude\"]\n",
    "bright_home_checkins_light = bright_home_checkins.drop(columns=TO_REMOVE)\n",
    "gowalla_home_checkins_light = gowalla_home_checkins.drop(columns=TO_REMOVE)\n",
    "\n",
    "#Mean the position for each user\n",
    "bright_home_pos = bright_home_checkins_light.groupby([\"user\"]).mean()\n",
    "gowalla_home_pos = gowalla_home_checkins_light.groupby([\"user\"]).mean()\n",
    "\n",
    "#Define the transformation from x-y to long-lat (return to long lat for more precise distance mesure) (linear trans -> ok)\n",
    "trans = Transformer.from_crs(\n",
    "    \"epsg:3857\", #pseudo mercator https://epsg.io/3857\n",
    "    \"epsg:4326\", #long lat notation\n",
    "    always_xy=True,\n",
    ")\n",
    "\n",
    "#Transform data from x-y to long-lat\n",
    "long, lat = trans.transform(bright_home_pos[\"x\"].values, bright_home_pos[\"y\"].values)\n",
    "bright_home_pos[\"lat\"] = lat\n",
    "bright_home_pos[\"long\"] = long\n",
    "long, lat = trans.transform(gowalla_home_pos[\"x\"].values, gowalla_home_pos[\"y\"].values)\n",
    "gowalla_home_pos[\"lat\"] = lat\n",
    "gowalla_home_pos[\"long\"] = long\n",
    "\n",
    "#Now we can drop x and y\n",
    "TO_REMOVE = [\"x\", \"y\"]\n",
    "bright_home_pos_light = bright_home_pos.drop(columns=TO_REMOVE)\n",
    "gowalla_home_pos_light = gowalla_home_pos.drop(columns=TO_REMOVE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lat</th>\n",
       "      <th>long</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>user</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>30.263758</td>\n",
       "      <td>-97.751843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>45.202484</td>\n",
       "      <td>-0.746831</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>37.858417</td>\n",
       "      <td>-122.268497</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>37.775651</td>\n",
       "      <td>-122.440133</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>30.288267</td>\n",
       "      <td>-97.736698</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>196544</th>\n",
       "      <td>-25.433409</td>\n",
       "      <td>-49.281533</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>196561</th>\n",
       "      <td>37.528673</td>\n",
       "      <td>-122.006051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>196577</th>\n",
       "      <td>51.748050</td>\n",
       "      <td>-0.478077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>196578</th>\n",
       "      <td>51.745909</td>\n",
       "      <td>-0.495867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>196585</th>\n",
       "      <td>50.105516</td>\n",
       "      <td>8.571525</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>107069 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "              lat        long\n",
       "user                         \n",
       "0       30.263758  -97.751843\n",
       "1       45.202484   -0.746831\n",
       "2       37.858417 -122.268497\n",
       "4       37.775651 -122.440133\n",
       "5       30.288267  -97.736698\n",
       "...           ...         ...\n",
       "196544 -25.433409  -49.281533\n",
       "196561  37.528673 -122.006051\n",
       "196577  51.748050   -0.478077\n",
       "196578  51.745909   -0.495867\n",
       "196585  50.105516    8.571525\n",
       "\n",
       "[107069 rows x 2 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gowalla_home_pos_light"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Compute distance between friends homes position\n",
    "Add to the list of friends their home position and compute the Haversine distance (https://en.wikipedia.org/wiki/Haversine_formula) between their home. (There is no need to check if pairs of friends appear one or two times because this will not affect the probability distribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Library\n",
    "import numpy as np\n",
    "\n",
    "#Constants\n",
    "R = 6371000 #radius of the earth, in meters\n",
    "\n",
    "#Create a structure with the friend relationships and their homes position\n",
    "bright_friend_homes = pd.merge(loc_brightkite_edges, bright_home_pos_light, left_on=\"friend1\", right_index=True, how=\"outer\")\n",
    "bright_friend_homes = bright_friend_homes.rename(columns={\"lat\":\"lat1\",\"long\":\"long1\"})\n",
    "bright_friend_homes = pd.merge(bright_friend_homes, bright_home_pos_light, left_on=\"friend2\", right_index=True, how=\"outer\")\n",
    "bright_friend_homes = bright_friend_homes.rename(columns={\"lat\":\"lat2\",\"long\":\"long2\"})\n",
    "gowalla_friend_homes = pd.merge(loc_gowalla_edges, gowalla_home_pos_light, left_on=\"friend1\", right_index=True, how=\"outer\")\n",
    "gowalla_friend_homes = gowalla_friend_homes.rename(columns={\"lat\":\"lat1\",\"long\":\"long1\"})\n",
    "gowalla_friend_homes = pd.merge(gowalla_friend_homes, gowalla_home_pos_light, left_on=\"friend2\", right_index=True, how=\"outer\")\n",
    "gowalla_friend_homes = gowalla_friend_homes.rename(columns={\"lat\":\"lat2\",\"long\":\"long2\"})\n",
    "\n",
    "#Because we removed checkout during the process their may exist user without home position so we have to remove them\n",
    "bright_friend_homes = bright_friend_homes.dropna()\n",
    "gowalla_friend_homes = gowalla_friend_homes.dropna()\n",
    "\n",
    "#Remove useless columns\n",
    "TO_REMOVE = [\"friend1\", \"friend2\"]\n",
    "bright_friend_homes_light = bright_friend_homes.drop(columns=TO_REMOVE)\n",
    "gowalla_friend_homes_light = gowalla_friend_homes.drop(columns=TO_REMOVE)\n",
    "\n",
    "#Compute the distance using the Haversine distance\n",
    "#Transform long/lat to radians, create delta values\n",
    "lat1_r = np.radians(bright_friend_homes_light[\"lat1\"])\n",
    "lat2_r = np.radians(bright_friend_homes_light[\"lat2\"])\n",
    "lon1_r = np.radians(bright_friend_homes_light[\"long1\"])\n",
    "lon2_r = np.radians(bright_friend_homes_light[\"long2\"])\n",
    "d_lat = lat2_r - lat1_r\n",
    "d_lon = lon2_r - lon1_r\n",
    "\n",
    "#Formula\n",
    "a = np.sin(d_lat / 2)**2 + np.cos(lat1_r) * np.cos(lat2_r) * np.sin(d_lon / 2)**2\n",
    "c = 2 * np.arcsin(np.sqrt(a))\n",
    "dist = R*c\n",
    "\n",
    "#Convert to km\n",
    "dist_km_bright = dist/1000\n",
    "\n",
    "#Transform long/lat to radians, create delta values\n",
    "lat1_r = np.radians(gowalla_friend_homes_light[\"lat1\"])\n",
    "lat2_r = np.radians(gowalla_friend_homes_light[\"lat2\"])\n",
    "lon1_r = np.radians(gowalla_friend_homes_light[\"long1\"])\n",
    "lon2_r = np.radians(gowalla_friend_homes_light[\"long2\"])\n",
    "d_lat = lat2_r - lat1_r\n",
    "d_lon = lon2_r - lon1_r\n",
    "\n",
    "#Formula\n",
    "a = np.sin(d_lat / 2)**2 + np.cos(lat1_r) * np.cos(lat2_r) * np.sin(d_lon / 2)**2\n",
    "c = 2 * np.arcsin(np.sqrt(a))\n",
    "dist = R*c\n",
    "\n",
    "#Convert to km\n",
    "dist_km_gowalla = dist/1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Plot probabilty distribution\n",
    "The probability distribution is estimated using a gaussian kernel density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAFQCAYAAACCgGiMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABJa0lEQVR4nO3deXiU5dX48e8hJCQCjQurIgaKGhYrmygtKKgVRHBp7WujFRK1iC1Y+PWtVVNrawWF1wIt2lorClGJu9a1YgWqqJRFUEEUNKBGwl4RlCXA+f1xP5M8M5lJJslMZjI5n+uaa2ae9Z5nljP3LqqKMcYYE2vNEp0AY4wxqckCjDHGmLiwAGOMMSYuLMAYY4yJCwswxhhj4sICjDHGmLho8gFGRHJEREWkuff8ZREZE4fz3Csit8T6uFGc9zoR2SIie0TkmCj3GSwiH8U7bfEmIteLyJ01bPM9EVnvXZ+LGyhptRbP90RE5ojI7VFst1FEzo2wLmk+M973uVscjz9EREqrWV/xXQ/dVkTWiMiQeKUt6ahqo7gBG4G9wB5gMzAHaBWD4+YACjSPYVrzgcVJcM3SvWt2ag2vfY/v9m6i0x2Sxl8Bq4HdwAbgV7XYNxMoBdpVs81rwC8S9NqS4vp736Xbo9huI3Buoj8TUaRTgW4R1i3y1p8asvxZb/mQKI4/BCiNMi0RtwV+Bzyc6OsVz1tjy8GMUtVWQG+gD3BTYpOT9NrjfmTX1LDdkarayrudGs8EiUhabXcBRgNHAcOB8SLy42h2VNV9wMve/pGcQITrI05DfEeqvf6B3LWJmXX4PhNezv4MYFvCUtSIVfedbmwBBgBV3Qy8ggs0AIjIGSLyloh8KSLv+rOhIrJIRO4QkaUisktE/iEiR4c7trftNb7nPxWRtSKyW0Q+EJG+3vIbReQT3/JLvOXdgXuBgV6xy5fe8qBiCO+4H4vIThF5TkSO9a1TERnnFd38V0TuERGJkN4WIjJTRDZ5t5nespOAQJHFlyKyINrrGyZb31dEVnqv9QkReSzwWkQkX0QWh+xfUUThve6/ishLIvI1MFREjhWRp0Rkm4hsEJHrI6VFVaep6juqelBVPwL+AXzPd64/icjnIvKViKwQkcEhh1gEXBDhdX4CdAWe996rFt77P1lE3gS+AbqKyHdFZJn32VkmIt/1HWORiNzuffb2iMjzInKMiDzipWmZiOREcdn96RoiIqUi8msR2Qw8GOY9iXgNReR3IvK4iBR579kaEenvW99HRN7x1j2G+xMSWNdGRF7wvkc7ReQNCQ6yvUXkPe9aPCYimf40+46zUURu8r4b/xWRBwPbhnm93xaRBSKyQ0S2e9fuyJBj/W+483rrfyUiZd7n/6ooLvEjwGVS+cOYBzwDHPAdM+z3KiTdN3vp3SgiV/iWRyxy9LY9V0SGAzd76dgj7jfrRyKyImT7X4rIsxGOVe3vmrjv6mZv3esi0jMkjfeKyKve5+DfInKCb32ut26niHwkIv8Tsm/QdzrilU50FiraG77sOdAJeB/4k/f8OGAHMAIXNL/vPW+rldniL4BeQEvgKbysKSFFZN6213iPf+Ttdxrun3Q34ATfumO9810GfA109NblE1JEhq8YAjgb2A70BVoAs4DXfdsq8AJwJNAZ989qeITrchuwBGgHtAXeAv4Q7rWF2TfsenzZeiAD+BT4Ba7I7Qe4L+Lt1bzWiiIK73XvwgWFZsARwArgt96xuwIlwLAoPgMCrATG+Zb9BDgGaA78Eld8mulb3xfYGc3nyvf+fwb09I7ZHvgvcKX3PM97foxv+4+BbwPZwAe4f8jnetsXAQ/W4fofBKZ6n4+skPekWXXXEFf0sg/3fUgD7gCWhLyfk7z381Kg3Pd+3oH7g5Tu3QYD4rtWS3Gf+6OBtYH3gpCiIG/b1cDx3rZvEqEYDve9+r73WtsCrwMzQ44V6bzDgS1UfrfnUXMR2TXAfOB8b9lSYCCuOHVIFN+rwPsz3UvzWbjv/8lhvuvhrsu5vvfpYd+6FsBOoLtv2Urgh9W8lrC/a976q4DW3nFnAqtCfo92A2d66/+E9z32jvU5UID7DPfF/V71jPCdzgyXPlVtdAFmj3dRFFd2fqS37tfAQyHbvwKM8b0Rd/rW9cD9SKZRfYB5hSjL54FVwEXe43yqDzCzgWm+da1wX/Ic77kCg3zrHwdujHDeT4ARvufDgI3e46DXFmbfwPovfbf/JfjH7EzvQyy+/RZTuwBT5Ft3OvBZyPY3EeFHOGS73wPvAi2q2ea/+MrXgROBQzV8rkIDzG2+51cCS0P2eRvI921f6Fv3R+Bl3/NR+L7Ytbj+BwgOlP73pNpriPvh+lfI532v7/3cFPJ+vuV7P2/D5RKr/EB71+onvufTgHtD0+fb1v9HYATwSZTfpYuBlVGe9wGCv9snEV2A+QlQDJwMrPPW+QNMdd+rIbgA0zLkO3pLmO96uOsSNsB4y/4KTPYe98R9nsN+3qnmdy3Mtkd61yXbl8ZHfetbAYdwfwguA94I2f9vwK2+fYvCpSn01tiKyC5W1da4Ny0XaOMtPwH4kZet/1JcsdQgoKNv3899jz/F/TtrQ/WOx33QqhCR0SKyyne+XlEcL+BYLw0AqOoeXI7rON82m32Pv8F9AGo8lvf42AjbRtJGVY/0bneFOf4X6n2yPJ9TO/7tTwCODXmvbsblFCISkfG4cvMLVHW/b/kvxRVh7vKOlU3w+9Aa92+rrukNvb54z/3v1Rbf471hnkd67wLCXf9t6uqQwonmGoZ+fjLF1eWEez/9r+//cDmy+SJSIiI3hpw72s8lVP3Ohf1cikg7EXlURL4Qka+Ah6n6XYp03mPDnCcaT+NKEiYAD4VZX9P36r+q+nU16+tqLnC5iAjuz83j/s97GGF/10QkTUTuFFeM/xUusEHwda3Y1/sN2um9hhOA00M+X1cAHSKcN6LGFmAAUNV/46Jo4Mv4OS4Hc6Tv1lJV/U1Uj/c97ozLMWyv4VSf44o+gnhllX8HxuOKSo7EFQcE6kk0dJ8Qm3BvYuB4LXHFPF/UsF+Nx8K9tk11OE4kZcBx3gc+wH8tv8YVewEgIv4PYUBocNoQ8l61VtURkRLglavfCJyjqv5y/sG43Ov/AEd578MuKt8HgO64XE9t+NMben3BXeO6vFd1TUOoWl9Dn3DvZ+eKk6ruVtVfqmpXXO7r/4nIOXV6BVW/c5E+l3fgXu93VPVbuNxF2DrHMMrCnKdGqvoNrgHIdYQPMDV9r47yvreR1keVjDDpWoLLhQwGLo+QNr9Iv2uXAxfhimqzcbllCL6uFfuKSCtc8eMm3Ofr3yGfr1aqel11aQ+nUQYYz0zg+yLSG/ePZ5SIDPMid6a4SsdOvu1/IiI9ROQIXDHAk6p6qIZz3A/8r4j0E6ebF1xa4i7wNgARKcDlYAK2AJ1EJCPCcecBBSLS26s4nAL8R1U31uL1BxQDvxGRtiLSBlcu/3AdjhPJ27is83gRaS4iFwEDfOvfBXp6ryUTl+2vzlLgK3EV2Fne+9VLRE4Lt7FXeToF+L6qloSsbo0rqtgGNBeR3wLfCtnmLNwPSV29BJwkIpd7r/8yXFHEC/U4Zn3V6hqGeBt3za73Xs8P8L2fIjLS+5wL8BXuva/pexLJz0Wkk1fxfDPwWITtWuOKv78UkeNwTdOj9TiQ7/tu31qLfW8GzorwvYvme/V7Ecnw/uiMBJ6oxbnB/U7kSNWWikXA3cBBVV1cdbcgkX7XWgP7cSUjR+C+Q6FGiMgg73fqD7jfoM9xn+2TRORKEUn3bqeJa8BUK402wKjqNtwbcYt3US7CfWC24SLwrwh+fQ/hcj2bca1mIrZc8p3jCWAyLiDsxrWVP1pVP8CVtb+N+5CcgqvEDFiAa/q6WUSq5JJU9TXgFlylXBkulxRV09swbgeWA+/hGj684y2LCVU9gKvYvxpXR/AT3Adwv7d+He6D/S9gPa5+prrjHcL9M+6N69eyHRfIsyPscjsud7dMXGubPSJyr7fuFVzwWIcrHtiHL+vuBbwRuGKHOlHVHbgfj1/ivqw3ACNVtabcb9zU4Rr69w28n/m48v3LcMVFASfi3ss9uM/3X1R1UR2TOg9XmV7i3SJ9Ln+Pq0jeBbwYkp5qqerLuD+bC3BFe1G3llTVTdX8gNf0vdqMu36bcK3Sxqnqh9Ge2xMISDtE5B3f8odwf1hryr0Etp1D1d+1Itx34gtcw5MlYfadhwvIO4F+uGIwVHU3cB7uN2mTd+xAg5NaCbQOSWkisghXmXZ/otOSCkTkP7hK1gcTnZbqiMgE4HhVvSHRaWlqRGQjrrHMvxKdlsZGRLKArUBfVV1fzXaLqOPvmojMwTU++E1d0xkN68BlaiQiZ+H61GzH/cv5DvDPhCYqCqo6K9FpMKYOrgOWVRdcGoukDzBeRdpfcBVfi1T1kQQnqSk6GVfW3QrXqu5SVS1LbJKMST1ezk9wTbUbvYQUkYnIA7hy7a2q2su3fDiuw08acL+q3ikiVwJfqurzIvKYql7W4Ak2xhhTa4mq5J+D64FbQdywDfcA5+Na6eSJSA9cr/1AxW1dW7MYY4xpYAkpIlPV16Xq+EwDgI8DTVFF5FFcy7BSXJBZRTUBUUTGAmMBWrZs2S83Nzf2CTfGmCZgxYoV21W1bX2Pk0x1MMcR3Du0FDckxp+Bu0XkAuD5SDur6n3AfQD9+/fX5cuXxzGpxhiTukQk2hERqpVMASZcz131hmMoiOoAIqOAUd26xW2uIWOMMVFKpo6WpQQPe9CJWg69oKrPq+rY7Owa+5sZY4yJs2QKMMuAE0Wkizd0wY+B5xKcJmOMMXWUkCIyESnGjYjcRtwkRbeq6mxxI+a+gmum/ICq1jQTY+hxrYjMmBRWXl5OaWkp+/ZFGmja1EZmZiadOnUiPT09LsdPyaFirJLfmNS0YcMGWrduzTHHHIOEn+TVRElV2bFjB7t376ZLly5B60Rkhar2j7Br1JKpiKzeRGSUiNy3a1dtp/8wxjQG+/bts+ASIyLCMcccE9fcYEoFGKvkNyb1WXCJnXhfy5QKMMYYY5JHSgUYKyIzxsRbWloavXv35tRTT6Vv37689dZbEbf97ne/W+PxcnJy2L696vRCixYtCjp2fn4+Tz75ZJXtrrnmGj744AMApkwJN69Y4qRUgLEiMmOMX3Ex9OoFaWnuvri4/sfMyspi1apVvPvuu9xxxx3cdNNNVbY5dMgNm1hd8KlJaICJ5P7776dHjx6ABRhjjGkQxcVQWAizZsG+fe6+sDA2QSbgq6++4qijjgJcQBg6dCiXX345p5xyCgCtWrUC4PDhw/zsZz+jZ8+ejBw5khEjRgTlRmbNmkXfvn055ZRT+PDDD9m4cSP33nsvM2bMoHfv3rzxxhtB573lllvIz8/n8OHDDBkyhOXLl3PjjTeyd+9eevfuzRVXXAHAww8/zIABA+jduzfXXnttReBrKBZgjDEpafJkmD0bhg6F9HR3P3u2W14fgR/x3NxcrrnmGm655ZaKdUuXLmXy5MkVRVYBTz/9NBs3buT999/n/vvv5+233w5a36ZNG9555x2uu+467rrrLnJychg3bhyTJk1i1apVDB48uGLbG264ga1bt/Lggw/SrFnlT/idd95Zkbt65JFHWLt2LY899hhvvvkmq1atIi0tjUceadjptJJpLLJ6s46WxpiAtWth0KDgZYMGueX1EfgRB3j77bcZPXo0q1evBmDAgAFV+pQALF68mB/96Ec0a9aMDh06MHTo0KD1P/jBDwDo168fTz/9dMRz/+EPf+D000/nvvvuqzGdr732GitWrOC0004DXGBs165dVK8xVlIqwKjq88Dz/fv3/2mi02KMSazu3WHxYpdzCVi82C2PlYEDB7J9+3a2bdsGQMuWLcNuV1OH9hYtWgCuAcHBgwcjbnfaaaexYsUKdu7cydFHH13tMVWVMWPGcMcdd1S7XTxZEZkxJiUVFsLVV8PChVBe7u6vvtotj5UPP/yQQ4cOccwxx1S73aBBg3jqqac4fPgwW7ZsYdGiRTUeu3Xr1uzevTto2fDhw7nxxhu54IILqqwDSE9Pp7y8HIBzzjmHJ598kq1btwKwc+dOPv00JqPwRy2lcjDGGBOQl+fuJ0xwxWLdu7v6l8DyugrUwYDLJcydO5e0tLRq9/nhD3/Ia6+9Rq9evTjppJM4/fTTqam166hRo7j00kv5xz/+waxZsyqW/+hHP2L37t1ceOGFvPTSS0H7jB07lu985zv07duXRx55hNtvv53zzjuPw4cPk56ezj333MMJJ5xQtxdeByk1FpmvDuan69evT3RyjDExtnbtWrrHsoyrAe3Zs4dWrVqxY8cOBgwYwJtvvkmHDh0Snayw1zRWY5GlVA7G6mCMMclq5MiRfPnllxw4cIBbbrklKYJLvKVUgDHGmGQVTb1LqrFKfmOMMXFhAcYYY0xcWIAxxhgTFykVYGw0ZWOMSR4pFWBsNGVjTLxt2bKFyy+/nK5du9KvXz8GDhzIM888E5dzBQbL3LhxI7169YrLOeIppQKMMcZUUVYGw4fD5s31PpSqcvHFF3PmmWdSUlLCihUrePTRRyktLY1BQlOPBRhjTGqbNg2WLoWpU+t9qAULFpCRkcG4ceMqlp1wwglMmDCBffv2UVBQwCmnnEKfPn1YuHAhACNGjOC9994DoE+fPtx2222AG3L//vvvZ8+ePZxzzjkVw/X/4x//qDYNGzduZPDgwfTt27fGCc8SzfrBGGNSU1aWmwgGYOVK6NMHZs6EzEzYu7dOh1yzZg19+/YNu+6ee+4B4P333+fDDz/kvPPOY926dZx55pm88cYb5OTk0Lx5c958803AjbD8k5/8hMzMTJ555hm+9a1vsX37ds444wwuvPBCRCTsedq1a8err75KZmYm69evJy8vj+XLl9fp9cSbBRhjTGoqKYGbb4YlS6B3b8jNhYEDIYazPv785z9n8eLFZGRk0KlTJyZMmABAbm4uJ5xwAuvWrWPw4MH8+c9/pkuXLlxwwQW8+uqrfPPNN2zcuJGTTz6Z8vJybr75Zl5//XWaNWvGF198wZYtWyL29C8vL2f8+PEVc7ysW7cuZq8n1izAGGNSU8eOkJ0NW7ZA+/ZuSOXsbKjHEC09e/bkqaeeqnh+zz33sH37dvr3789xxx0Xdp/TTjuN5cuX07VrV77//e+zfft2/v73v9OvXz8AHnnkEbZt28aKFStIT08nJyeHfYGcVxgzZsygffv2vPvuuxw+fJjMzMw6v554szoYY0zq2rkT8vNh/nx3v2NHvQ539tlns2/fPv76179WLPvmm28AOPPMMytmjFy3bh2fffYZJ598MhkZGRx//PE8/vjjnHHGGQwePJi77rqrYpbKXbt20a5dO9LT01m4cGGNQ+rv2rWLjh070qxZMx566KEGnwa5NlIqB2MzWhpjghQVVT6ePr3ehxMRnn32WSZNmsS0adNo27YtLVu2ZOrUqVx00UWMGzeOU045hebNmzNnzpyKicQGDx7Ma6+9xhFHHMHgwYMpLS2tCDBXXHEFo0aNon///hVTMVfnZz/7GT/84Q954oknGDp0aMRJzpJBSg3XH9C/f39N1kovY0zdNebh+pNVPIfrtyIyY4wxcWEBxhhjTFxYgDHGNCqpWKyfKPG+lhZgjDGNRmZmJjt27LAgEwOqyo4dO+LazDmlWpEZY1Jbp06dKC0tZdu2bYlOSkrIzMykU6dOcTu+BRhjTKORnp5Oly5dEp0MEyUrIjPGGBMXFmCMMcbERdIHGBHpKiKzReTJRKfFGGNM9OIaYETkARHZKiKrQ5YPF5GPRORjEbmxumOoaomqXh3PdBpjjIm9eFfyzwHuBioGBBKRNOAe4PtAKbBMRJ4D0oA7Qva/SlW3xjmNxhhj4iCuAUZVXxeRnJDFA4CPVbUEQEQeBS5S1TuAkXU9l4iMBcYCdO7cua6HMcYYEyOJqIM5Dvjc97zUWxaWiBwjIvcCfUTkpkjbqep9qtpfVfu3bds2dqk1xhhTJ4noBxNuHtCI3XJVdQcwLtJ6Y4wxySkROZhS4Hjf807AplgcWERGich9u3btisXhjDHG1EMiAswy4EQR6SIiGcCPgedicWBVfV5Vx2ZnZ8ficMYYY+oh3s2Ui4G3gZNFpFRErlbVg8B44BVgLfC4qq6J0fksB2OMMUnCZrQ0xhgTxGa0NMYYk9RSKsBYEZkxxiSPlAowVslvjDHJI6UCjDHGmOSRUgHGisiMMSZ5pFSAsSIyY4xJHikVYIwxxiSPlAowtSoiKyuD4cNh8+bIy8JtY4wxJiopFWAqisi2bq0+cABMmwZLl8LUqZGXhdsmmsBkjDEGVDXlbv3S0lQnTtQKEyeqHnWUu8/MVAV3W7my8nGkZf7nmZnBxwp3fGOMaeSA5RqD3+LUHCqmRw9dvnZt8MKVK6FPH/c4LQ1OPBHWroXu3eHUU93yd991y048sXK/9evdNuvXw6FDVY8V7viZmVBSAgUFMGcOdOgQ89dojDHxEquhYhIxH0z8HXEE5OYGB47evd2ygQNdgHnqKWjfHsrLoWNHl0eZP79y2QknwKefVj6/6ip3v2RJ5bEiHX/KFFesFihemzEjQRfCGGMSJ6XqYAKV/LpqFWzZ4gJHhw7ucfv27j47G/bvh/x8F1Dy82HHDti5M3jZV18FP9+3z+3rP1a44xcVueUzZ8KCBe5eBLKyEnZdjDEmEVKziKxnT10+bBhs3+4WtGkDY8bA3LluWVFR3Q48enTVY4Ue/7PPoHVrl9MJFMEFcjWqVmxmjEl6sSoiS80Ak+jh+idOdEEsPd0Vq40Z44rJJk1yQSjw3BhjkpAN1x9jxcXQq5ernunVCyZMCH5eXFx1m0jLqhS3zZrlisms2MwY05TEoilast369etX0dxu3jzVnj1VmzVz9/PmVV02frxqly6qCxaoHjigWlio2ry5uz9wwC1v29bdAttEWtalizue//hP37NJNT9fNTfXJSo3V7WgQLWsTHXTJtVhw9xjY4xJAsSomXLCg0Esb8Ao4D7oFzZwRAoK6ekumAT07Kn6xz+6+4CcHHfzC7essNAdLzTofDj8F66vTLt2wX1mrA+NMSbJxCrApGQdTL9+/fWuu5YzbBjccAPcfnvlui5d3P2GDZXL0tJcC+M1ayqf797t6uoDXV/S0lwd/eHDwfuFLuvZ09Xt+5ctXAg7Ro7m66w2zNw5holHz+XK//6ZZofD9KvJzIS9e2NzIYwxpg6sDqYaIjB0qAsOzzwTvO6zz1z3Fr/cXBcUArp3h3vvdfcBnTu7rjF+4ZZ9+GHwfgClpfCjb4ro/MR0lu4/lc5PTGfgcZ9Tcma+O3mgD01BQXDkM8aYRiwlA0xAaOCA8EHhkkugeXOX0ygvh4svhl//2t2Xl7vlX3/tboFtIi1LS3PH8/vtbyEnxwW99HR3f+fcjry5OrtqHx1rvmyMSREpHWBCA0ekoDBvHlx7rWs5lpkJzz4L48a5+8xMt/xPf3K3wDaRll17rTue//gbN8JttwWnbdAg0J07mdssnz7b5jO3WT4blu9wK23wTGNMKohFRU6y3fr27RexRVe4VmTz5tW9Miyc0ON36uQq+/0iNQaYN0+t4t8Yk1BYJX9kIv21Z8/lFBZCXl6iU+P6xhQWwuzZLueyeDFhGyAcapFF2oF97olV/BtjEsQGu6xGv36QyI78oQJBbsKEytFjDh6EW28N3u7w+hKKTriZszKWcGK/3qzPcINndn10SsMn2hhj6iml6mBqNaNlA8vLg9WrXcu21auhRw+Xk/H7/X0d2d0sm+NbbOFgm/Yc32ILi1ZmU7zQKv6NMY1PSgUYDcxomZ2d6KTUqLAQrr46uDHAtGkwuOdO0q/JR+bPJ/2afEacvoPJk7GKf2NMo5OSRWSNQaRis14riiDd22j6dI4ph7WZBE/fbANlGmMagZSs5E/4aMp11KuXGxdz6NDKZf6K/76yknfUKv6NMfFlPflTULhis1NblVDcIp+vO+fyn/29+bpzLk+0KuCZ6dbj3xiT3KyILImEKzbbdURHBg3JpuXCLdCpPenl5ZwxLJvz7+nAJdclNr3GGFMdy8EkmdDWZps2wXFZwfPLHJe5o8oQOMYYk2wswCS57t3h31cXwfTpcOqpMH06/766iGOPhSEnl/GKDOeskze7ic6MMSaJWIBJcuHqZS67DPbvh4dOmcZ5Ry3l4VOmUliIBRljTFKxOpgkF65e5rNtWWSyD54CVq7k+D59KGEm+67IhDxrWWaMSQ6Wg2kEQutlvi0lHB6dHzSXzMeDC8jRDaSluebOlpsxxiRao8jBiMjFwAVAO+AeVZ2f2BQl1lE9OvLFnmyO9+aS2f91OQs+zyYrpwP71rkhaK6+2m2bDIN9GmOaprjnYETkARHZKiKrQ5YPF5GPRORjEbmxumOo6rOq+lMgH7gsjsltFAoLYdn8nXx+Tj7lL87nobR8jjq0gylTKic0mz0bN8SMMcYkSEMUkc0BhvsXiEgacA9wPtADyBORHiJyioi8EHJr59v1N95+TVpeHuy/r4jz104n8/RT+elX09l3XxF5QyrHKxs0CD74wBWXWbGZMSYR4l5Epqqvi0hOyOIBwMeqWgIgIo8CF6nqHcDI0GOIiAB3Ai+r6jvhziMiY4GxAJ07d47dC0hSeXmVxV+9ekGnTgSNV/b7ljNo3twNPROYg8aKzYwxDalBxiLzAswLqtrLe34pMFxVr/GeXwmcrqrjI+x/PTAGWAasUtV7qztfYx2LrK4OZmTRvNyNV1a+dCXpA9x4ZeXNM0kvr2xVtnAhjB4N2dmVLdKSZVI2Y0zyaOxjkUmYZREjnar+WVX7qeq46oJLMs8HE0/NPy2h5Mx8SjJyyTyjNyUZuTxIAZQEj1dWWupus2bBvn3uvrDQNYG2ojRjTKwlKsCUAsf7nncCNtX3oI1pPpiY6tiRrn2y6dpyC4fatKdryy2kHZ3N4o87BM0j89vfQk6OawQQaAxw+eXwt79VDToWZIwx9ZWoALMMOFFEuohIBvBj4Ln6HrSp5mAA2Bk8XtngHju4+mr4fMI0dOlSPh8/lY0b4bbbgnd75hk3D40/6MyeDTfcYLkaY0z9xL0ORkSKgSFAG2ALcKuqzhaREcBMIA14QFVj1qi2qdXBhJWV5bIkBM8jcygjk7T9lfUyaWmuv+aaNZW7PvSQq6tZsCC4gcDkyVZfY0xT0GjqYFQ1T1U7qmq6qnZS1dne8pdU9SRV/XYsg4vxlJS4HE1uLu8c7g25uZScVcDA9ht488kyDg8bzptPbSYtDS65JHjXcEVp1q/GGFNbKTVUTJMuIgvVsaNrLub19mfLFrr2yWbS1A58cu00vpy/lE/GTuXaa2HePIKCTriitEGDsCkCjDG1klIBpslW8kcSUi/DrFnkXS6M3jmTo1cuYPTOmcy6W1hXmhUUdDp18vrV+CxeDMcea/UyxpjoNUg/mIZmdTARlJXBzTfDkiWVHWHWr3ejaAKsXAl9XF3NPslkxeMlDPx7AW+PncMl13UA4LHHrF7GmFQXqzqYRjHYZbREZBQwqlu3bolOSnIKLTYrL4errnL3S5ZUjMzMwIG8fNoUdl87le47l/LJ8qm0OGIGRUWuPgYq62UmTLAAY4wJz3IwTc3o0dCmDYwZA3PnwvbtcPTRUFTkavTLy+Grr8LmajQzE9lb2QKtvBwyMys3NcakBsvBmLopKqp8PH26ux892tXRBILOZ59B69ZBuZpnNg+k3f1T+J7vUP56GRt6xhgTKqUCjBWR1VG4oDNxYlBRWveB2Yz4VQce0rKKepnLvHoZG1DTGBOOtSIz4YW0QMtts4PJkwlqbdaihav0t/4yxphwrA7GRMc3MkBN9TItWkCPHlZsZkxj1Wh68psU4RsZoKJe5sgC3no4eMTm3/+einlobPBMY5o2CzAmOmFGBug+MJsrf9WBhQtdzmXhQjfn2Q03WLGZMSbFAowNFRNnEeplJkxwzZUnTHAjM996a/BuNsyMMU2T1cGYmOrVyxWLBTpkgs2kaUxjY3UwJnn4JjUrLHRNlf3FZpddBvv3W72MMU1NVAFGREaKiAUjE960abB0KUydSl4eVYrNwjVnvvxyKCiwgTONSWXRdrT8MfAnEXkKeFBVrUTdhG+6PHMmeZmZ5O0NntRs0KDK3YqL3RQBBw64nI110DQmNUWVK1HVnwB9gE+AB0XkbREZKyKt45q6WrJK/gYWpukyBQWwIbjpcvfuLogETJ4M48e7vjLW0syY1BV1sZeqfgU8BTwKdAQuAd4RkQlxSlutWU/+Bham6TLZ2dChQ9BmofUya9dW1sMEWEszY1JPVEVkInIhUAB8G3gIGKCqW0XkCGAtMCt+STRJLdB02T86c4hAsdeECS6IpKfDFVcEF4ctXuxyOsaY1BFVM2URKQLuV9XXw6w7R1Vfi0fi6sqaKSeBsjJXXDZnTpUcTXGxy73Mnm2TlxmTjBq6mXJZaHARkakAyRZcTJLwtSwLFa6lmQUXY1JPtDmYd1S1b8iy91T1O3FLWT1YDiaBIgyKSWYm+FqWGWOSV4PkYETkOhF5H8gVkfd8tw3Ae/U9uUlBUbYsC6e42PWJsb4xxqSGmir55wEvA3cAN/qW71bVnXFLlWm8QluWlZdXtiyrQ70MWNGZMY1VTXUwqqobgZ8Du303ROTo+Cat9qwfTJIIGRSTHTvc8mrqZSZPdsHFRmE2JnVUWwcjIi+o6kivSEwB8a1WVe0a7wTWhdXBJJko6mXS0twm6emVu5WXu00OHWrg9BrTxDVIHYyqjvTuu6hqV+8+cEvK4GKSUKR6mSVLKgbJDO3tD9Y3xpjGrto6GBHpW916VX0ntskxKSlSvcycORVFZoWFM7j66vB9Y4wxjVNNRWQLq9lXVfXs2Cep/qyILAmNHg1t2rge//36VZZ7+YrMDqZn0vukvTZnjDEJFqsismpzMKo6tLr1xkStqKjy8eefw803uyKyQJHZwIE0nzKF1cGNyygudrkYCzrGND41FZGdraoLROQH4dar6tPxSZZJaVE2ZS5e2MGaLhvTiNXUTPks735UmNvIOKbLpLoomjJb02VjGreohoppbKwOphGK0JRZMzORkpKKXE35MR2s6bIxcdagg12KyDEi8mcReUdEVojIn0TkmPqe3JgKYZoyP3NkAW89vCEoV2NNl41pPKId7PJV4HXgYW/RFcAQVT03jmmrM8vBNFITJ7rGAOnpUF7O4V1f0eywy6qUL11J+oDK1mbND9jAmcbES4O0IvM5WlX/4Ht+u4hcXN+TR0NEugO/ANoAr6nqXxvivCYBQiYva/bZZ5Rsaw1LlnDiGb1Zn+Fam3V9dEqiU2qMiUK088EsFJEfi0gz7/Y/wIs17SQiD4jIVhFZHbJ8uIh8JCIfi8iNkfYHUNW1qjoO+B+g3hHVJLGiIpg+HU491d0/+SRd+2TTteUWDrVpT9eWW+jaJ5vihR1s1GVjGoGahuvfLSJfAdfiRlY+4N0eBSZFcfw5wPCQY6YB9wDnAz2APBHpISKniMgLIbd23j4XAosBm9ysqQlpbbZh+Q4KC+Fvvyuj/Nzh/O33mykstCBjTDKKeysyEckBXlDVXt7zgcDvVHWY9/wmAFW9I4pjvaiqF0RYNxYYC9C5c+d+n376aWxegEkqvXrBrFkw9LlJMHcujBnDwgtnMGECrF5d8/7GmJo1dB0MInIUcCKQGVgWOo1ylI4DPvc9LwVOr+a8Q4AfAC2AlyJtp6r3AfeBq+SvQ7pMI7BsTRZZZwc3Zx46cybLyKRXLxtmxphkElWAEZFrcBXtnYBVwBnA20BdxiKTMMsiBgRVXQQsiurAIqOAUd26datDskxjcP5JJRR3uZmOn1YOM7MiYyAXfzCFolnW49+YZBJtJf8vgNOAT73xyfoA2+p4zlLgeN/zTsCmOh4riKo+r6pjs7OzY3E4k4Su/V1H/vl2NuVfbEHbt6f8iy0sXp3NmF93sB7/xiSZaAPMPlXdByAiLVT1Q+DkOp5zGXCiiHQRkQzgx8BzdTxWEJvRMvXl5cGQ7+xkXkY+fbfNZ15GPkcd3sGtt+LGMfPmlxk0yA2QaYxJnGgDTKmIHAk8C7wqIv8gilyHiBTjitJOFpFSEblaVQ8C44FXgLXA46q6pi6JD2U5mKahyxtFjNk+nZWHT2XM9ulM61nkJisL6fF/7LFYc2ZjEqjWrchE5CwgG/inqh6IS6rqyXryNy0HM7JoXu4q/v09/veSyZIFe6tMYGb1MsZUr0HHIvNO2FdErge+A5QmY3CxIrKmqfmnJZScmU9JRi6ZZ/SmJCOXx44o4J0nN1i9jDEJFO1gl78F5gLH4IZseVBEfhPPhNWFFZE1UR07VunxX/ZNNgMu7GD1MsYkULQ5mDzgNFW9VVVvxTVTviJ+yTKmlkJ6/HfN3hG2XsZGYjam4UQ7mvLLQJ6qfuk9PxJ4WFWTatIxXz+Yn65fvz7RyTEJFKlexkZiNqZmDVIHIyKzROTPwH5gjYjMEZEHgdXAnvqePNasiMwEhKuXKTmrgOafbUh00oxpMmrqyR9oirUCeMa3fFFcUmNMrHj1Mry/hUPp7aG8HPpkQ4cOiU6ZMU1GtQFGVecGHnudIk/ynn6kquXxTJgx9RYyvwzbtyc6RcY0KdG2IhsCrMcNs/8XYJ2InBm/ZNWNNVM2QULnlykqorg4QZ0vfa3ZjGkqom1F9kfgPFU9S1XPBIYBM+KXrLqxOhhTneJiN8ryrFmwb5+7DzeXTFyCkK81mzFNRbQBJl1VPwo8UdV1QHp8kmRMfEye7DpbDs0tI33UcIZ238zs2XDDDZUB5fjj4Re/iBCEosiFhAangxlZIAIzZ9L3ywUwc6Z7npXVYK/bmESJNsCsEJHZIjLEu/0dV/FvTKOxdq0bzt+fmygthdLSyoDSvLkLDps3hxkBoIZcSLgcUu/WJRS3yOfrzrn8Z39vvu6cyxOtCnhmurVmM02AqtZ4w0329f+Ap3GtySYBLaLZtyFvwCjgvm7duqkxofZKpiq428qVFY/3Sqbqpk2qw4ZpRynT+fNVe/as3O9wZvj9NDMz6Pg9e6ouWBB8zpwc1dmtf6F61FGq7dqpHnWUfvbDie743jm1rCzOr9yY2gGWawx+k2vMwYhIM2CFqk5X1R+o6iWqOkNV98cr6NWVWh2MqcbLd5fwZCuXmyjv6XITD1DA83/aUJE7ufMolzvxDynz1kMlPH6E61OT1q+yTw0bgnMhFTkkn88+g+a7g0cZOC5zhzu+1cuYFBdtT/5HgJtU9bP4J6n+bDRlE8lH50+kw6tF7DuUTmZaOS0PfUVzDrmV3hTMAPskk7T9e1m8GC67DP6wZyJXNS+ieVY6B/eW85CMIeveGUEjM/fq5YrFhg6tXNali7v3x6JDLbJIOxA87TMAmZmw10YZMInX0KMpd8T15H9NRJ4L3Op7cmMa2sltd5J9fT7tV80n+/p8Ng24uEquZl5mARe2XcKCjOH8btxmWrSAiwbtJP2afGT+fNKvyWfE6TuqjMxcWOimBFi40PXrXLgQvv7a3fzLvtvejTJAbm7FtM8UVM0RGdPoRVOOBpwV7haLMrp43Pr161e3gkfTJH04/Bf6ZdpRupl2+mXaUbr2/ImqEye6epOJE7VZM9UDB4L3OXBAtVmzqseaN0/1rJM26T8ZpmeeVKbz5rllPXu67Xv2dM/1F78IqpfRiRPdAULqZcLua0ycEaM6mGp78otIJjAO6Aa8D8xWNyOlMSnj5LY74fp8sseMgX79yH55JrxMRfHVIWZyqFUm7K8svoo0MnNeHuQtnQbbljJsxFTIm1GxPMjLVUcZKC6G8vHTGLlzKS/0nMqyy2fw4ovw0LQyBv69gLfHzuHKX3UIfzxjklF10Qd4DHgYuBY3XfKfYhHV4nXDWpGZ+tq0STU/XzU31z3PzdVPzirQ044v0wULXM5lwQLVLl1Ux48Pzl2Up0fX2iwc/74Hllbue6B5ZlBuasGC4BZuxsQDMcrB1PSD/b7vcXPgnVicNN43KyIz9RKm+Cq0qGr8eBdk/EFnwPGb9JMzg4OTFhSEbYYcery+HTfppmHB+5aTFjZgfUPNAcuY+ohVgKmpkr9iQEu1ojHTVIRMXsaOHeTlwerVcOiQu1+40BsVwDcl851zO/Lm6mzYsgXat3f32VVHcA7XIfOdso58qcH7PnPkVTxIflBjgLLhBQw/yRoDmMahpuH6TxWRr7zHAmR5zwVQVf1WXFNnTCIUFVU+nj497CYVfV7KylwLsDlzGDSoA4/u3AmT8qsdwbliyBqvOfPQoZCTA2vf3En3sZX7dvvXdt746mjKv9hC8/btObi3nJc3ZTPuXptywDQOUfWDaWysH4yJt4o+L89NcoFkzBgWXjiDCRNcDqc6aWku55LuG83voYdg9GhYsMAFrsWLXZPnRzNGs3Z7G2buHMPEo+dyZvftdHmjKPLBjYmBWPWDSXh9STxuVgdj4i1SpXx5es31I+GGlFmwQLVTp7o1SW6Ipszhml+b1EVDVPI31psFGBN3m1yF/icZudqsmeonGa61WTTjis2bV7WBQJcudQsonTqptm2ruviJTXrovGG6+MmyWh0r2vN16aL62Q8n6mFvLLVYn8MkFwswFmBMokXqLBmFuuY6QoNTTo5qhw7qOofGqSlztYOEmpQUqwCTUnUwIjIKGNWtW7efrl+/PtHJMalu9Gho0ya4Qr8ovvUjoeOd7ZUssqg6rtleMsnS2IxrdlyzMj6/8maaLV3iWjd0787h0wfSqWgKmw5bg4NUFKs6mJpakTUqqvo88Hz//v1/mui0mCYgitZmsRY6YvP5J5XwSM7NfDV/Cd0DTZlzBvLjkin8O0bnPKpHR77Yk83xgSbU5eV8sSebo3tYcDHVi3awS2NMEuje3bUwC7j2dx15aXE2HZttQdu3p/yLLbz8VjbfOa9DzKZ9LiyEZfN38vk5+ZS/OJ/Pz8ln2Ss7KCys/+sxqc0CjDGNSOiIzR06QPahnTxxRD59t81nXkY+vTru4MUXI0z7XAd5ebD/viLOXzudzNNP5fy109l/X5GNh2ZqlFJ1MAHWD8aksuJi11nTqw6hsJAa56VZuJCo+ugYA7Grg7EAY0xj5BtBIHQomnAdOcvL3Xxmhw41bDJN49TQE44ZY5JJNdMth9bTgHt+7LHErF7GmGikVCsyY1JeVpbLnkBls+SZM4OmWw7U08yeXTnszGWXuV1mzQoeigZsbhkTP5aDMaYxKSlxIzxXM91yXp6ro5kwwcWdCROgRQt47LHg0Z9nz6bKtM/GxJIFGGMak44d3RQANUwJEDq9wKZNwf1nwD3/4AMrNjPx0ygCjIi0FJEVIjIy0WkxJuHCzFdTk3D1Mr//PTRvHrvmzMaEimsrMhF5ABgJbFXVXr7lw4E/AWnA/ap6Zw3HuQ34Glijqi/UdF5rRWZMsMAkZ/56mWHD4IYb4PbbK7ez5swGGk8rsjnAcP8CEUkD7gHOB3oAeSLSQ0ROEZEXQm7tRORc4ANgS5zTakzKCtTL/H5cGQsyhvO7cZs5eBBuvTV4u0GDXP+aaBQXR1e8Fu12JvXEtRWZqr4uIjkhiwcAH6tqCYCIPApcpKp34HI7QURkKNASF4z2ishLqno4zHZjgbEAnTt3junrMCYV5OVB3tJpsG0pw0ZMpderM1i8OLhD5uLFrjitJuFyROFapUW7nUlRsRiSubobkAOs9j2/FFcsFnh+JXB3FMfJB0ZGc04brt8YVd20SXXYMDdHTWbkIffrMi9NpEnTQqcJiHY7k1yI0XD9iajklzDLaqwIUtU5WkP9i4iMEpH7du3aVefEGZMy/J0xIzRvfvmeDUHFZpMnezmLsjIYPhw2bw576NBRnSF88Vq025nUlIgAUwoc73veCdgUiwOr6vOqOjY7OzsWhzOmccrKAhHXAXPBAnd/7LHw0ENVmjdfcl0HFo2YxrCjlvLvEVMri62qGSkAIo8WEFq8Fu12JkXFIhtU3Y2qRWTNgRKgC5ABvAv0jOU5rYjMNGmbNqnm56vm5rrnubmqBQWql16qOmmS6qpV7j4tLWyxWeDWR3zLMoNnr4x22uf6Tg9tEoPGMGUyUAyUAeW4nMvV3vIRwDrgE6AwhucbBdzXrVu3GF1mYxqpaKZzDhOINp5xmT53xGW6p3OuHjiguqdzrj7eqkCf/ktZld2jnfa5rtNDm8SJVYCx0ZSNSUXRTuc8caJbnp4O5eUUyRiGDlGOX1i57POzx3D+hzOsb0wTYlMmhyEio4BR3bp1S3RSjEmsaKdzDowK4AUimbGd47IIWnbc1u31r5SvZnoBk7osB2OMqRC3ycomTXI5qTFjYMaMeqfTxFdj6clvjGlEQqdkXrjQPS8sjLBDDc2Zw7ZoE3HLTcpLqQBj/WCMqZ9wQ/1X9I0JJ6Q5c+iwMM/8sebpBUzqsiIyY0zthZv4DNgnmbz92t6gYWFePnkiJ/+nstGAFZMlPysiM8YkTpiRAZ45soAVT2yoMqnZR2/XfnoBkxpSqhWZMaaBhE58Vl7Op19mM/LC4BZigwZB5u4iDgUaslXXos2knJTKwVgdjDENKGTis67ZO2xYGBPE6mCMMTERaWj+ahsJmKRkHS2NMUklEEQmTHCjJXfvbsGlqbMAY4yJmbw8CyimktXBGGNiq6bOl6bJSKkAozYfjDGJV8NcMqbpSKkAY4xJoPoOC1PHnE/o6AHFxbVOuYkTCzDGmNiIMC1zNMPCFBdDUa9p7HxlKUU9p0YdJAIt12bNcgMLzJrlnluQSQ5WyW+MiY0wnS/Jzq5xeP6DGVnklbthZ8qXrmT0gD5w+UwOjsmk+YG91e47ebJrFh0Y/TkwesCECdbYIBlYDsYYEzs7az8szLldSigblg+5uaSf1htycykbXsCANhtqLPpau9b1ufEbNIj6z19jYiKlcjA24ZgxCRbtRGc+b3zckXbDsmFpZc7nyxOyWVnWgQWPBHfahOCcSffubp1//hobPSB5pFQOxlqRGdP4dO8O2z4MzvmsXbyDnByqDJw5eXLwvrWev8Y0qJQKMMaY5Bfa6mvoUPjux0UsHDWd8h6nsnDUdH74dRG33UZQy7JwRV+1nr/GNKiUKiIzxiS3SOOVXXBB8BAznTq5m79PzeILZ4Qt+rLRA5KX5WCMMQ2motVXbhnpo4YztPtmZs92RVurV8OhQ+5+w5Yshp7t+tSUv+L61Aw9W1i1LrZTLVsfmviyAGOMaTAVrb58OZNwRV/NPy2h5Mx8SjJyyTyjNyUZuZScVUDzz+o51bKvyM360MSfFZEZYxrM15pFekbwVMvpM2fytWQCvj4vHTvStU82vL+FQ+len5o+NfepqZEvsE1+dYb1oYkzCzDGmAbz8t0lHPr1zZx/9BIyevbmQOdcXto5kObTpnBJ6MaBPjVjxsDcubB9e91PnJXlsilQEdhWMxMdkQl7KwOb9aGJrZQKMNYPxpjkdsnPOvLR89kcfHULezLak5lWzinnZZN7XZicSR361ERUUgI33wxLllQMY/PM5oG0u38K3/NtZn1oYiul6mCsH4wxye/ktjvJvj6f9qvmk319Prltau7tX6OaBsoMHcZmyxa6D8zmyl91sD40cZRSORhjTCMQy5xJgH+KgBkzwm8TUuSWu317RR8am4EzPkRVE52GmOvfv78uX7480ckwxsRbmLoVwPW63Fv9QJkmMhFZoar963uclCoiM8Y0MfWYIqCCzcAZNxZgjDGNV5i6lWimCAhiM3DGjdXBGGMat7o2Zw5XvDZzphWvxZAFGGNM41bXRgNhmi4zcCBMmRLzJDZVVkRmjElNdWi6XOviNVMtCzDGmNQUTd1KHWbgNNGzZsrGmNRiTZfrrck0UxaRISLyhojcKyJDEp0eY0ySi9B0+ZnpG2xo/gYW1wAjIg+IyFYRWR2yfLiIfCQiH4vIjTUcRoE9QCZQGq+0GmNSRJi6lQ83Z/PL/+tgQ/M3sLgWkYnImbjgUKSqvbxlacA64Pu4gLEMyAPSgDtCDnEVsF1VD4tIe2C6ql5R03mtiMyYpm3D4NG8vrYNM3eOYeLRc2m5dzvHvFBUMTQ/uLHHJkxwE5yZYLEqIotrM2VVfV1EckIWDwA+VtUSABF5FLhIVe8ARlZzuP8CLSKtFJGxwFiAzp071yfZxphGrLgYCr8oYvYTsHQQLF48nbPPhqKQ8g8bmj/+ElEHcxzwue95qbcsLBH5gYj8DXgIuDvSdqp6n6r2V9X+bdu2jVlijTGNS8W0zEMhPd3d5+TAb38bvJ0NzR9/iQgwEmZZxHI6VX1aVa9V1ctUdVG1BxYZJSL37dq1q75pNMY0UhXTMvvcdhts3IgNzd/AEhFgSoHjfc87AZticWCbD8YY0727y534derkbhMmuNbKEybY0PwNIREBZhlwooh0EZEM4MfAc7E4sOVgjDGFhS53EppbmTbNVegfOuTuLbjEX1wr+UWkGBgCtBGRUuBWVZ0tIuOBV3Atxx5Q1TWxOJ+qPg88379//5/G4njGmMYnEDhsIrHEs578xhhjgjSZnvy1YUVkxhiTPFIqwFglvzHGJI+UCjDGGBNrxcXYGGZ1lFITjonIKGBUt27dEp0UY0wKKC52rdJmz3Z9axYvdi3SwBoNRCOlcjBWRGaMqS9/jqWgAC6/PHhUgNmzXas0U7OUCjDGGFMfgRxLYNTl8nJ45JHgYjEbwyx6KRVgrBWZMaY+Qscx6969std/QMUYZjVNyWxSK8BYEZkxpj5CxzErLIS774YPPggzhlk0UzI3cSlVyW+MMfURGMcsMG9MXh6sWQN33eXGMOveHdaVZtH88pApmWfOtCmZw0ipHIwxxtRHuHHM5s2DBx+sHMOs+afhp2Rmw4ZEJz/pWA7GGGM8UY1jFjolc3m5e96hQ0LSnMxSKsBYPxhjTH3l5UXRx2XnTpeLGTMG5s6F7dsbImmNjg12aYwxJogNdmmMMQ3FmiTXiQUYY4ypSUiTZBufLDopVQdjjDExUVbmWob9+9+uSz8ENUm+RDLp8NpeG5+sBimVg7Ge/MaYmAjkWK68skqT5GeOLGDFExtsfLIopFQOxqZMNsbUS1ZW1RxLgNck+dMvsxl5YXCTZBufLLyUysEYY0y9lITpRPntb8PYsTB/PuTn0zV7B4sXB+9WMT6ZCZJSORhjjKmXcJ0ox4yBGTPc+unT+fo0mHh11TlirIisKgswxpjGJVABP2dOfHrP19CJMqre/gawjpbGmMZm0iT3w+/PWZiYilVHS8vBGGMah3AV8DaKcVKzSn5jTOMQrgLeRjFOaimVg7HBLo1JYTaKcaOTUjkYm9HSmBQXqID3mgyzY0eiU2SqkVI5GGNMiisqqnw8fXri0mGiklI5GGOMMcnDAowxxpi4sABjjDEmLizAGGOMiQsLMMYYY+LCAowxxpi4sABjjDEmLizAGGOMiYuk72gpIs2APwDfApar6twEJ8kYY0wU4pqDEZEHRGSriKwOWT5cRD4SkY9F5MYaDnMRcBxQDpTGK63GGGNiK945mDnA3UDF+A4ikgbcA3wfFzCWichzQBpwR8j+VwEnA2+r6t9E5EngtTin2RhjTAzENcCo6usikhOyeADwsaqWAIjIo8BFqnoHMDL0GCJSChzwnh6KdC4RGQuM9Z7uD8011VM2sCvG20faJtzyaJb5n/sftwGCp+Srv1hfj9pci3DLa/M81tejoa5FpHXJ9NlI9Pck3HL7bFS/LNJn4+SaEhsVVY3rDcgBVvueXwrc73t+JXB3NfsfAcwGZgE/j/Kcy2P8Gu6L9faRtgm3PJpl/uchj2N6LeJxPWpzLWp67VFcm6T+bFS3Ptk/G4n+nthnI/k+G4mo5JcwyyLO26yq3wBXxy85UXk+DttH2ibc8miWPV/NuliL9fWozbUIt7y2z2Opoa5FpHXJ9NlI9Pck3HL7bFS/LK6fDfGiVdx4RWQvqGov7/lA4HeqOsx7fhOAuiKyWJ1zucZgPulUYNcimF2PSnYtgtn1qBSra5GIfjDLgBNFpIuIZAA/Bp6L8Tnui/HxGjO7FsHselSyaxHMrkelmFyLuOZgRKQYGIKrPNsC3Kqqs0VkBDAT13LsAVWdHLdEGGOMSYi4F5EZY4xpmmyoGGOMMXFhAcYYY0xcWIAxxhgTFykfYESkpYjMFZG/i8gViU5PoolIVxGZ7Q270+SJyMXeZ+MfInJeotOTSCLSXUTuFZEnReS6RKcn0bzfjhUiUmWEkaZGRIaIyBve52NItPs1ygBTy0E0fwA8qao/BS5s8MQ2gNpcD1UtUdVEd1yNq1pej2e9z0Y+cFkCkhtXtbwWa1V1HPA/QMr1B6nD4Lu/Bh5v2FQ2nFpeDwX2AJnUZtDhWA6N0FA34EygL8FD0KQBnwBdgQzgXaAHcBPQ29tmXqLTnujr4Vv/ZKLTnWTX449A30SnPdHXAvcn7C3g8kSnPZHXAjgX10cvHxiZ6LQnwfVo5q1vDzwS7TkaZQ5GVV8HdoYsrhhEU1UPAI/ihvovBTp52zTK11uTWl6PlFeb6yHOVOBlVX2nodMab7X9bKjqc6r6XSDlipNreS2GAmcAlwM/9ealSim1uR6qethb/1+gRbTnSPoJx2rhOOBz3/NS4HTgz8DdInIB8R+jK5mEvR4icgwwGegjIjdpDIfoSXKRPh8TcP9Ws0Wkm6rem4jENbBIn40huCLlFsBLDZ+shAh7LVR1PICI5APbfT+wqS7SZ+MHwDDgSNwULFFJpQATdhBNVf0aKGjoxCSBSNdjBzCuoROTBCJdjz/j/oQ0JZGuxSJgUcMmJeGqHXxXVec0XFKSQqTPxtPA07U9WCpl+0qB433POwGbEpSWZGDXI5hdj0p2LSrZtQgW0+uRSgGmIQbRbEzsegSz61HJrkUluxbBYno9GmWA8QbRfBs4WURKReRqVT0IjAdeAdYCj6vqmkSms6HY9Qhm16OSXYtKdi2CNcT1sMEujTHGxEWjzMEYY4xJfhZgjDHGxIUFGGOMMXFhAcYYY0xcWIAxxhgTFxZgjDHGxIUFGBMTInJIRFaJyBoReVdE/l9ggEAR6S8iEYdjEZEcEbm84VJb5dyra94yaJ98ETk2XmmqRTr2JDoNobz5dX7rPZ4jIpfW8ThtReSfsU2daWgWYEys7FXV3qraE/g+MAK4FUBVl6vq9dXsm4MbtbaxyAcSHmCS1A3AX+p7EFXdBpSJyPfqnySTKBZgTMyp6lZgLDDeGw5/iIi8ACAiZ3k5nVUislJEWgN3AoO9ZZO8XMUbIvKOd/uut+8QEVkkbsbFD0XkERERb91pIvKWl3taKiKtRSRNRP5PRJaJyHsicm2EJDcXN+vpe96xj/CO2U9E/i1uVsNXRKSj94+8P/CIl96zRORpb/uLRGSviGSISKaIlHjLvy0i//SO84aI5HrL24rIU176lgV+TEXkd+Img1okIiUiEjE4i8hk7zUvEZH23rITROQ17/W8JiKdveVzROSvIrLQO+5Z3nnWisgc3zHPE5G3vWv/hIi08pbfKSIfeMe9K0xaTgL2q+r2MOv+4J2/mYhsFJEp3jmWi0hf7/p+IiL+gVifJQWnDWhSEj3pjd1S4wbsCbPsv7gJioYAL3jLnge+5z1uhRvRu2K9t/wIINN7fCKw3Hs8BNiFG4CvGW6Yi0G4iZFKgNO87b7lHXcs8BtvWQtgOdAlJI05uNFzA2l6APhfIB038VZbb/llwAPe40VAf+9xc2CD9/gu3FhO3wPOAoq95a8BJ3qPTwcWeI/nAYO8x52Btd7j33nnbgG0AXYA6WGurwKjvMfTfK/1eWCM9/gq4Fnv8Rzc/B6Cm/PkK+AU71quAHp753sdaOnt82vgt8DRwEdUjv5xZJj0FAB/9D2fA1zqpe1vvn03Atd5j2cA7wGtgbbAVt/+xwHvJ/qzbbe631JpuH6TfMIN/f0mMF1EHgGeVtVSLxPil46bw6c3cAg4ybduqaqWAojIKlyA2AWUqeoyAFX9ylt/HvAdXz1ANi5gbQg53+eq+qb3+GHgeuCfQC/gVS99aUBZaEJV9aC4qWW74yZrmo6bKTANeMP79/9d4Anf6wxM2HQu0MO3/Ftejg7gRVXdD+wXka24QB06Ve0B4AXv8Qpc0STAQNy8LgAP4X7gA55XVRWR94Etqvq+d63W4K5lJ9wMhm966crABfKvgH3A/SLyou+8fh2BbSHLbgH+o6pjQ5YHBlB8H2ilqruB3SKyT0SOVNUvga1YUWSjZgHGxIWIdMUFh61A98ByVb3T+4EaASwRkXPD7D4J2AKcivt3vc+3br/v8SHcZ1jwzeHhTwYwQVVfqSG5ofuqt+8aVR1Yw74AbwDnA+XAv3D/3NNwOaFmwJeq2jvMfs2Agaq6NyjR7oc93OsMVa6qWsM2EPz6Asc9HHKOw97+h4BXVTUv9CAiMgA4BzfC7njg7JBN9uKCuN8yoJ+IHK2q/tkTa0oHuPnfg66NaVysDsbEnIi0Be4F7vb9AAbWfVtV31fVqbgiq1xgN66IJCAblyM5DFyJ+7GuzofAsSJymneO1iLSHDci7HUiku4tP0lEWobZv7OIBAJJHrAYVxzUNrBcRNJFpKe3TWh6XwcmAm+rq5w+xntda7zc1AYR+ZF3HBGRU7395uN+qAPXpncNrzNab+GCALg6jMW12HcJ8D0R6eal6QjvurUCslX1JdxrDZfWtUC3kGX/xNWxvejLnUXrJKBWLfxMcrEcjImVLK/IKh04iCuamR5mu4kiMhT3T/kD4GXcv9aDIvIu7t//X4CnvB/lhcDX1Z1YVQ+IyGXALBHJwv3rPRe4H1fs8464bME24OIwh1gLjBGRvwHrgb96x7wU+LOIZOO+KzOBNV4a7xWRvbjiqP/girBe9473Hq4uIRBcrwD+KiK/8a7Po8C7uKK4e0TkPe/4rxOb2UavBx4QkV95rznqGV1VdZu4aYKLRSRQlPcbXFD9h4hk4nJ3k8Ls/jrwRxER/x8LVX3CCy7PiciIWryOocCLtdjeJBkbrt8YEzMi8idcPc+/YnCs14GLVPW/9U+ZSQQrIjPGxNIUXCvAevGKWadbcGncLAdjjDEmLiwHY4wxJi4swBhjjIkLCzDGGGPiwgKMMcaYuLAAY4wxJi7+P4hc5CYW2BTCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 435x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Library\n",
    "from sklearn.neighbors import KernelDensity\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Constants\n",
    "N_PROB_POINTS = 70\n",
    "PAPER_FIG_ASPECT_RATION = 2.9/2.4\n",
    "\n",
    "#Create and fit the kernel density\n",
    "kd = KernelDensity(bandwidth=2, kernel='gaussian')\n",
    "data = dist_km_bright.values.reshape((len(dist_km_bright.values), 1))\n",
    "kd.fit(data)\n",
    "\n",
    "#Compute probabilty for some values\n",
    "dist_space = np.logspace(0, 5, N_PROB_POINTS)\n",
    "dist_space = dist_space.reshape((len(dist_space), 1))\n",
    "density_bright = np.exp(kd.score_samples(dist_space))\n",
    "\n",
    "#Create and fit the kernel density\n",
    "kd_2 = KernelDensity(bandwidth=2, kernel='gaussian')\n",
    "data = dist_km_gowalla.values.reshape((len(dist_km_gowalla.values), 1))\n",
    "kd_2.fit(data)\n",
    "\n",
    "#Compute probabilty for some values\n",
    "density_gowalla = np.exp(kd_2.score_samples(dist_space))\n",
    "\n",
    "#Plot\n",
    "plt.figure(figsize=(5*PAPER_FIG_ASPECT_RATION,5))\n",
    "plt.title(\"Replication of Figure 2a) from Friendship and Mobility paper\")\n",
    "plt.plot(dist_space[:], density_bright, \"o\", fillstyle = \"none\", color = \"blue\", label = \"Brightkite\")\n",
    "plt.plot(dist_space[:], density_gowalla, \"*\", fillstyle = \"none\", color = \"red\", label = \"Gowalla\")\n",
    "plt.ylim((10**-6, 10**0))\n",
    "plt.xlim((10**0, 10**5))\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"Distance between homes (km)\")\n",
    "plt.ylabel(\"Probability\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this milestone, I replicated the figure 2a of *Friendship and Mobility*. To do so, I did the following:\n",
    "\n",
    "0. Apply some preprocessing to clean data and remove useless information;\n",
    "1. Convert long-lat check-ins' positions to x-y position (Mercator);\n",
    "2. Assign to each check-in its square in the map;\n",
    "3. Find for each user the square that has more check-ins;\n",
    "4. Mean all the check-ins in that square to estimate home position;\n",
    "5. Convert homes position from x-y to long-lat;\n",
    "6. Compute the distance (Haversine) between two friends homes;\n",
    "7. Estimate the probability distribution of the distance with a Gaussian Kernel Density Estimator;\n",
    "8. Plot the result.\n",
    "\n",
    "The plot obtained is quite similar to the one presented in the paper. Qualitatively the results are identical: the behaviour of the plot is the same, and we can notice the kink around 100km. Quantitatively the result is also really near to what the authors presented. However, there may be a little difference around the 1km values where the probability seems higher in their plot (around ~0.1 wherein my plot the probability is more about ~0.05). We cannot be sure that there are differences or not as it is hard to say due to the size of their plot, the fact that they don't display the 0.1 probabilty value and that there is no interpolation between the points. However, the difference may come from some difference in assumptions or some errors.\n",
    "\n",
    "To replicate their data, some assumption had to be done due to a lack of information about the method that they used:\n",
    "1. The preprocessing cleaning step may be different;\n",
    "2. The square position assignment: I projected the data to Mercator in order to find the home square and then projected back its position;\n",
    "3. During this projection the squares were created around a reference position that is the centre of long-lat system (Null Island);\n",
    "4. When there was a tie between two possible square home we kept the first one in our order, but they may have done differently;\n",
    "5. I computed the distance between two points with Haversine distance, they may have used an other metric,\n",
    "6. A Gaussian Kernel Density was used as probability distribution estimator, they did not mentionned how they computed the distribution.\n",
    "\n",
    "Some numerical error may have come from:\n",
    "1. Conversion from long-lat to x-y and reverse (some extremes values had to be drop out and there exist imprecisions);\n",
    "2. Computation of the distance between the points.\n",
    "\n",
    "In terms of computation efficiency, on my computer, the code takes about a little more than 1 minute to run (1 minute 10 seconds). The most consuming blocks of code are the data loading (due to the size of the data to load) and the probabilty estimator (due to the complexity of the task with gaussian kernel and the size of the dataset). This duration seems reasonable to me knowing the size of the dataset.\n",
    "\n",
    "To conclude, I was able to reproduce quite reliably their results, asserting the existence of the kink around 100km."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}