Recursive_Search_geolife_private_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_private_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_private_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": 11,
   "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_private_500tessellation.csv')\n",
    "\n",
    "    # Write to csv (append)\n",
    "    if not file_exists:\n",
    "        df.to_csv('results_recursive_search_geolife_private_500tessellation.csv', mode='a', header=True, index=False)\n",
    "    else:\n",
    "        df.to_csv('results_recursive_search_geolife_private_500tessellation.csv', mode='a', header=False, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAy0klEQVR4nO3de1xVZb7H8e8GBBQF1OSWilSa4rU0dY92UxKV08mki2VF5eTJwVLJMuY4XkscO5lZqNWYqJPHsskuliJe0kpMpSzTXmRlYSrgGROEEgTW+aOXe2aH1+2GhQ+f9+u1Xq/Ws5691u/ZjfmdZz1rL4dlWZYAAAAM5WN3AQAAADWJsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wA9RTBQUFuu2229S8eXM5HA7NmTPngs85ZcoUORwOt7aKigo98cQTatWqlXx8fDRkyJDzPm+bNm10//33u/Y//PBDORwOffjhhxdW8Dk41ZgcDodGjx5d49eWpIyMDDkcDv3www+1cj3ARIQdwBAOh+OctpMBYdy4ccrMzFRqaqqWLl2qgQMHnvK8J4PFm2++6VFdr776qp555hnddtttWrx4scaNG+fpEC/YjBkz9Pbbb9t2/TOpy7UBFzs/uwsA4B1Lly5121+yZImysrKqtXfo0EGStGHDBt1yyy0aP358jda1YcMGXXrppXruuee8ds7rrrtOv/76q/z9/c/rczNmzNBtt912XrNLEydO1JNPPnmeFZ6/09V27733atiwYQoICKjxGgBTEXYAQ9xzzz1u+1u3blVWVla19pMKCwsVGhpa43XVxHV8fHwUGBjo1XP+XmlpqYKCguTn5yc/P/v+U+nr6ytfX1/brg+YgNtYQD1zcg2IZVlKT0933d46Xx9//LGuueYaBQYG6vLLL9dLL73kdvyHH36Qw+HQxo0btXv37mq30U7Fsiw99dRTatmypRo1aqQbb7xRu3fvrtbvVGt29u7dq8TEREVERCgwMFAtW7bUsGHDVFRUJOm323ylpaVavHixq5aT64BOrsvZs2eP7r77bjVt2lR9+/Z1O3Yqr732mq688koFBgaqe/fu2rx5s9vx+++/X23atKn2ud+f80y1nW7Nzrx589SxY0cFBAQoKipKycnJOnr0qFufG264QZ06ddKePXt04403qlGjRrr00ks1a9asU44HMBUzO0A9c91112np0qW69957ddNNN+m+++4773Ps2rVLAwYMUIsWLTRlyhRVVFRo8uTJCg8Pd/Vp0aKFli5dqqefflolJSVKS0uT9K/baKcyadIkPfXUUxo8eLAGDx6szz77TAMGDFB5efkZ6ykvL1d8fLzKysr0yCOPKCIiQgcOHNCqVat09OhRhYSEaOnSpfrjH/+onj17auTIkZKkyy+/3O08t99+u9q2basZM2bIsqwzXnPTpk16/fXX9eijjyogIEDz5s3TwIEDtW3bNnXq1OmMn/29c6nt302ZMkVTp05VXFycRo0apdzcXM2fP1/bt2/XJ598ogYNGrj6/vzzzxo4cKCGDh2qO+64Q2+++aYmTJigzp07a9CgQedVJ3DRsgAYKTk52TrTH3FJVnJy8lnPs3HjRkuStWLFClfbkCFDrMDAQOvHH390te3Zs8fy9fWtds3rr7/e6tix41mvU1hYaPn7+1sJCQlWVVWVq/3Pf/6zJclKSkqqVtPGjRsty7Kszz//vFqNpxIUFOR2npMmT55sSbLuuuuu0x77d5IsSdaOHTtcbT/++KMVGBho3Xrrra62pKQkKzo6+pzOebraFi1aZEmy9u3bZ1nWv76nAQMGWJWVla5+L774oiXJevXVV11t119/vSXJWrJkiautrKzMioiIsBITE6tdCzAVt7EAnJfKykplZmZqyJAhat26tau9Q4cOio+P9/i869atU3l5uR555BG3Wzxjx44962dDQkIkSZmZmfrll188ruHhhx8+575Op1Pdu3d37bdu3Vq33HKLMjMzVVlZ6XENZ3Pyexo7dqx8fP71n/CHHnpIwcHBev/99936N27c2G3dlr+/v3r27Knvv/++xmoE6hrCDgBJUn5+vtv266+/nrLf4cOH9euvv6pt27bVjl155ZUeX//HH3+UpGrnbdGihZo2bXrGz8bExCglJUV/+9vfdMkllyg+Pl7p6emu9TrnKiYm5pz7nmr87dq10y+//KLDhw+f13XPx8nv6ffftb+/vy677DLX8ZNatmxZbc1R06ZN9fPPP9dYjUBdQ9gBIEmKjIx0215//XW7Szovzz77rL788kv9+c9/1q+//qpHH31UHTt21E8//XTO52jYsKFXazrdwuaanPn5vdM9yWWdZU0SYBLCDgBJUlZWltt2ultSLVq0UMOGDbV3795qx3Jzcz2+fnR0tCRVO+/hw4fPeRaic+fOmjhxojZv3qyPPvpIBw4c0IIFC1zHPXnq7HRONf5vvvlGjRo1UosWLST9NoPy+yekJFWbfTmf2k5+T7//rsvLy7Vv3z7XcQD/QtgBIEmKi4tz2yIjI0/Zz9fXV/Hx8Xr77beVl5fnav/666+VmZl5Qddv0KCBXnjhBbdZh3N5jUVxcbEqKirc2jp37iwfHx+VlZW52oKCgk4ZPjyRnZ2tzz77zLW/f/9+vfPOOxowYIBrNuXyyy9XUVGRvvzyS1e/Q4cOaeXKldXOd661xcXFyd/fX3PnznX7nhYuXKiioiIlJCRcwKgAM/HoOYDzNnXqVK1Zs0bXXnut/vSnP6miokIvvPCCOnbs6PYX+/lo0aKFxo8fr7S0NP3Hf/yHBg8erM8//1yrV6/WJZdccsbPbtiwQaNHj9btt9+udu3aqaKiQkuXLpWvr68SExNd/bp3765169Zp9uzZioqKUkxMjHr16uVRvZ06dVJ8fLzbo+fSb9/NScOGDdOECRN066236tFHH9Uvv/yi+fPnq127dm5B6Xxqa9GihVJTUzV16lQNHDhQ//mf/6nc3FzNmzdP11xzzWl/RBKozwg7AM5bly5dlJmZqZSUFE2aNEktW7bU1KlTdejQIY/DjiQ99dRTCgwM1IIFC7Rx40b16tVLa9euPetsRdeuXRUfH6/33ntPBw4cUKNGjdS1a1etXr1avXv3dvWbPXu2Ro4cqYkTJ+rXX39VUlKSx2Hn+uuvl9Pp1NSpU5WXl6fY2FhlZGSoS5curj7NmzfXypUrlZKSoieeeEIxMTFKS0vT3r17q4Wd86ltypQpatGihV588UWNGzdOzZo108iRIzVjxgy339gB8BuHxSo1AABgMNbsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjd/ZkVRVVaWDBw+qSZMmXv05eQAAUHMsy9KxY8cUFRUlH5/Tz98QdiQdPHhQrVq1srsMAADggf3796tly5anPU7YkdSkSRNJv31ZwcHBNlcDAADORXFxsVq1auX6e/x0CDv619uGg4ODCTsAAFxkzrYEhQXKAADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRbA87Bw4c0D333KPmzZurYcOG6ty5s3bs2OE6blmWJk2apMjISDVs2FBxcXHau3ev2zmOHDmi4cOHKzg4WKGhoRoxYoRKSkpqeygAAKAOsjXs/Pzzz+rTp48aNGig1atXa8+ePXr22WfVtGlTV59Zs2Zp7ty5WrBggT799FMFBQUpPj5ex48fd/UZPny4du/eraysLK1atUqbN2/WyJEj7RgSAACoYxyWZVl2XfzJJ5/UJ598oo8++uiUxy3LUlRUlB577DGNHz9eklRUVKTw8HBlZGRo2LBh+vrrrxUbG6vt27erR48ekqQ1a9Zo8ODB+umnnxQVFXXWOoqLixUSEqKioiJeFwEAwEXiXP/+tnVm591331WPHj10++23KywsTFdddZVeeeUV1/F9+/YpPz9fcXFxrraQkBD16tVL2dnZkqTs7GyFhoa6go4kxcXFycfHR59++ukpr1tWVqbi4mK3DQAAmMnWsPP9999r/vz5atu2rTIzMzVq1Cg9+uijWrx4sSQpPz9fkhQeHu72ufDwcNex/Px8hYWFuR338/NTs2bNXH1+Ly0tTSEhIa6tVatW3h4aAACoI2wNO1VVVbr66qs1Y8YMXXXVVRo5cqQeeughLViwoEavm5qaqqKiIte2f//+Gr0eAACwj5+dF4+MjFRsbKxbW4cOHfSPf/xDkhQRESFJKigoUGRkpKtPQUGBunXr5upTWFjodo6KigodOXLE9fnfCwgIUEBAgLeGAXikzZPvn7XPDzMTaqESADCbrTM7ffr0UW5urlvbN998o+joaElSTEyMIiIitH79etfx4uJiffrpp3I6nZIkp9Opo0ePKicnx9Vnw4YNqqqqUq9evWphFAAAoC6zdWZn3Lhx+sMf/qAZM2bojjvu0LZt2/Tyyy/r5ZdfliQ5HA6NHTtWTz31lNq2bauYmBj95S9/UVRUlIYMGSLpt5mggQMHum5/nThxQqNHj9awYcPO6UksAABgNlvDzjXXXKOVK1cqNTVV06ZNU0xMjObMmaPhw4e7+jzxxBMqLS3VyJEjdfToUfXt21dr1qxRYGCgq89rr72m0aNHq3///vLx8VFiYqLmzp1rx5AAAEAdY+vv7NQV/M4O7MCaHQC4MBfF7+wAAADUNMIOAAAwGmEHAAAYjbADAACMZuvTWKibWDgLADAJMzsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0fzsLsB0bZ58/6x9fpiZUAuVAABQPzGzAwAAjEbYAQAARrM17EyZMkUOh8Nta9++vev48ePHlZycrObNm6tx48ZKTExUQUGB2zny8vKUkJCgRo0aKSwsTI8//rgqKipqeygAAKCOsn3NTseOHbVu3TrXvp/fv0oaN26c3n//fa1YsUIhISEaPXq0hg4dqk8++USSVFlZqYSEBEVERGjLli06dOiQ7rvvPjVo0EAzZsyo9bEAAIC6x/aw4+fnp4iIiGrtRUVFWrhwoZYtW6Z+/fpJkhYtWqQOHTpo69at6t27t9auXas9e/Zo3bp1Cg8PV7du3TR9+nRNmDBBU6ZMkb+/f20PBwAA1DG2r9nZu3evoqKidNlll2n48OHKy8uTJOXk5OjEiROKi4tz9W3fvr1at26t7OxsSVJ2drY6d+6s8PBwV5/4+HgVFxdr9+7dp71mWVmZiouL3TYAAGAmW8NOr169lJGRoTVr1mj+/Pnat2+frr32Wh07dkz5+fny9/dXaGio22fCw8OVn58vScrPz3cLOiePnzx2OmlpaQoJCXFtrVq18u7AAABAnWHrbaxBgwa5/rlLly7q1auXoqOj9cYbb6hhw4Y1dt3U1FSlpKS49ouLiwk8AAAYyvbbWP8uNDRU7dq107fffquIiAiVl5fr6NGjbn0KCgpca3wiIiKqPZ11cv9U64BOCggIUHBwsNsGAADMVKfCTklJib777jtFRkaqe/fuatCggdavX+86npubq7y8PDmdTkmS0+nUrl27VFhY6OqTlZWl4OBgxcbG1nr9AACg7rH1Ntb48eN18803Kzo6WgcPHtTkyZPl6+uru+66SyEhIRoxYoRSUlLUrFkzBQcH65FHHpHT6VTv3r0lSQMGDFBsbKzuvfdezZo1S/n5+Zo4caKSk5MVEBBg59AAAEAdYWvY+emnn3TXXXfpn//8p1q0aKG+fftq69atatGihSTpueeek4+PjxITE1VWVqb4+HjNmzfP9XlfX1+tWrVKo0aNktPpVFBQkJKSkjRt2jS7hgQAAOoYW8PO8uXLz3g8MDBQ6enpSk9PP22f6OhoffDBB94uDQAAGKJOrdkBAADwNsIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaHUm7MycOVMOh0Njx451tR0/flzJyclq3ry5GjdurMTERBUUFLh9Li8vTwkJCWrUqJHCwsL0+OOPq6KioparBwAAdVWdCDvbt2/XSy+9pC5duri1jxs3Tu+9955WrFihTZs26eDBgxo6dKjreGVlpRISElReXq4tW7Zo8eLFysjI0KRJk2p7CAAAoI6yPeyUlJRo+PDheuWVV9S0aVNXe1FRkRYuXKjZs2erX79+6t69uxYtWqQtW7Zo69atkqS1a9dqz549+vvf/65u3bpp0KBBmj59utLT01VeXm7XkAAAQB1ie9hJTk5WQkKC4uLi3NpzcnJ04sQJt/b27durdevWys7OliRlZ2erc+fOCg8Pd/WJj49XcXGxdu/efdprlpWVqbi42G0DAABm8rPz4suXL9dnn32m7du3VzuWn58vf39/hYaGurWHh4crPz/f1effg87J4yePnU5aWpqmTp16gdUDAICLgW0zO/v379eYMWP02muvKTAwsFavnZqaqqKiIte2f//+Wr0+AACoPbaFnZycHBUWFurqq6+Wn5+f/Pz8tGnTJs2dO1d+fn4KDw9XeXm5jh496va5goICRURESJIiIiKqPZ11cv9kn1MJCAhQcHCw2wYAAMxk222s/v37a9euXW5tDzzwgNq3b68JEyaoVatWatCggdavX6/ExERJUm5urvLy8uR0OiVJTqdTTz/9tAoLCxUWFiZJysrKUnBwsGJjY2t3QLBVmyffP2ufH2Ym1EIlAIC6xraw06RJE3Xq1MmtLSgoSM2bN3e1jxgxQikpKWrWrJmCg4P1yCOPyOl0qnfv3pKkAQMGKDY2Vvfee69mzZql/Px8TZw4UcnJyQoICKj1MQEAgLrH1gXKZ/Pcc8/Jx8dHiYmJKisrU3x8vObNm+c67uvrq1WrVmnUqFFyOp0KCgpSUlKSpk2bZmPVAACgLqlTYefDDz902w8MDFR6errS09NP+5no6Gh98MEHNVwZAAC4WNn+OzsAAAA1qU7N7ADA+WJxOoCzYWYHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIzm0esivv/+e1122WXergUwxrm8wgAAUDs8CjtXXHGFrr/+eo0YMUK33XabAgMDvV0XAC/jHVIA6iuPbmN99tln6tKli1JSUhQREaH/+q//0rZt27xdGwAAwAXzKOx069ZNzz//vA4ePKhXX31Vhw4dUt++fdWpUyfNnj1bhw8f9nadAAAAHrmgBcp+fn4aOnSoVqxYob/+9a/69ttvNX78eLVq1Ur33XefDh065K06AQAAPHJBYWfHjh3605/+pMjISM2ePVvjx4/Xd999p6ysLB08eFC33HKLt+oEAADwiEcLlGfPnq1FixYpNzdXgwcP1pIlSzR48GD5+PyWnWJiYpSRkaE2bdp4s1YAAIDz5lHYmT9/vh588EHdf//9ioyMPGWfsLAwLVy48IKKAyQe4wYAXBiPws7evXvP2sff319JSUmenB4AAMBrPFqzs2jRIq1YsaJa+4oVK7R48eILLgoAAMBbPAo7aWlpuuSSS6q1h4WFacaMGRdcFAAAgLd4FHby8vIUExNTrT06Olp5eXkXXBQAAIC3eBR2wsLC9OWXX1Zr/+KLL9S8efMLLgoAAMBbPAo7d911lx599FFt3LhRlZWVqqys1IYNGzRmzBgNGzbM2zUCAAB4zKOnsaZPn64ffvhB/fv3l5/fb6eoqqrSfffdx5odAABQp3gUdvz9/fX6669r+vTp+uKLL9SwYUN17txZ0dHR3q4PAADggngUdk5q166d2rVr561aAHiIH14EgNPzKOxUVlYqIyND69evV2FhoaqqqtyOb9iwwSvFAQAAXCiPws6YMWOUkZGhhIQEderUSQ6Hw9t1AQAAeIVHYWf58uV64403NHjwYG/XAwAA4FUePXru7++vK664wtu1AAAAeJ1HYeexxx7T888/L8uyvF0PAACAV3l0G+vjjz/Wxo0btXr1anXs2FENGjRwO/7WW295pTgAAIAL5VHYCQ0N1a233urtWgAAALzOo7CzaNEib9cBAABQIzxasyNJFRUVWrdunV566SUdO3ZMknTw4EGVlJR4rTgAAIAL5dHMzo8//qiBAwcqLy9PZWVluummm9SkSRP99a9/VVlZmRYsWODtOgEAADzi0czOmDFj1KNHD/38889q2LChq/3WW2/V+vXrvVYcAADAhfJoZuejjz7Sli1b5O/v79bepk0bHThwwCuFoX7gnU4AgJrm0cxOVVWVKisrq7X/9NNPatKkyQUXBQAA4C0ehZ0BAwZozpw5rn2Hw6GSkhJNnjyZV0gAAIA6xaPbWM8++6zi4+MVGxur48eP6+6779bevXt1ySWX6H//93+9XSMAAIDHPAo7LVu21BdffKHly5fryy+/VElJiUaMGKHhw4e7LVgGAACwm0dhR5L8/Px0zz33eLMWAAAAr/Mo7CxZsuSMx++77z6PigEAAPA2j8LOmDFj3PZPnDihX375Rf7+/mrUqBFhBwAA1BkePY31888/u20lJSXKzc1V3759z2uB8vz589WlSxcFBwcrODhYTqdTq1evdh0/fvy4kpOT1bx5czVu3FiJiYkqKChwO0deXp4SEhLUqFEjhYWF6fHHH1dFRYUnwwIAAAby+N1Yv9e2bVvNnDmz2qzPmbRs2VIzZ85UTk6OduzYoX79+umWW27R7t27JUnjxo3Te++9pxUrVmjTpk06ePCghg4d6vp8ZWWlEhISVF5eri1btmjx4sXKyMjQpEmTvDUsAABwkfN4gfIpT+bnp4MHD55z/5tvvtlt/+mnn9b8+fO1detWtWzZUgsXLtSyZcvUr18/Sb+9bb1Dhw7aunWrevfurbVr12rPnj1at26dwsPD1a1bN02fPl0TJkzQlClTqv3CM3A25/KLzj/MTKiFSgAA3uJR2Hn33Xfd9i3L0qFDh/Tiiy+qT58+HhVSWVmpFStWqLS0VE6nUzk5OTpx4oTi4uJcfdq3b6/WrVsrOztbvXv3VnZ2tjp37qzw8HBXn/j4eI0aNUq7d+/WVVdddcprlZWVqayszLVfXFzsUc0AAKDu8yjsDBkyxG3f4XCoRYsW6tevn5599tnzOteuXbvkdDp1/PhxNW7cWCtXrlRsbKx27twpf39/hYaGuvUPDw9Xfn6+JCk/P98t6Jw8fvLY6aSlpWnq1KnnVSdgB94dBgAXzqOwU1VV5bUCrrzySu3cuVNFRUV68803lZSUpE2bNnnt/KeSmpqqlJQU135xcbFatWpVo9cEAAD28OqaHU/4+/vriiuukCR1795d27dv1/PPP68777xT5eXlOnr0qNvsTkFBgSIiIiRJERER2rZtm9v5Tj6tdbLPqQQEBCggIMDLIwEAAHWRR2Hn32dFzmb27Nnnde6qqiqVlZWpe/fuatCggdavX6/ExERJUm5urvLy8uR0OiVJTqdTTz/9tAoLCxUWFiZJysrKUnBwsGJjY8/rugAAwEwehZ3PP/9cn3/+uU6cOKErr7xSkvTNN9/I19dXV199taufw+E443lSU1M1aNAgtW7dWseOHdOyZcv04YcfKjMzUyEhIRoxYoRSUlLUrFkzBQcH65FHHpHT6VTv3r0l/fb29djYWN17772aNWuW8vPzNXHiRCUnJzNzAwAAJHkYdm6++WY1adJEixcvVtOmTSX99kODDzzwgK699lo99thj53SewsJC3XfffTp06JBCQkLUpUsXZWZm6qabbpIkPffcc/Lx8VFiYqLKysoUHx+vefPmuT7v6+urVatWadSoUXI6nQoKClJSUpKmTZvmybAAAICBPAo7zz77rNauXesKOpLUtGlTPfXUUxowYMA5h52FCxee8XhgYKDS09OVnp5+2j7R0dH64IMPzq1wAABQ73j0C8rFxcU6fPhwtfbDhw/r2LFjF1wUAACAt3gUdm699VY98MADeuutt/TTTz/pp59+0j/+8Q+NGDHC7XUOAAAAdvPoNtaCBQs0fvx43X333Tpx4sRvJ/Lz04gRI/TMM894tUAAAIAL4VHYadSokebNm6dnnnlG3333nSTp8ssvV1BQkFeLAwAAuFAX9NbzQ4cO6dChQ2rbtq2CgoJkWZa36gIAAPAKj8LOP//5T/Xv31/t2rXT4MGDdejQIUnSiBEjzvlJLAAAgNrgUdgZN26cGjRooLy8PDVq1MjVfuedd2rNmjVeKw4AAOBCebRmZ+3atcrMzFTLli3d2tu2basff/zRK4UBAAB4g0dhp7S01G1G56QjR47wmgYYr82T79tdAgDgPHh0G+vaa6/VkiVLXPsOh0NVVVWaNWuWbrzxRq8VBwAAcKE8mtmZNWuW+vfvrx07dqi8vFxPPPGEdu/erSNHjuiTTz7xdo0AAAAe8yjsdOrUSd98841efPFFNWnSRCUlJRo6dKiSk5MVGRnp7RoBXGTO5VbfDzMTaqESAPAg7Jw4cUIDBw7UggUL9N///d81URMAAIDXnPeanQYNGujLL7+siVoAAAC8zqMFyvfcc48WLlzo7VoAAAC8zqM1OxUVFXr11Ve1bt06de/evdo7sWbPnu2V4gAAAC7UeYWd77//Xm3atNFXX32lq6++WpL0zTffuPVxOBzeqw4AAOACnVfYadu2rQ4dOqSNGzdK+u31EHPnzlV4eHiNFAcAAHChzmvNzu/far569WqVlpZ6tSAAAABv8miB8km/Dz8AAAB1zXmFHYfDUW1NDmt0AABAXXZea3Ysy9L999/vetnn8ePH9fDDD1d7Guutt97yXoUAAAAX4LzCTlJSktv+Pffc49ViAAAAvO28ws6iRYtqqg4AAIAacUELlAEAAOo6wg4AADCaR6+LAIDa0ObJ9+0uAYABmNkBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGi8CBSAy7m8ePOHmQm1UAkAeA8zOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0WwNO2lpabrmmmvUpEkThYWFaciQIcrNzXXrc/z4cSUnJ6t58+Zq3LixEhMTVVBQ4NYnLy9PCQkJatSokcLCwvT444+roqKiNocCAADqKFvDzqZNm5ScnKytW7cqKytLJ06c0IABA1RaWurqM27cOL333ntasWKFNm3apIMHD2ro0KGu45WVlUpISFB5ebm2bNmixYsXKyMjQ5MmTbJjSAAAoI6x9XURa9ascdvPyMhQWFiYcnJydN1116moqEgLFy7UsmXL1K9fP0nSokWL1KFDB23dulW9e/fW2rVrtWfPHq1bt07h4eHq1q2bpk+frgkTJmjKlCny9/e3Y2gAAKCOqFNrdoqKiiRJzZo1kyTl5OToxIkTiouLc/Vp3769WrdurezsbElSdna2OnfurPDwcFef+Ph4FRcXa/fu3ae8TllZmYqLi902AABgpjoTdqqqqjR27Fj16dNHnTp1kiTl5+fL399foaGhbn3Dw8OVn5/v6vPvQefk8ZPHTiUtLU0hISGurVWrVl4eDQAAqCvqTNhJTk7WV199peXLl9f4tVJTU1VUVOTa9u/fX+PXBAAA9rB1zc5Jo0eP1qpVq7R582a1bNnS1R4REaHy8nIdPXrUbXanoKBAERERrj7btm1zO9/Jp7VO9vm9gIAABQQEeHkUAACgLrJ1ZseyLI0ePVorV67Uhg0bFBMT43a8e/fuatCggdavX+9qy83NVV5enpxOpyTJ6XRq165dKiwsdPXJyspScHCwYmNja2cgAACgzrJ1Zic5OVnLli3TO++8oyZNmrjW2ISEhKhhw4YKCQnRiBEjlJKSombNmik4OFiPPPKInE6nevfuLUkaMGCAYmNjde+992rWrFnKz8/XxIkTlZyczOwNAACwN+zMnz9fknTDDTe4tS9atEj333+/JOm5556Tj4+PEhMTVVZWpvj4eM2bN8/V19fXV6tWrdKoUaPkdDoVFBSkpKQkTZs2rbaGAQAA6jBbw45lWWftExgYqPT0dKWnp5+2T3R0tD744ANvlgYAAAxRZ57GAgAAqAmEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARqsTr4sAcPFo8+T7deo8AHA2zOwAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDReBIp6gxdPAkD9xMwOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADAaYQcAABiNsAMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAoxF2AACA0Qg7AADAaIQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGszXsbN68WTfffLOioqLkcDj09ttvux23LEuTJk1SZGSkGjZsqLi4OO3du9etz5EjRzR8+HAFBwcrNDRUI0aMUElJSS2OAgAA1GW2hp3S0lJ17dpV6enppzw+a9YszZ07VwsWLNCnn36qoKAgxcfH6/jx464+w4cP1+7du5WVlaVVq1Zp8+bNGjlyZG0NAQAA1HF+dl580KBBGjRo0CmPWZalOXPmaOLEibrlllskSUuWLFF4eLjefvttDRs2TF9//bXWrFmj7du3q0ePHpKkF154QYMHD9b//M//KCoqqtbGAgAA6iZbw86Z7Nu3T/n5+YqLi3O1hYSEqFevXsrOztawYcOUnZ2t0NBQV9CRpLi4OPn4+OjTTz/Vrbfeespzl5WVqayszLVfXFxccwOpx9o8+b7dJQAAUHcXKOfn50uSwsPD3drDw8Ndx/Lz8xUWFuZ23M/PT82aNXP1OZW0tDSFhIS4tlatWnm5egAAUFfU2bBTk1JTU1VUVOTa9u/fb3dJAACghtTZsBMRESFJKigocGsvKChwHYuIiFBhYaHb8YqKCh05csTV51QCAgIUHBzstgEAADPV2bATExOjiIgIrV+/3tVWXFysTz/9VE6nU5LkdDp19OhR5eTkuPps2LBBVVVV6tWrV63XDAAA6h5bFyiXlJTo22+/de3v27dPO3fuVLNmzdS6dWuNHTtWTz31lNq2bauYmBj95S9/UVRUlIYMGSJJ6tChgwYOHKiHHnpICxYs0IkTJzR69GgNGzaMJ7EAAIAkm8POjh07dOONN7r2U1JSJElJSUnKyMjQE088odLSUo0cOVJHjx5V3759tWbNGgUGBro+89prr2n06NHq37+/fHx8lJiYqLlz59b6WAAAQN1ka9i54YYbZFnWaY87HA5NmzZN06ZNO22fZs2aadmyZTVRHgAAMECdXbMDAADgDYQdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBohB0AAGA0wg4AADCare/GwsWrzZPv210CAADnhJkdAABgNMIOAAAwGmEHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDR+J2dOuBcfrPmh5kJtVAJAADmIezUM/wYIACgvuE2FgAAMBphBwAAGI2wAwAAjMaaHYOwHgcAgOqY2QEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDSexrpI8KQVAACeYWYHAAAYjbADAACMRtgBAABGI+wAAACjEXYAAIDRCDsAAMBoPHoOABepc/lJih9mJtRCJUDdRtgBgDqI39YCvIfbWAAAwGiEHQAAYDTCDgAAMBphBwAAGI0FygCMx1NLQP3GzA4AADAaMzsA4EXMIp0Z3w/swMwOAAAwGjM7AHCO+KE/4OJkTNhJT0/XM888o/z8fHXt2lUvvPCCevbsaXdZAFDncWsJpjMi7Lz++utKSUnRggUL1KtXL82ZM0fx8fHKzc1VWFiY3eUBQL3AzBfqKiPCzuzZs/XQQw/pgQcekCQtWLBA77//vl599VU9+eSTNlcHALiYeWvmixk0+1z0Yae8vFw5OTlKTU11tfn4+CguLk7Z2dmn/ExZWZnKyspc+0VFRZKk4uJir9dXVfaL188JwPvO5c+/t/48tx63wivnORfeGlddq7k2ncv3463vuTbH3mly5ln7fDU1vtbO44mT35dlWWfuaF3kDhw4YEmytmzZ4tb++OOPWz179jzlZyZPnmxJYmNjY2NjYzNg279//xmzwkU/s+OJ1NRUpaSkuParqqp05MgRNW/eXA6Hw2vXKS4uVqtWrbR//34FBwd77bx1WX0cs1Q/x10fxyzVz3HXxzFL9XPcF9uYLcvSsWPHFBUVdcZ+F33YueSSS+Tr66uCggK39oKCAkVERJzyMwEBAQoICHBrCw0NrakSFRwcfFH8j8ab6uOYpfo57vo4Zql+jrs+jlmqn+O+mMYcEhJy1j4X/Y8K+vv7q3v37lq/fr2rraqqSuvXr5fT6bSxMgAAUBdc9DM7kpSSkqKkpCT16NFDPXv21Jw5c1RaWup6OgsAANRfRoSdO++8U4cPH9akSZOUn5+vbt26ac2aNQoPD7e1roCAAE2ePLnaLTOT1ccxS/Vz3PVxzFL9HHd9HLNUP8dt6pgdlnW257UAAAAuXhf9mh0AAIAzIewAAACjEXYAAIDRCDsAAMBohB0AAGA0wk4NSU9PV5s2bRQYGKhevXpp27ZtdpdUozZv3qybb75ZUVFRcjgcevvtt+0uqcalpaXpmmuuUZMmTRQWFqYhQ4YoNzfX7rJq3Pz589WlSxfXL6w6nU6tXr3a7rJq1cyZM+VwODR27Fi7S6lRU6ZMkcPhcNvat29vd1k17sCBA7rnnnvUvHlzNWzYUJ07d9aOHTvsLqtGtWnTptq/a4fDoeTkZLtL8wrCTg14/fXXlZKSosmTJ+uzzz5T165dFR8fr8LCQrtLqzGlpaXq2rWr0tPT7S6l1mzatEnJycnaunWrsrKydOLECQ0YMEClpaV2l1ajWrZsqZkzZyonJ0c7duxQv379dMstt2j37t12l1Yrtm/frpdeekldunSxu5Ra0bFjRx06dMi1ffzxx3aXVKN+/vln9enTRw0aNNDq1au1Z88ePfvss2ratKndpdWo7du3u/17zsrKkiTdfvvtNlfmJd559zj+Xc+ePa3k5GTXfmVlpRUVFWWlpaXZWFXtkWStXLnS7jJqXWFhoSXJ2rRpk92l1LqmTZtaf/vb3+wuo8YdO3bMatu2rZWVlWVdf/311pgxY+wuqUZNnjzZ6tq1q91l1KoJEyZYffv2tbsM240ZM8a6/PLLraqqKrtL8QpmdrysvLxcOTk5iouLc7X5+PgoLi5O2dnZNlaGmlZUVCRJatasmc2V1J7KykotX75cpaWl9eJddMnJyUpISHD78226vXv3KioqSpdddpmGDx+uvLw8u0uqUe+++6569Oih22+/XWFhYbrqqqv0yiuv2F1WrSovL9ff//53Pfjgg3I4HHaX4xWEHS/7v//7P1VWVlZ7VUV4eLjy8/Ntqgo1raqqSmPHjlWfPn3UqVMnu8upcbt27VLjxo0VEBCghx9+WCtXrlRsbKzdZdWo5cuX67PPPlNaWprdpdSaXr16KSMjQ2vWrNH8+fO1b98+XXvttTp27JjdpdWY77//XvPnz1fbtm2VmZmpUaNG6dFHH9XixYvtLq3WvP322zp69Kjuv/9+u0vxGiPejQXYLTk5WV999ZXx6xlOuvLKK7Vz504VFRXpzTffVFJSkjZt2mRs4Nm/f7/GjBmjrKwsBQYG2l1OrRk0aJDrn7t06aJevXopOjpab7zxhkaMGGFjZTWnqqpKPXr00IwZMyRJV111lb766istWLBASUlJNldXOxYuXKhBgwYpKirK7lK8hpkdL7vkkkvk6+urgoICt/aCggJFRETYVBVq0ujRo7Vq1Spt3LhRLVu2tLucWuHv768rrrhC3bt3V1pamrp27arnn3/e7rJqTE5OjgoLC3X11VfLz89Pfn5+2rRpk+bOnSs/Pz9VVlbaXWKtCA0NVbt27fTtt9/aXUqNiYyMrBbaO3ToYPztu5N+/PFHrVu3Tn/84x/tLsWrCDte5u/vr+7du2v9+vWutqqqKq1fv75erGmoTyzL0ujRo7Vy5Upt2LBBMTExdpdkm6qqKpWVldldRo3p37+/du3apZ07d7q2Hj16aPjw4dq5c6d8fX3tLrFWlJSU6LvvvlNkZKTdpdSYPn36VPsJiW+++UbR0dE2VVS7Fi1apLCwMCUkJNhdildxG6sGpKSkKCkpST169FDPnj01Z84clZaW6oEHHrC7tBpTUlLi9v/29u3bp507d6pZs2Zq3bq1jZXVnOTkZC1btkzvvPOOmjRp4lqTFRISooYNG9pcXc1JTU3VoEGD1Lp1ax07dkzLli3Thx9+qMzMTLtLqzFNmjSpthYrKChIzZs3N3qN1vjx43XzzTcrOjpaBw8e1OTJk+Xr66u77rrL7tJqzLhx4/SHP/xBM2bM0B133KFt27bp5Zdf1ssvv2x3aTWuqqpKixYtUlJSkvz8DIsHdj8OZqoXXnjBat26teXv72/17NnT2rp1q90l1aiNGzdakqptSUlJdpdWY041XknWokWL7C6tRj344INWdHS05e/vb7Vo0cLq37+/tXbtWrvLqnX14dHzO++804qMjLT8/f2tSy+91Lrzzjutb7/91u6yatx7771nderUyQoICLDat29vvfzyy3aXVCsyMzMtSVZubq7dpXidw7Isy56YBQAAUPNYswMAAIxG2AEAAEYj7AAAAKMRdgAAgNEIOwAAwGiEHQAAYDTCDgAAMBphBwAAGI2wAwAAjEbYAQAARiPsAAAAo/0/GJiD9GLVCicAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 32%|███████████████████████▉                                                    | 974/3092 [1:07:54<2:21:41,  4.01s/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_private_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_private_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
}