Recursive_Search_geolife_500tessellation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4ac95068",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')\n",
    "\n",
    "import attack\n",
    "import data_loader as dl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1ff41bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gp\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "import libpysal\n",
    "import itertools\n",
    "from joblib import Parallel, delayed\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely import Point\n",
    "import math\n",
    "from sklearn import metrics\n",
    "from tqdm import tqdm\n",
    "import os\n",
    "import datetime\n",
    "tqdm.pandas()\n",
    "\n",
    "import sys\n",
    "sys.stdout = open('logfile_recursive_geolife_500tessellation', 'w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e6af6820",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "gdf = gp.read_file('data/geolife/geolife_clustering_HL_500tessellation.geojson').to_crs(4326)\n",
    "\n",
    "# format clusering_HL columns as object\n",
    "gdf['clustering_HL'] = gdf['clustering_HL'].astype('object')\n",
    "gdf['START_POINT'] = gdf.geometry.apply(lambda x: Point(x.coords[0]))\n",
    "gdf['END_POINT'] = gdf.geometry.apply(lambda x: Point(x.coords[-1]))\n",
    "\n",
    "gdf = gdf.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "be91c235",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.8/dist-packages/skmob/tessellation/tilers.py:117: ShapelyDeprecationWarning: The 'cascaded_union()' function is deprecated. Use 'unary_union()' instead.\n",
      "  base_shape = gpd.GeoSeries(cascaded_union(polygons), crs=base_shape.crs)\n"
     ]
    }
   ],
   "source": [
    "from skmob.tessellation import tilers\n",
    "import skmob\n",
    "from shapely.geometry import Polygon\n",
    "\n",
    "\n",
    "MIN_LNG, MIN_LAT, MAX_LNG, MAX_LAT = 116.08, 39.66, 116.69, 40.27  # raw_points_gdf.geometry.total_bounds\n",
    "\n",
    "TILE_DIAMETER_IN_METERS = (\n",
    "    500  # approximately. For h3 the most appropriate resolution is found\n",
    ")\n",
    "TILE_TYPE = \"squared\"  # other option: \"squared\", \"h3_tessellation\"\n",
    "####################################################################\n",
    "\n",
    "base_shape = gp.GeoDataFrame(\n",
    "    index=[0],\n",
    "    crs=4326,\n",
    "    geometry=[\n",
    "        Polygon(\n",
    "            zip(\n",
    "                [MIN_LNG, MAX_LNG, MAX_LNG, MIN_LNG],\n",
    "                [MIN_LAT, MIN_LAT, MAX_LAT, MAX_LAT],\n",
    "            )\n",
    "        )\n",
    "    ],\n",
    ")\n",
    "\n",
    "tessellation = tilers.tiler.get(\n",
    "    TILE_TYPE,\n",
    "    base_shape=base_shape,\n",
    "    meters=TILE_DIAMETER_IN_METERS,\n",
    ")#.drop('H3_INDEX', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8259521f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_points(gdf, tessellation):\n",
    "    sp = gdf[['PERSON_ID', 'TRIP_ID', 'TRIP_START', 'START_POINT', 'clustering_HL']].copy()\n",
    "    ep = gdf[['PERSON_ID', 'TRIP_ID', 'TRIP_END', 'END_POINT', 'clustering_HL']].copy()\n",
    "\n",
    "\n",
    "    sp['lat'] = sp.START_POINT.apply(lambda x: x.y)\n",
    "    sp['lng'] = sp.START_POINT.apply(lambda x: x.x)\n",
    "    ep['lat'] = ep.END_POINT.apply(lambda x: x.y)\n",
    "    ep['lng'] = ep.END_POINT.apply(lambda x: x.x)\n",
    "\n",
    "    sp = sp.rename(columns={'TRIP_START': 'datetime'})\n",
    "    ep = ep.rename(columns={'TRIP_END': 'datetime'})\n",
    "\n",
    "    sp.drop('START_POINT', axis=1, inplace=True)\n",
    "    ep.drop('END_POINT', axis=1, inplace=True)\n",
    "\n",
    "    points = pd.concat([sp, ep])\n",
    "    points = gp.GeoDataFrame(points, geometry=gp.points_from_xy(points.lng, points.lat, crs='epsg:4326'))\n",
    "    \n",
    "    points = gp.sjoin(points, tessellation, predicate='within', how='left').drop('index_right', axis=1)\n",
    "    \n",
    "    return points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "23a4dd5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_location_entropy(tile_id, mapped_points_gdf):\n",
    "    assert isinstance(tile_id, str)\n",
    "    t_trips = mapped_points_gdf.query(\"tile_ID == @tile_id\")\n",
    "    le = 0\n",
    "    c_l = len(t_trips)\n",
    "    for p in t_trips.clustering_HL:\n",
    "        c_l_u = len(t_trips.query('clustering_HL == @p'))\n",
    "        p_u_l =  c_l_u/ c_l\n",
    "        le += p_u_l * math.log(p_u_l)\n",
    "    return -le"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2befdb64",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_inverse_loc_frequency(row, loc_frequencies, points):\n",
    "    \n",
    "    tl_ID = str(row.tile_ID)\n",
    "    cl_ID = row.clustering_HL\n",
    "    \n",
    "    cl_loc_freqs = loc_frequencies.query('clustering_HL == @cl_ID')\n",
    "        \n",
    "    tf = row.freq/cl_loc_freqs.freq.sum()\n",
    "    \n",
    "    idf = math.log(points.clustering_HL.nunique()/(points.query('tile_ID == @tl_ID').clustering_HL.nunique() + 1))\n",
    "    \n",
    "    tf_idf = tf * idf\n",
    "    \n",
    "    return tf_idf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b6283a89",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_location_similarity(cl_i, cl_j, loc_frequencies, points):\n",
    "    overlap_tile_ids = set(loc_frequencies.query(\"clustering_HL == @cl_i\").tile_ID).intersection(set(loc_frequencies.query(\"clustering_HL == @cl_j\").tile_ID))\n",
    "    nr_overlapping_tiles = len(overlap_tile_ids)\n",
    "    loc_sim = 0\n",
    "    \n",
    "    if nr_overlapping_tiles == 0:\n",
    "        return loc_sim\n",
    "    for t in overlap_tile_ids:\n",
    "        cl_i_freq = loc_frequencies[(loc_frequencies['clustering_HL'] == cl_i) & (loc_frequencies['tile_ID'] == t)].tf_idf.iloc[0]\n",
    "        cl_j_freq = loc_frequencies[(loc_frequencies['clustering_HL'] == cl_j) & (loc_frequencies['tile_ID'] == t)].tf_idf.iloc[0]\n",
    "        \n",
    "        # Alternative 1: Entropy times avg tfidfs\n",
    "#         loc_sim += points.query(\"tile_ID == @t\").location_entropy.iloc[0] * (cl_i_freq + cl_j_freq)/2\n",
    "\n",
    "        # Alternative 2: Interaction between tfidfs\n",
    "        loc_sim += cl_i_freq * cl_j_freq\n",
    "    return loc_sim/nr_overlapping_tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8f5561fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_max_n(M, n=10):\n",
    "    # (score, i, j)\n",
    "    max_n = [(0,0,0)] * n\n",
    "    \n",
    "    for i in range(M.shape[0]):\n",
    "        for j in range(M.shape[1]):\n",
    "            if not all(x[0] > M[i,j] for x in max_n):\n",
    "                index_smallest_val = max_n.index(sorted(max_n, reverse=False)[0])\n",
    "                max_n[index_smallest_val] = (M[i,j], i, j) # score must be first for sorting later\n",
    "    \n",
    "    print('Nr of links that are equal to one of max values: ', sum([sum(M == max_val[0]).sum() for max_val in max_n]))\n",
    "    return max_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6d0af3a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_nr_overlaps(cl_i, cl_j, gdf):\n",
    "    nr_overlaps = 0\n",
    "    \n",
    "    gdf['TRIP_START'] = pd.to_datetime(gdf['TRIP_START'], format='%Y-%m-%d %H:%M:%S')\n",
    "    gdf['TRIP_END'] = pd.to_datetime(gdf['TRIP_END'], format='%Y-%m-%d %H:%M:%S')\n",
    "    \n",
    "    cl_i_trips = gdf.query(\"clustering_HL == @cl_i\").copy()\n",
    "    cl_j_trips = gdf.query(\"clustering_HL == @cl_j\").copy()\n",
    "    \n",
    "    for i, t_i in cl_i_trips.iterrows():\n",
    "        t_i_s_x = t_i\n",
    "        t_i_e_x = pd.to_datetime(t_i['TRIP_END'], format='%Y-%m-%d %H:%M:%S')\n",
    "        \n",
    "        for i, t_j in cl_j_trips.iterrows():\n",
    "            if ((t_i.TRIP_START <= t_j.TRIP_START and t_i.TRIP_END >= t_j.TRIP_START) or (t_j.TRIP_START <= t_i.TRIP_START and t_j.TRIP_END >= t_i.TRIP_START) or (t_i.TRIP_START <= t_j.TRIP_START and t_i.TRIP_END >= t_j.TRIP_END) or (t_j.TRIP_START <= t_i.TRIP_START and t_j.TRIP_START >= t_i.TRIP_START)) and (t_i.TRIP_ID != t_j.TRIP_ID):\n",
    "                print('Overlapping trips in two clusters:', cl_i, cl_j, 'Trips: ', t_i.TRIP_ID, t_j.TRIP_ID)\n",
    "                nr_overlaps += 1\n",
    "                \n",
    "    return nr_overlaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8f0adacd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def store_results_recursive_search(clustering, gdf, best_match_locsim):\n",
    "    ground_truth = attack.getGroundTruth(gdf)\n",
    "    \n",
    "    # Write all clustering metrics of evaluate() to csv and add columns for parameters\n",
    "    result_dicts = []\n",
    "    result_dict = {}\n",
    "    result_dict['Homogeneity'] = metrics.homogeneity_score(ground_truth, clustering)\n",
    "    result_dict['Completeness'] = metrics.completeness_score(ground_truth, clustering)\n",
    "    result_dict['V-measure'] = metrics.v_measure_score(ground_truth, clustering)\n",
    "    result_dict['Rand index'] = metrics.rand_score(ground_truth, clustering)\n",
    "    result_dict['ARI'] = metrics.adjusted_rand_score(ground_truth, clustering)\n",
    "    result_dict['MI'] = metrics.mutual_info_score(ground_truth, clustering)\n",
    "    result_dict['NMI'] = metrics.normalized_mutual_info_score(ground_truth, clustering)\n",
    "    result_dict['AMI'] = metrics.adjusted_mutual_info_score(ground_truth, clustering)\n",
    "    result_dict['Cluster accuracy'] = attack.cluster_acc(ground_truth, clustering)\n",
    "    result_dict['Best match locsim'] = best_match_locsim\n",
    "    result_dicts.append(result_dict)\n",
    "\n",
    "    df = pd.DataFrame(result_dicts)\n",
    "\n",
    "    # Add column with date and time\n",
    "    df['Date'] = datetime.datetime.now().strftime(\"%Y-%m-%d\")\n",
    "    df['Time'] = datetime.datetime.now().strftime(\"%H:%M:%S\")\n",
    "\n",
    "    # Check if file exists\n",
    "    file_exists = os.path.isfile('results_recursive_search_geolife_500tessellation.csv')\n",
    "\n",
    "    # Write to csv (append)\n",
    "    if not file_exists:\n",
    "        df.to_csv('results_recursive_search_geolife_500tessellation.csv', mode='a', header=True, index=False)\n",
    "    else:\n",
    "        df.to_csv('results_recursive_search_geolife_500tessellation.csv', mode='a', header=False, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "430c6a7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_clustering_result(gdf, data_name):\n",
    "    # Write the cluster to a geojson file\n",
    "    gdf.drop(['START_POINT', 'END_POINT'], axis=1).to_file(f'data/geolife/{data_name}_clustering_tfidf.geojson', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45a8fe7c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGzCAYAAADJ3dZzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3yUlEQVR4nO3de3hU1b3/8U/uCYEkBEkmKRBSBSFcFTSMYq2SEiDliEQrNkJUjlSaIBBRSA9ylyCtiGgg2iLghYPSI7aCXEJAqBJuQRTBAlYkIJmEFpOQKLnO7w9/TDsmSBgm7GH7fj3Pfh5nrTV7f/d+FD6uvfYeL7vdbhcAAIBJeRtdAAAAQHMi7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AAAAFMj7AA/UsXFxbrnnnvUpk0beXl5aeHChZe9zxkzZsjLy8uprba2Vk8++aTat28vb29vDRs27JL327FjRz344IOOz++//768vLz0/vvvX17BTdDYOXl5eSk9Pb3Zjy1Jy5cvl5eXl7788ssrcjzAjAg7gEl4eXk1aTsfECZOnKiNGzcqMzNTr732mgYNGtTofs8Hiz//+c8u1fXKK6/o97//ve655x6tWLFCEydOdPUUL9vcuXP1zjvvGHb8H+LJtQFXO1+jCwDgHq+99prT51dffVW5ubkN2rt27SpJ2rJli+666y5NmjSpWevasmWLfvKTn+i5555z2z5/9rOf6dtvv5W/v/8lfW/u3Lm65557Lml2aerUqZoyZcolVnjpLlTbyJEjNWLECAUEBDR7DYBZEXYAk3jggQecPu/cuVO5ubkN2s8rKSlRWFhYs9fVHMfx9vZWYGCgW/f5fZWVlQoODpavr698fY37o9LHx0c+Pj6GHR8wA25jAT8y59eA2O12ZWdnO25vXaoPPvhAN910kwIDA3XttdfqpZdecur/8ssv5eXlpa1bt+rgwYMNbqM1xm63a86cOWrXrp1atGihO+64QwcPHmwwrrE1O0ePHlVycrIsFosCAwPVrl07jRgxQmVlZZK+u81XWVmpFStWOGo5vw7o/LqcQ4cO6de//rVat26t/v37O/U15o033tD111+vwMBA9enTR9u3b3fqf/DBB9WxY8cG3/v+Pn+otgut2Vm8eLG6deumgIAARUdHKy0tTaWlpU5jfv7zn6t79+46dOiQ7rjjDrVo0UI/+clPNH/+/EbPBzArZnaAH5mf/exneu211zRy5Ej94he/0KhRoy55HwcOHNDAgQPVtm1bzZgxQ7W1tZo+fboiIyMdY9q2bavXXntNTz/9tCoqKpSVlSXp37fRGjNt2jTNmTNHQ4YM0ZAhQ7Rv3z4NHDhQ1dXVP1hPdXW1EhMTVVVVpXHjxsliseirr77S2rVrVVpaqtDQUL322mv67//+b918880aM2aMJOnaa6912s+9996rTp06ae7cubLb7T94zG3btunNN9/UY489poCAAC1evFiDBg3S7t271b179x/87vc1pbb/NGPGDM2cOVMJCQkaO3asDh8+rCVLlmjPnj368MMP5efn5xj79ddfa9CgQRo+fLh+9atf6c9//rMmT56sHj16aPDgwZdUJ3DVsgMwpbS0NPsP/ScuyZ6WlnbR/WzdutUuyb569WpH27Bhw+yBgYH248ePO9oOHTpk9/HxaXDM22+/3d6tW7eLHqekpMTu7+9vT0pKstfX1zvaf/e739kl2VNTUxvUtHXrVrvdbrd/9NFHDWpsTHBwsNN+zps+fbpdkv3++++/YN9/kmSXZN+7d6+j7fjx4/bAwED73Xff7WhLTU21x8TENGmfF6pt2bJldkn2Y8eO2e32f1+ngQMH2uvq6hzjXnzxRbsk+yuvvOJou/322+2S7K+++qqjraqqym6xWOzJyckNjgWYFbexAFySuro6bdy4UcOGDVOHDh0c7V27dlViYqLL+928ebOqq6s1btw4p1s8EyZMuOh3Q0NDJUkbN27UN99843INjz76aJPHWq1W9enTx/G5Q4cOuuuuu7Rx40bV1dW5XMPFnL9OEyZMkLf3v/8If+SRRxQSEqJ169Y5jW/ZsqXTui1/f3/dfPPN+uKLL5qtRsDTEHYASJJsNpvT9u233zY67vTp0/r222/VqVOnBn3XX3+9y8c/fvy4JDXYb9u2bdW6desf/G5sbKwyMjL0pz/9Sddcc40SExOVnZ3tWK/TVLGxsU0e29j5d+7cWd98841Onz59Sce9FOev0/evtb+/v3760586+s9r165dgzVHrVu31tdff91sNQKehrADQJIUFRXltL355ptGl3RJnn32WX3yySf63e9+p2+//VaPPfaYunXrppMnTzZ5H0FBQW6t6UILm5tz5uf7LvQkl/0ia5IAMyHsAJAk5ebmOm0XuiXVtm1bBQUF6ejRow36Dh8+7PLxY2JiJKnBfk+fPt3kWYgePXpo6tSp2r59u/72t7/pq6++Uk5OjqPflafOLqSx8z9y5IhatGihtm3bSvpuBuX7T0hJajD7cim1nb9O37/W1dXVOnbsmKMfwL8RdgBIkhISEpy2qKioRsf5+PgoMTFR77zzjgoLCx3tn332mTZu3HhZx/fz89MLL7zgNOvQlJ+xKC8vV21trVNbjx495O3traqqKkdbcHBwo+HDFfn5+dq3b5/j84kTJ/SXv/xFAwcOdMymXHvttSorK9Mnn3ziGFdUVKQ1a9Y02F9Ta0tISJC/v78WLVrkdJ2WLl2qsrIyJSUlXcZZAebEo+cALtnMmTO1YcMG3Xbbbfrtb3+r2tpavfDCC+rWrZvTX+yXom3btpo0aZKysrL0y1/+UkOGDNFHH32k9evX65prrvnB727ZskXp6em699571blzZ9XW1uq1116Tj4+PkpOTHeP69OmjzZs3a8GCBYqOjlZsbKzi4+Ndqrd79+5KTEx0evRc+u7anDdixAhNnjxZd999tx577DF98803WrJkiTp37uwUlC6ltrZt2yozM1MzZ87UoEGD9F//9V86fPiwFi9erJtuuumCL5EEfswIOwAuWc+ePbVx40ZlZGRo2rRpateunWbOnKmioiKXw44kzZkzR4GBgcrJydHWrVsVHx+vTZs2XXS2olevXkpMTNS7776rr776Si1atFCvXr20fv169evXzzFuwYIFGjNmjKZOnapvv/1WqampLoed22+/XVarVTNnzlRhYaHi4uK0fPly9ezZ0zGmTZs2WrNmjTIyMvTkk08qNjZWWVlZOnr0aIOwcym1zZgxQ23bttWLL76oiRMnKjw8XGPGjNHcuXOd3rED4DtedlapAQAAE2PNDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXesyOpvr5ep06dUqtWrdz6OnkAANB87Ha7zp49q+joaHl7X3j+hrAj6dSpU2rfvr3RZQAAABecOHFC7dq1u2A/YUdSq1atJH13sUJCQgyuBgAANEV5ebnat2/v+Hv8Qgg7+vevDYeEhBB2AAC4ylxsCQoLlAEAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKn5Gl0APE/HKesuOubLeUlXoBIAAC4fMzsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUDA07dXV1euqppxQbG6ugoCBde+21mj17tux2u2OM3W7XtGnTFBUVpaCgICUkJOjo0aNO+zlz5oxSUlIUEhKisLAwjR49WhUVFVf6dAAAgAcyNOw888wzWrJkiV588UV99tlneuaZZzR//ny98MILjjHz58/XokWLlJOTo127dik4OFiJiYk6d+6cY0xKSooOHjyo3NxcrV27Vtu3b9eYMWOMOCUAAOBhvOz/OY1yhf3yl79UZGSkli5d6mhLTk5WUFCQXn/9ddntdkVHR+vxxx/XpEmTJEllZWWKjIzU8uXLNWLECH322WeKi4vTnj171LdvX0nShg0bNGTIEJ08eVLR0dEXraO8vFyhoaEqKytTSEhI85zsVaTjlHUXHfPlvKQrUAkAABfW1L+/DZ3ZueWWW5SXl6cjR45Ikj7++GN98MEHGjx4sCTp2LFjstlsSkhIcHwnNDRU8fHxys/PlyTl5+crLCzMEXQkKSEhQd7e3tq1a1ejx62qqlJ5ebnTBgAAzMnXyINPmTJF5eXl6tKli3x8fFRXV6enn35aKSkpkiSbzSZJioyMdPpeZGSko89msykiIsKp39fXV+Hh4Y4x35eVlaWZM2e6+3QAAIAHMnRm56233tIbb7yhlStXat++fVqxYoX+8Ic/aMWKFc163MzMTJWVlTm2EydONOvxAACAcQyd2XniiSc0ZcoUjRgxQpLUo0cPHT9+XFlZWUpNTZXFYpEkFRcXKyoqyvG94uJi9e7dW5JksVhUUlLitN/a2lqdOXPG8f3vCwgIUEBAQDOcEQAA8DSGzux888038vZ2LsHHx0f19fWSpNjYWFksFuXl5Tn6y8vLtWvXLlmtVkmS1WpVaWmpCgoKHGO2bNmi+vp6xcfHX4GzAAAAnszQmZ2hQ4fq6aefVocOHdStWzd99NFHWrBggR5++GFJkpeXlyZMmKA5c+aoU6dOio2N1VNPPaXo6GgNGzZMktS1a1cNGjRIjzzyiHJyclRTU6P09HSNGDGiSU9iAQAAczM07Lzwwgt66qmn9Nvf/lYlJSWKjo7Wb37zG02bNs0x5sknn1RlZaXGjBmj0tJS9e/fXxs2bFBgYKBjzBtvvKH09HQNGDBA3t7eSk5O1qJFi4w4JQAA4GEMfc+Op+A9O854zw4A4GpwVbxnBwAAoLkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkZGnY6duwoLy+vBltaWpok6dy5c0pLS1ObNm3UsmVLJScnq7i42GkfhYWFSkpKUosWLRQREaEnnnhCtbW1RpwOAADwQIaGnT179qioqMix5ebmSpLuvfdeSdLEiRP17rvvavXq1dq2bZtOnTql4cOHO75fV1enpKQkVVdXa8eOHVqxYoWWL1+uadOmGXI+AADA83jZ7Xa70UWcN2HCBK1du1ZHjx5VeXm52rZtq5UrV+qee+6RJP39739X165dlZ+fr379+mn9+vX65S9/qVOnTikyMlKSlJOTo8mTJ+v06dPy9/dv0nHLy8sVGhqqsrIyhYSENNv5XS06Tll30TFfzku6ApUAAHBhTf3722PW7FRXV+v111/Xww8/LC8vLxUUFKimpkYJCQmOMV26dFGHDh2Un58vScrPz1ePHj0cQUeSEhMTVV5eroMHD17wWFVVVSovL3faAACAOXlM2HnnnXdUWlqqBx98UJJks9nk7++vsLAwp3GRkZGy2WyOMf8ZdM73n++7kKysLIWGhjq29u3bu+9EAACAR/GYsLN06VINHjxY0dHRzX6szMxMlZWVObYTJ040+zEBAIAxfI0uQJKOHz+uzZs36+2333a0WSwWVVdXq7S01Gl2p7i4WBaLxTFm9+7dTvs6/7TW+TGNCQgIUEBAgBvPAAAAeCqPmNlZtmyZIiIilJT070Wvffr0kZ+fn/Ly8hxthw8fVmFhoaxWqyTJarXqwIEDKikpcYzJzc1VSEiI4uLirtwJAAAAj2X4zE59fb2WLVum1NRU+fr+u5zQ0FCNHj1aGRkZCg8PV0hIiMaNGyer1ap+/fpJkgYOHKi4uDiNHDlS8+fPl81m09SpU5WWlsbMDQAAkOQBYWfz5s0qLCzUww8/3KDvueeek7e3t5KTk1VVVaXExEQtXrzY0e/j46O1a9dq7NixslqtCg4OVmpqqmbNmnUlTwEAAHgwj3rPjlF4z44z3rMDALgaXHXv2QEAAGgOhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhB0AAGBqhoedr776Sg888IDatGmjoKAg9ejRQ3v37nX02+12TZs2TVFRUQoKClJCQoKOHj3qtI8zZ84oJSVFISEhCgsL0+jRo1VRUXGlTwUAAHggQ8PO119/rVtvvVV+fn5av369Dh06pGeffVatW7d2jJk/f74WLVqknJwc7dq1S8HBwUpMTNS5c+ccY1JSUnTw4EHl5uZq7dq12r59u8aMGWPEKQEAAA/jZbfb7UYdfMqUKfrwww/1t7/9rdF+u92u6OhoPf7445o0aZIkqaysTJGRkVq+fLlGjBihzz77THFxcdqzZ4/69u0rSdqwYYOGDBmikydPKjo6+qJ1lJeXKzQ0VGVlZQoJCXHfCV6lOk5Zd9ExX85LugKVAABwYU39+9vQmZ2//vWv6tu3r+69915FRETohhtu0B//+EdH/7Fjx2Sz2ZSQkOBoCw0NVXx8vPLz8yVJ+fn5CgsLcwQdSUpISJC3t7d27drV6HGrqqpUXl7utAEAAHMyNOx88cUXWrJkiTp16qSNGzdq7Nixeuyxx7RixQpJks1mkyRFRkY6fS8yMtLRZ7PZFBER4dTv6+ur8PBwx5jvy8rKUmhoqGNr3769u08NAAB4CEPDTn19vW688UbNnTtXN9xwg8aMGaNHHnlEOTk5zXrczMxMlZWVObYTJ0406/EAAIBxDA07UVFRiouLc2rr2rWrCgsLJUkWi0WSVFxc7DSmuLjY0WexWFRSUuLUX1tbqzNnzjjGfF9AQIBCQkKcNgAAYE6Ghp1bb71Vhw8fdmo7cuSIYmJiJEmxsbGyWCzKy8tz9JeXl2vXrl2yWq2SJKvVqtLSUhUUFDjGbNmyRfX19YqPj78CZwEAADyZr5EHnzhxom655RbNnTtXv/rVr7R79269/PLLevnllyVJXl5emjBhgubMmaNOnTopNjZWTz31lKKjozVs2DBJ380EDRo0yHH7q6amRunp6RoxYkSTnsQCAADmZmjYuemmm7RmzRplZmZq1qxZio2N1cKFC5WSkuIY8+STT6qyslJjxoxRaWmp+vfvrw0bNigwMNAx5o033lB6eroGDBggb29vJScna9GiRUacEgAA8DCGvmfHU/CeHWe8ZwcAcDW4Kt6zAwAA0NwIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNR8jS4AcIeOU9ZddMyX85KuQCUAAE9D2Glm/CUMAICxuI0FAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADAABMzaWw88UXX7i7DgAAgGbhUti57rrrdMcdd+j111/XuXPn3F0TAACA27gUdvbt26eePXsqIyNDFotFv/nNb7R792531wYAAHDZXAo7vXv31vPPP69Tp07plVdeUVFRkfr376/u3btrwYIFOn36tLvrBAAAcMllLVD29fXV8OHDtXr1aj3zzDP6/PPPNWnSJLVv316jRo1SUVGRu+oEAABwyWWFnb179+q3v/2toqKitGDBAk2aNEn/+Mc/lJubq1OnTumuu+5yV50AAAAucSnsLFiwQD169NAtt9yiU6dO6dVXX9Xx48c1Z84cxcbG6rbbbtPy5cu1b9++H9zPjBkz5OXl5bR16dLF0X/u3DmlpaWpTZs2atmypZKTk1VcXOy0j8LCQiUlJalFixaKiIjQE088odraWldOCwAAmJCvK19asmSJHn74YT344IOKiopqdExERISWLl160X1169ZNmzdv/ndBvv8uaeLEiVq3bp1Wr16t0NBQpaena/jw4frwww8lSXV1dUpKSpLFYtGOHTtUVFSkUaNGyc/PT3PnznXl1AAAgMm4FHaOHj160TH+/v5KTU29eAG+vrJYLA3ay8rKtHTpUq1cuVJ33nmnJGnZsmXq2rWrdu7cqX79+mnTpk06dOiQNm/erMjISPXu3VuzZ8/W5MmTNWPGDPn7+1/6yQEAAFNx6TbWsmXLtHr16gbtq1ev1ooVKy5pX0ePHlV0dLR++tOfKiUlRYWFhZKkgoIC1dTUKCEhwTG2S5cu6tChg/Lz8yVJ+fn56tGjhyIjIx1jEhMTVV5eroMHD17wmFVVVSovL3faAACAObkUdrKysnTNNdc0aI+IiLik20fx8fFavny5NmzYoCVLlujYsWO67bbbdPbsWdlsNvn7+yssLMzpO5GRkbLZbJIkm83mFHTO95/v+6H6Q0NDHVv79u2bXDMAALi6uHQbq7CwULGxsQ3aY2JiHDMzTTF48GDHP/fs2VPx8fGKiYnRW2+9paCgIFdKa5LMzExlZGQ4PpeXlxN4AAAwKZdmdiIiIvTJJ580aP/444/Vpk0bl4sJCwtT586d9fnnn8tisai6ulqlpaVOY4qLix1rfCwWS4Ons85/bmwd0HkBAQEKCQlx2gAAgDm5FHbuv/9+PfbYY9q6davq6upUV1enLVu2aPz48RoxYoTLxVRUVOgf//iHoqKi1KdPH/n5+SkvL8/Rf/jwYRUWFspqtUqSrFarDhw4oJKSEseY3NxchYSEKC4uzuU6AACAebh0G2v27Nn68ssvNWDAAMej4vX19Ro1atQlrdmZNGmShg4dqpiYGJ06dUrTp0+Xj4+P7r//foWGhmr06NHKyMhQeHi4QkJCNG7cOFmtVvXr10+SNHDgQMXFxWnkyJGaP3++bDabpk6dqrS0NAUEBLhyagAAwGRcCjv+/v568803NXv2bH388ccKCgpSjx49FBMTc0n7OXnypO6//37961//Utu2bdW/f3/t3LlTbdu2lSQ999xz8vb2VnJysqqqqpSYmKjFixc7vu/j46O1a9dq7NixslqtCg4OVmpqqmbNmuXKaQEAABNyKeyc17lzZ3Xu3Nnl769ateoH+wMDA5Wdna3s7OwLjomJidF7773ncg0AAMDcXAo7dXV1Wr58ufLy8lRSUqL6+nqn/i1btrilOAAAgMvlUtgZP368li9frqSkJHXv3l1eXl7urgsAAMAtXAo7q1at0ltvvaUhQ4a4ux4AAAC3cunRc39/f1133XXurgUAAMDtXAo7jz/+uJ5//nnZ7XZ31wMAAOBWLt3G+uCDD7R161atX79e3bp1k5+fn1P/22+/7ZbiAAAALpdLYScsLEx33323u2sBAABwO5fCzrJly9xdBwAAQLNwac2OJNXW1mrz5s166aWXdPbsWUnSqVOnVFFR4bbiAAAALpdLMzvHjx/XoEGDVFhYqKqqKv3iF79Qq1at9Mwzz6iqqko5OTnurhMAAMAlLs3sjB8/Xn379tXXX3+toKAgR/vdd9/t9CvlAAAARnNpZudvf/ubduzYIX9/f6f2jh076quvvnJLYQAAAO7g0sxOfX296urqGrSfPHlSrVq1uuyiAAAA3MWlsDNw4EAtXLjQ8dnLy0sVFRWaPn06PyEBAAA8iku3sZ599lklJiYqLi5O586d069//WsdPXpU11xzjf73f//X3TUCAAC4zKWw065dO3388cdatWqVPvnkE1VUVGj06NFKSUlxWrAMAABgNJfCjiT5+vrqgQcecGctAAAAbudS2Hn11Vd/sH/UqFEuFQMAAOBuLoWd8ePHO32uqanRN998I39/f7Vo0YKwAwAAPIZLT2N9/fXXTltFRYUOHz6s/v37s0AZAAB4FJd/G+v7OnXqpHnz5jWY9QEAADCS28KO9N2i5VOnTrlzlwAAAJfFpTU7f/3rX50+2+12FRUV6cUXX9Stt97qlsIAAADcwaWwM2zYMKfPXl5eatu2re688049++yz7qgLAADALVwKO/X19e6uAwAAoFm4dc0OAACAp3FpZicjI6PJYxcsWODKIQAAANzCpbDz0Ucf6aOPPlJNTY2uv/56SdKRI0fk4+OjG2+80THOy8vLPVUCAAC4yKWwM3ToULVq1UorVqxQ69atJX33osGHHnpIt912mx5//HG3FgkAAOAql9bsPPvss8rKynIEHUlq3bq15syZw9NYAADAo7gUdsrLy3X69OkG7adPn9bZs2cvuygAAAB3cSns3H333XrooYf09ttv6+TJkzp58qT+7//+T6NHj9bw4cNdKmTevHny8vLShAkTHG3nzp1TWlqa2rRpo5YtWyo5OVnFxcVO3yssLFRSUpJatGihiIgIPfHEE6qtrXWpBgAAYD4urdnJycnRpEmT9Otf/1o1NTXf7cjXV6NHj9bvf//7S97fnj179NJLL6lnz55O7RMnTtS6deu0evVqhYaGKj09XcOHD9eHH34oSaqrq1NSUpIsFot27NihoqIijRo1Sn5+fpo7d64rpwYAAEzGpZmdFi1aaPHixfrXv/7leDLrzJkzWrx4sYKDgy9pXxUVFUpJSdEf//hHpzVAZWVlWrp0qRYsWKA777xTffr00bJly7Rjxw7t3LlTkrRp0yYdOnRIr7/+unr37q3Bgwdr9uzZys7OVnV19QWPWVVVpfLycqcNAACY02W9VLCoqEhFRUXq1KmTgoODZbfbL3kfaWlpSkpKUkJCglN7QUGBampqnNq7dOmiDh06KD8/X5KUn5+vHj16KDIy0jEmMTFR5eXlOnjw4AWPmZWVpdDQUMfWvn37S64bAABcHVwKO//61780YMAAde7cWUOGDFFRUZEkafTo0Zf02PmqVau0b98+ZWVlNeiz2Wzy9/dXWFiYU3tkZKRsNptjzH8GnfP95/suJDMzU2VlZY7txIkTTa4ZAABcXVwKOxMnTpSfn58KCwvVokULR/t9992nDRs2NGkfJ06c0Pjx4/XGG28oMDDQlTJcFhAQoJCQEKcNAACYk0thZ9OmTXrmmWfUrl07p/ZOnTrp+PHjTdpHQUGBSkpKdOONN8rX11e+vr7atm2bFi1aJF9fX0VGRqq6ulqlpaVO3ysuLpbFYpEkWSyWBk9nnf98fgwAAPhxcynsVFZWOs3onHfmzBkFBAQ0aR8DBgzQgQMHtH//fsfWt29fpaSkOP7Zz89PeXl5ju8cPnxYhYWFslqtkiSr1aoDBw6opKTEMSY3N1chISGKi4tz5dQAAIDJuPTo+W233aZXX31Vs2fPlvTdb2DV19dr/vz5uuOOO5q0j1atWql79+5ObcHBwWrTpo2jffTo0crIyFB4eLhCQkI0btw4Wa1W9evXT5I0cOBAxcXFaeTIkZo/f75sNpumTp2qtLS0JocuAABgbi6Fnfnz52vAgAHau3evqqur9eSTT+rgwYM6c+aM4x047vDcc8/J29tbycnJqqqqUmJiohYvXuzo9/Hx0dq1azV27FhZrVYFBwcrNTVVs2bNclsNAADg6uZS2OnevbuOHDmiF198Ua1atVJFRYWGDx+utLQ0RUVFuVzM+++/7/Q5MDBQ2dnZys7OvuB3YmJi9N5777l8TAAAYG6XHHZqamo0aNAg5eTk6H/+53+aoyYAAAC3ueQFyn5+fvrkk0+aoxYAAAC3c+lprAceeEBLly51dy0AAABu59KandraWr3yyivavHmz+vTp0+D3sBYsWOCW4gAAAC7XJYWdL774Qh07dtSnn36qG2+8UZJ05MgRpzFeXl7uqw4AAOAyXVLY6dSpk4qKirR161ZJ3/08xKJFixr8PhUAAICnuKQ1O9//VfP169ersrLSrQUBAAC4k0trds77fvgB4F4dp6xr0rgv5yU1cyUAcPW6pJkdLy+vBmtyWKMDAAA82SXN7Njtdj344IOO3506d+6cHn300QZPY7399tvuqxC4gpoyk8IsCgBcXS4p7KSmpjp9fuCBB9xaDAAAgLtdUthZtmxZc9UBAADQLFx6gzIAAMDVgrADAABM7bIePQfgGVhYDQAXxswOAAAwNcIOAAAwNcIOAAAwNcIOAAAwNcIOAAAwNcIOAAAwNR49B34keDwdwI8VMzsAAMDUCDsAAMDUCDsAAMDUCDsAAMDUWKAMj9eUhbUAAFwIYQcA3Iin3gDPw20sAABgaoQdAABgaoaGnSVLlqhnz54KCQlRSEiIrFar1q9f7+g/d+6c0tLS1KZNG7Vs2VLJyckqLi522kdhYaGSkpLUokULRURE6IknnlBtbe2VPhUAAOChDF2z065dO82bN0+dOnWS3W7XihUrdNddd+mjjz5St27dNHHiRK1bt06rV69WaGio0tPTNXz4cH344YeSpLq6OiUlJclisWjHjh0qKirSqFGj5Ofnp7lz5xp5avBAnrbQ2dPqAQCzMjTsDB061Onz008/rSVLlmjnzp1q166dli5dqpUrV+rOO++UJC1btkxdu3bVzp071a9fP23atEmHDh3S5s2bFRkZqd69e2v27NmaPHmyZsyYIX9/fyNOCwAAeBCPWbNTV1enVatWqbKyUlarVQUFBaqpqVFCQoJjTJcuXdShQwfl5+dLkvLz89WjRw9FRkY6xiQmJqq8vFwHDx684LGqqqpUXl7utAEAAHMyPOwcOHBALVu2VEBAgB599FGtWbNGcXFxstls8vf3V1hYmNP4yMhI2Ww2SZLNZnMKOuf7z/ddSFZWlkJDQx1b+/bt3XtSAADAYxgedq6//nrt379fu3bt0tixY5WamqpDhw416zEzMzNVVlbm2E6cONGsxwMAAMYx/KWC/v7+uu666yRJffr00Z49e/T888/rvvvuU3V1tUpLS51md4qLi2WxWCRJFotFu3fvdtrf+ae1zo9pTEBAgAICAtx8JgAAwBMZPrPzffX19aqqqlKfPn3k5+envLw8R9/hw4dVWFgoq9UqSbJarTpw4IBKSkocY3JzcxUSEqK4uLgrXjsAAPA8hs7sZGZmavDgwerQoYPOnj2rlStX6v3339fGjRsVGhqq0aNHKyMjQ+Hh4QoJCdG4ceNktVrVr18/SdLAgQMVFxenkSNHav78+bLZbJo6darS0tKYuQEAAJIMDjslJSUaNWqUioqKFBoaqp49e2rjxo36xS9+IUl67rnn5O3treTkZFVVVSkxMVGLFy92fN/Hx0dr167V2LFjZbVaFRwcrNTUVM2aNcuoUwIAAB7G0LCzdOnSH+wPDAxUdna2srOzLzgmJiZG7733nrtLAwAAJuFxa3YAAADcibADAABMjbADAABMzfD37AAwn6b8yOmX85KuQCUAwMwOAAAwOcIOAAAwNcIOAAAwNcIOAAAwNcIOAAAwNcIOAAAwNcIOAAAwNd6zA1wi3iEDAFcXwg6AS9KUsAcAnoTbWAAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNR49ByAw4/5sXLenwSYFzM7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1Ag7AADA1AwNO1lZWbrpppvUqlUrRUREaNiwYTp8+LDTmHPnziktLU1t2rRRy5YtlZycrOLiYqcxhYWFSkpKUosWLRQREaEnnnhCtbW1V/JUAACAhzI07Gzbtk1paWnauXOncnNzVVNTo4EDB6qystIxZuLEiXr33Xe1evVqbdu2TadOndLw4cMd/XV1dUpKSlJ1dbV27NihFStWaPny5Zo2bZoRpwQAADyMob96vmHDBqfPy5cvV0REhAoKCvSzn/1MZWVlWrp0qVauXKk777xTkrRs2TJ17dpVO3fuVL9+/bRp0yYdOnRImzdvVmRkpHr37q3Zs2dr8uTJmjFjhvz9/Rsct6qqSlVVVY7P5eXlzXuiAADAMB61ZqesrEySFB4eLkkqKChQTU2NEhISHGO6dOmiDh06KD8/X5KUn5+vHj16KDIy0jEmMTFR5eXlOnjwYKPHycrKUmhoqGNr3759c50SAAAwmMeEnfr6ek2YMEG33nqrunfvLkmy2Wzy9/dXWFiY09jIyEjZbDbHmP8MOuf7z/c1JjMzU2VlZY7txIkTbj4bAADgKQy9jfWf0tLS9Omnn+qDDz5o9mMFBAQoICCg2Y8DAACM5xFhJz09XWvXrtX27dvVrl07R7vFYlF1dbVKS0udZneKi4tlsVgcY3bv3u20v/NPa50fA8C8Ok5Zd9ExX85LugKVAPBUht7GstvtSk9P15o1a7RlyxbFxsY69ffp00d+fn7Ky8tztB0+fFiFhYWyWq2SJKvVqgMHDqikpMQxJjc3VyEhIYqLi7syJwIAADyWoTM7aWlpWrlypf7yl7+oVatWjjU2oaGhCgoKUmhoqEaPHq2MjAyFh4crJCRE48aNk9VqVb9+/SRJAwcOVFxcnEaOHKn58+fLZrNp6tSpSktL41YVAAAwNuwsWbJEkvTzn//cqX3ZsmV68MEHJUnPPfecvL29lZycrKqqKiUmJmrx4sWOsT4+Plq7dq3Gjh0rq9Wq4OBgpaamatasWVfqNAB4uKbc6gJgXoaGHbvdftExgYGBys7OVnZ29gXHxMTE6L333nNnaQAAwCQ85tFzAACA5kDYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApkbYAQAApuZrdAEAcCEdp6wzugQAJsDMDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDXCDgAAMDV+CBQuacoPNH45L+kKVAIAwA8j7MBQ/Ko1AKC5cRsLAACYmqFhZ/v27Ro6dKiio6Pl5eWld955x6nfbrdr2rRpioqKUlBQkBISEnT06FGnMWfOnFFKSopCQkIUFham0aNHq6Ki4gqeBQAA8GSGhp3Kykr16tVL2dnZjfbPnz9fixYtUk5Ojnbt2qXg4GAlJibq3LlzjjEpKSk6ePCgcnNztXbtWm3fvl1jxoy5UqcAAAA8nKFrdgYPHqzBgwc32me327Vw4UJNnTpVd911lyTp1VdfVWRkpN555x2NGDFCn332mTZs2KA9e/aob9++kqQXXnhBQ4YM0R/+8AdFR0dfsXMBAACeyWPX7Bw7dkw2m00JCQmOttDQUMXHxys/P1+SlJ+fr7CwMEfQkaSEhAR5e3tr165dF9x3VVWVysvLnTYAAGBOHht2bDabJCkyMtKpPTIy0tFns9kUERHh1O/r66vw8HDHmMZkZWUpNDTUsbVv397N1QMAAE/hsWGnOWVmZqqsrMyxnThxwuiSAABAM/HYsGOxWCRJxcXFTu3FxcWOPovFopKSEqf+2tpanTlzxjGmMQEBAQoJCXHaAACAOXls2ImNjZXFYlFeXp6jrby8XLt27ZLVapUkWa1WlZaWqqCgwDFmy5Ytqq+vV3x8/BWvGQAAeB5Dn8aqqKjQ559/7vh87Ngx7d+/X+Hh4erQoYMmTJigOXPmqFOnToqNjdVTTz2l6OhoDRs2TJLUtWtXDRo0SI888ohycnJUU1Oj9PR0jRgxgiexAACAJIPDzt69e3XHHXc4PmdkZEiSUlNTtXz5cj355JOqrKzUmDFjVFpaqv79+2vDhg0KDAx0fOeNN95Qenq6BgwYIG9vbyUnJ2vRokVX/FwuB78zBQBA8zE07Pz85z+X3W6/YL+Xl5dmzZqlWbNmXXBMeHi4Vq5c2Rzl4TLxu1f4Ifz7AeBK8dg1OwAAAO5A2AEAAKZG2AEAAKZm6JodwKxYjwIAnoOZHQAAYGqEHQAAYGqEHQAAYGqEHQAAYGqEHQAAYGo8jXWV4CclAABwDTM7AADA1Ag7AADA1LiNZSLc6gIAoCHCzo8Mb/YFAPzYcBsLAACYGmEHAACYGrexAOBHjvV+MDtmdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKkRdgAAgKnxcxEA4IH4CQfAfZjZAQAApsbMDgDALZiNgqcyzcxOdna2OnbsqMDAQMXHx2v37t1GlwQAADyAKWZ23nzzTWVkZCgnJ0fx8fFauHChEhMTdfjwYUVERBhdHgCTaMrMhae5GmtuCmaRcClMEXYWLFigRx55RA899JAkKScnR+vWrdMrr7yiKVOmGFwdAOA8QgqMcNWHnerqahUUFCgzM9PR5u3trYSEBOXn5zf6naqqKlVVVTk+l5WVSZLKy8vdXl991Tdu3yeAq1tT/qxpyp8d7tpPU3SYuNot+7mSx2rKfj6dmeiWYzVF9+kbLzrmStbTFJ5e8/n/Bux2+w+Ou+rDzj//+U/V1dUpMjLSqT0yMlJ///vfG/1OVlaWZs6c2aC9ffv2zVIjAPyn0IWetZ8fM0+7hp5WT1N4Qs1nz55VaGjoBfuv+rDjiszMTGVkZDg+19fX68yZM2rTpo28vLzcdpzy8nK1b99eJ06cUEhIiNv2e7XjujSO69IQ16RxXJfGcV0aZ+brYrfbdfbsWUVHR//guKs+7FxzzTXy8fFRcXGxU3txcbEsFkuj3wkICFBAQIBTW1hYWHOVqJCQENP9C+YOXJfGcV0a4po0juvSOK5L48x6XX5oRue8q/7Rc39/f/Xp00d5eXmOtvr6euXl5clqtRpYGQAA8ARX/cyOJGVkZCg1NVV9+/bVzTffrIULF6qystLxdBYAAPjxMkXYue+++3T69GlNmzZNNptNvXv31oYNGxosWr7SAgICNH369Aa3zH7suC6N47o0xDVpHNelcVyXxnFdJC/7xZ7XAgAAuIpd9Wt2AAAAfghhBwAAmBphBwAAmBphBwAAmBphBwAAmBphp5lkZ2erY8eOCgwMVHx8vHbv3m10SYbbvn27hg4dqujoaHl5eemdd94xuiTDZWVl6aabblKrVq0UERGhYcOG6fDhw0aXZbglS5aoZ8+ejje+Wq1WrV+/3uiyPM68efPk5eWlCRMmGF2KoWbMmCEvLy+nrUuXLkaX5RG++uorPfDAA2rTpo2CgoLUo0cP7d271+iyrjjCTjN48803lZGRoenTp2vfvn3q1auXEhMTVVJSYnRphqqsrFSvXr2UnZ1tdCkeY9u2bUpLS9POnTuVm5urmpoaDRw4UJWVlUaXZqh27dpp3rx5Kigo0N69e3XnnXfqrrvu0sGDB40uzWPs2bNHL730knr27Gl0KR6hW7duKioqcmwffPCB0SUZ7uuvv9att94qPz8/rV+/XocOHdKzzz6r1q1bG13alWeH29188832tLQ0x+e6ujp7dHS0PSsry8CqPIsk+5o1a4wuw+OUlJTYJdm3bdtmdCkep3Xr1vY//elPRpfhEc6ePWvv1KmTPTc313777bfbx48fb3RJhpo+fbq9V69eRpfhcSZPnmzv37+/0WV4BGZ23Ky6uloFBQVKSEhwtHl7eyshIUH5+fkGVoarQVlZmSQpPDzc4Eo8R11dnVatWqXKykp+7+7/S0tLU1JSktOfMz92R48eVXR0tH76058qJSVFhYWFRpdkuL/+9a/q27ev7r33XkVEROiGG27QH//4R6PLMgRhx83++c9/qq6ursFPVURGRspmsxlUFa4G9fX1mjBhgm699VZ1797d6HIMd+DAAbVs2VIBAQF69NFHtWbNGsXFxRldluFWrVqlffv2KSsry+hSPEZ8fLyWL1+uDRs2aMmSJTp27Jhuu+02nT171ujSDPXFF19oyZIl6tSpkzZu3KixY8fqscce04oVK4wu7YozxW9jAWaQlpamTz/9lLUG/9/111+v/fv3q6ysTH/+85+Vmpqqbdu2/agDz4kTJzR+/Hjl5uYqMDDQ6HI8xuDBgx3/3LNnT8XHxysmJkZvvfWWRo8ebWBlxqqvr1ffvn01d+5cSdINN9ygTz/9VDk5OUpNTTW4uiuLmR03u+aaa+Tj46Pi4mKn9uLiYlksFoOqgqdLT0/X2rVrtXXrVrVr187ocjyCv7+/rrvuOvXp00dZWVnq1auXnn/+eaPLMlRBQYFKSkp04403ytfXV76+vtq2bZsWLVokX19f1dXVGV2iRwgLC1Pnzp31+eefG12KoaKiohr8z0HXrl1/lLf4CDtu5u/vrz59+igvL8/RVl9fr7y8PNYboAG73a709HStWbNGW7ZsUWxsrNEleaz6+npVVVUZXYahBgwYoAMHDmj//v2OrW/fvkpJSdH+/fvl4+NjdIkeoaKiQv/4xz8UFRVldCmGuvXWWxu8yuLIkSOKiYkxqCLjcBurGWRkZCg1NVV9+/bVzTffrIULF6qyslIPPfSQ0aUZqqKiwun/tI4dO6b9+/crPDxcHTp0MLAy46SlpWnlypX6y1/+olatWjnWdYWGhiooKMjg6oyTmZmpwYMHq0OHDjp79qxWrlyp999/Xxs3bjS6NEO1atWqwXqu4OBgtWnT5ke9zmvSpEkaOnSoYmJidOrUKU2fPl0+Pj66//77jS7NUBMnTtQtt9yiuXPn6le/+pV2796tl19+WS+//LLRpV15Rj8OZlYvvPCCvUOHDnZ/f3/7zTffbN+5c6fRJRlu69atdkkNttTUVKNLM0xj10OSfdmyZUaXZqiHH37YHhMTY/f397e3bdvWPmDAAPumTZuMLssj8ei53X7ffffZo6Ki7P7+/vaf/OQn9vvuu8/++eefG12WR3j33Xft3bt3twcEBNi7dOlif/nll40uyRBedrvdblDOAgAAaHas2QEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKb2/wBVq9z6Ts1I7gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 29%|██████████████████████                                                      | 922/3183 [1:09:07<2:51:46,  4.56s/it]"
     ]
    }
   ],
   "source": [
    "iteration = 1\n",
    "\n",
    "# Define stopping criteria (either k is reached or min_sim is triggered)\n",
    "K_target = 73\n",
    "\n",
    "points = get_points(gdf, tessellation)\n",
    "loc_frequencies = points.groupby(['clustering_HL', 'tile_ID'])[['TRIP_ID']].nunique().reset_index().sort_values('clustering_HL').rename(columns={'TRIP_ID': 'freq'})\n",
    "loc_frequencies['tf_idf'] = loc_frequencies.apply(lambda x: get_inverse_loc_frequency(x, loc_frequencies, points), axis=1)\n",
    "Min_sim = loc_frequencies.tf_idf.quantile(0.75)**2\n",
    "print(f\"Min Sim threshold: {Min_sim}\")\n",
    "\n",
    "while True:\n",
    "    print('\\nNext Iteration: ', iteration)\n",
    "    # Get new points gdf matched with tessellation that contains updated clustering column\n",
    "    points = get_points(gdf, tessellation)\n",
    "    \n",
    "    # Compute location entropy for latest clustering (inverse normalized entropy)\n",
    "#     points['location_entropy'] = points.tile_ID.progress_apply(lambda x: get_location_entropy(x, points))\n",
    "#     points['location_entropy'] = points.location_entropy.apply(lambda x: 1 - (x - points.location_entropy.min())/(points.location_entropy.max() - points.location_entropy.min()))\n",
    "#     points.location_entropy.plot.hist(bins=30, title=\"Inverse location_entropy distribution\")\n",
    "#     plt.show()\n",
    "    \n",
    "    # Update location frequencies and loc frequency-inverse cluster frequency\n",
    "    loc_frequencies = points.groupby(['clustering_HL', 'tile_ID'])[['TRIP_ID']].nunique().reset_index().sort_values('clustering_HL').rename(columns={'TRIP_ID': 'freq'})\n",
    "    loc_frequencies['tf_idf'] = loc_frequencies.apply(lambda x: get_inverse_loc_frequency(x, loc_frequencies, points), axis=1)\n",
    "    loc_frequencies.tf_idf.plot.hist(bins=50, title=\"Tf-Idf distribution\")\n",
    "    plt.show()\n",
    "    \n",
    "    # Recalculate similarities between latest clusters\n",
    "    M = np.zeros((loc_frequencies.clustering_HL.nunique(), loc_frequencies.clustering_HL.nunique()))\n",
    "    for i, cl_i in tqdm(enumerate(loc_frequencies.clustering_HL.unique()), total=loc_frequencies.clustering_HL.nunique()):\n",
    "        for j, cl_j in enumerate(loc_frequencies.clustering_HL.unique()):\n",
    "            if i == j:\n",
    "                continue\n",
    "            M[i, j] = get_location_similarity(cl_i, cl_j, loc_frequencies, points)\n",
    "    np.fill_diagonal(M, 0)\n",
    "    \n",
    "    # Select n best links and assign new clustering IDs\n",
    "    print('New links of clusters:')\n",
    "    max_n =  get_max_n(M, 100)\n",
    "    # Adjust scores based on overlapping trips in the two candidate clusters\n",
    "    print(max_n)\n",
    "    max_n = [(t[0] * 1/(find_nr_overlaps(loc_frequencies.clustering_HL.unique()[t[1]], loc_frequencies.clustering_HL.unique()[t[2]], gdf) + 1), t[1], t[2]) for t in max_n]\n",
    "    print(max_n)\n",
    "    \n",
    "    # Take the best fitting clusters\n",
    "    for i, tup in enumerate(sorted(max_n, reverse=True)[:100]):\n",
    "        cl_1 = loc_frequencies.clustering_HL.unique()[tup[1]]\n",
    "        cl_2 = loc_frequencies.clustering_HL.unique()[tup[2]]\n",
    "        print(cl_1, cl_2)\n",
    "        # Assign new clustering ids (merge clustering ids)\n",
    "        gdf['clustering_HL'] = gdf.clustering_HL.apply(lambda x: cl_1 if x == cl_2 else x)\n",
    "    \n",
    "    # Evaluate result of iteration\n",
    "    attack.evaluate(gdf.clustering_HL.tolist(), gdf)\n",
    "    best_match_locsim = sorted(max_n, key=lambda x: x[0], reverse=True)[0][0]\n",
    "    store_results_recursive_search(gdf.clustering_HL.tolist(), gdf, best_match_locsim)\n",
    "    \n",
    "    print('Iteration', iteration, 'done. ', gdf.clustering_HL.nunique(), 'clusters left.')\n",
    "    \n",
    "    # Check break conditions\n",
    "    if gdf.clustering_HL.nunique() <= K_target:\n",
    "        print(\"K clusters reached! Done!\")\n",
    "        save_clustering_result(gdf, data_name=\"geolife_500tessellation\")\n",
    "        break\n",
    "    if sorted(max_n, reverse=True)[0][0] < Min_sim:\n",
    "        print(\"No good matches to make left! Done!\")\n",
    "        save_clustering_result(gdf, data_name=\"geolife_500tessellation\")\n",
    "        break\n",
    "                \n",
    "    iteration += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d8ebe1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.drop(['START_POINT', 'END_POINT'], axis=1).to_file(\"data/geolife/geolife_clustering_tfidf.geojson\", driver=\"GeoJSON\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}