Random-Sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random Sampling\n",
    "\n",
    "In this notebook you'll work with the obstacle's polygon representation itself.\n",
    "\n",
    "Your tasks will be:\n",
    "\n",
    "1. Create polygons.\n",
    "2. Sample random 3D points.\n",
    "3. Remove points contained by an obstacle polygon.\n",
    "\n",
    "Recall, a point $(x, y, z)$ collides with a polygon if the $(x, y)$ coordinates are contained by the polygon and the $z$ coordinate (height) is less than the height of the polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from shapely.geometry import Polygon, Point\n",
    "%matplotlib inline \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = 12, 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-310.2389   -439.2315     85.5         5.          5.         85.5     ]\n",
      " [-300.2389   -439.2315     85.5         5.          5.         85.5     ]\n",
      " [-290.2389   -439.2315     85.5         5.          5.         85.5     ]\n",
      " ...\n",
      " [ 257.8061    425.1645      1.75852     1.292725    1.292725    1.944791]\n",
      " [ 293.9967    368.3391      3.557666    1.129456    1.129456    3.667319]\n",
      " [ 281.5162    354.4156      4.999351    1.053772    1.053772    4.950246]]\n"
     ]
    }
   ],
   "source": [
    "# This is the same obstacle data from the previous lesson.\n",
    "filename = 'colliders.csv'\n",
    "data = np.loadtxt(filename, delimiter=',', dtype='Float64', skiprows=2)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.003809022184238454\n"
     ]
    }
   ],
   "source": [
    "test_start = [-122.397450 , 37.792480]\n",
    "test_goal = [ -122.396525 ,37.796175 ]\n",
    "dist = np.linalg.norm(np.array(test_goal) - np.array(test_start))\n",
    "print(dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_polygons(data):\n",
    "\n",
    "    polygons = []\n",
    "    for i in range(data.shape[0]):\n",
    "        north, east, alt, d_north, d_east, d_alt = data[i, :]\n",
    "        obstacle = [north-d_north, east-d_east, north+d_north, east+d_east]\n",
    "        # TODO: Extract the 4 corners of the obstacle\n",
    "        # \n",
    "        # NOTE: The order of the points matters since\n",
    "        # `shapely` draws the sequentially from point to point.\n",
    "        #\n",
    "        # If the area of the polygon is 0 you've likely got a weird\n",
    "        # order.\n",
    "        corners = [(obstacle[0], obstacle[1]), (obstacle[0], obstacle[3]),\n",
    "                   (obstacle[2], obstacle[3]), (obstacle[2], obstacle[1])]\n",
    "        \n",
    "        # TODO: Compute the height of the polygon\n",
    "        height = alt + d_alt\n",
    "\n",
    "        # TODO: Once you've defined corners, define polygons\n",
    "        #p = Polygon(corners)\n",
    "        p = Polygon(corners)\n",
    "        #polygons.append((p, height))\n",
    "        polygons.append((p, height))\n",
    "    return polygons\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "polygons = extract_polygons(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling 3D Points\n",
    "\n",
    "Now that we have the extracted the polygons, we need to sample random 3D points. Currently we don't know suitable ranges for x, y, and z. Let's figure out the max and min values for each dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X\n",
      "min = -315.2389, max = 604.7611\n",
      "\n",
      "Y\n",
      "min = -444.2315, max = 475.7685\n",
      "\n",
      "Z\n",
      "min = 0, max = 10\n"
     ]
    }
   ],
   "source": [
    "xmin = np.min(data[:, 0] - data[:, 3])\n",
    "xmax = np.max(data[:, 0] + data[:, 3])\n",
    "\n",
    "ymin = np.min(data[:, 1] - data[:, 4])\n",
    "ymax = np.max(data[:, 1] + data[:, 4])\n",
    "\n",
    "zmin = 0\n",
    "# Limit the z axis for the visualization\n",
    "zmax = 10\n",
    "\n",
    "print(\"X\")\n",
    "print(\"min = {0}, max = {1}\\n\".format(xmin, xmax))\n",
    "\n",
    "print(\"Y\")\n",
    "print(\"min = {0}, max = {1}\\n\".format(ymin, ymax))\n",
    "\n",
    "print(\"Z\")\n",
    "print(\"min = {0}, max = {1}\".format(zmin, zmax))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, it's time to sample points. All that's left is picking the distribution and number of samples. The uniform distribution makes sense in this situation since we we'd like to encourage searching the whole space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = 100\n",
    "\n",
    "xvals = np.random.uniform(xmin, xmax, num_samples)\n",
    "yvals = np.random.uniform(ymin, ymax, num_samples)\n",
    "zvals = np.random.uniform(zmin, zmax, num_samples)\n",
    "\n",
    "samples = list(zip(xvals, yvals, zvals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(278.08914319123738, 100.72655493328205, 9.5652225168559113),\n",
       " (-101.11321172176099, 250.26857863350403, 7.667133185649524),\n",
       " (-269.73873028190769, -252.62700052121644, 3.8193894885327695),\n",
       " (439.77525704105807, 67.052049616493264, 4.2507356570503232),\n",
       " (48.184308383967505, -387.93290576395611, 0.96827150319743427),\n",
       " (538.76766703383805, 305.25370048687239, 3.8731951018025024),\n",
       " (-83.363288113443815, -153.73576233688698, 6.7323118789843486),\n",
       " (-42.429755362920162, -401.5611019254917, 8.2977369663217502),\n",
       " (413.42553545056836, 165.93303736801215, 1.9117243004478701),\n",
       " (-152.85830067030213, -169.17076862528364, 0.61634326721839772)]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samples[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Removing Points Colliding With Obstacles\n",
    "\n",
    "Prior to remove a point we must determine whether it collides with any obstacle. Complete the `collides` function below. It should return `True` if the point collides with *any* obstacle and `False` if no collision is detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def collides(polygons, point):   \n",
    "    # TODO: Determine whether the point collides\n",
    "    # with any obstacles.\n",
    "    #print(point)\n",
    "    for (p, height) in polygons:\n",
    "        if p.contains(Point(point)) and height >= point[2]:\n",
    "            return True\n",
    "        \n",
    "    return False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `collides` for all points in the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time taken {0} seconds ... 6.301304340362549\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "to_keep = []\n",
    "for point in samples:\n",
    "    if not collides(polygons, point):\n",
    "        to_keep.append(point)\n",
    "time_taken = time.time() - t0\n",
    "print(\"Time taken {0} seconds ...\", time_taken)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64\n"
     ]
    }
   ],
   "source": [
    "print(len(to_keep))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Points Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from grid import create_grid\n",
    "grid = create_grid(data, zmax, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAK9CAYAAADScyEHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3W2QbNtZH/b/OlcImGswb4LAlWZGBBU2fkVMsGxSDkH+ANixqATKum5ihWBNXCYxtnHZsiYVJ2UfJZQxsqlgnBGCCGhf7AhcYEIcu3hxJR8s+wgIBmQXKunM0TUyXIyQjca2kO7Kh91zz5w5s2e6Z3r3funfr6qrT6/pM7N677V3P/30s9YutdYAAACPu9N3BwAAYKgEywAA0EKwDAAALQTLAADQQrAMAAAtBMsAANBCsAwAAC0EywAA0EKwDAAALV7Udwdu49M+7dPq/v5+391gjd75znd28nu/4Au+YGN/a9V+3MSY+35Tt33NY972Y9b3uGH6hnJuYFzu37+fX/7lXy7LPLeM+XLXBwcH9d69e313gzUqZalxu7LLxnlXf2vVftzEmPt+U7d9zWPe9mPW97hh+oZybmBcDg4Ocu/evaUGjzIMAABoIVgGAIAWgmUAAGghWAY2Yz5P9veTO3ea+/m87x4BwLVGvRoGq1tmIoTJDqzdfJ4cHianp83jk5PmcZLMZv31C9h6be+L3gs5I7MMdO/o6GGgfOb0tGkHgAGTWYYrDDWzcFW/Brm02YMHq7VfYqj7AvuGfo3ufMjoyCwD3dvdXa0dAAZCsAx07+7dZGfn0badnaYdAAZMsAx0bzZLjo+Tvb2klOb++NjkPgAGT80ysBmzmeAYgNERLK/ZbSYTXDZJweSE8bAsHwBMjzIMAABoIVimGyO4Wlut9drbGC3zum5yoz9d7VNjg6kzrlkHZRisn6u1AQATIbPM+rlaGwAwEYJl1m8NV2sDABgCwTLr52ptANMygnko0BXBMo8ppVx6W5qrtY1S235f5gZM2Nk8lJOTpNaH81AEzGwJwTLr52ptANNhHgpbzmoYdMPV2pgoS0qxdcxDYcvJLAMA7cxDYcsJlgGAduahsOUEywBAO/NQ2HJqloHB6GNlDTXIsATzUNhiMssAANBCsLxmtdYb39b1+wCAJbjYCktQhgEAbJ+zi62crSF9drGVRMkJj5BZBoAt4eqb57jYCksSLAMA28fFVoZhBKUwgmUAYPu42Er/zkphTk6SWh+WwgwsYBYsD9zFr8yWufE426tbNxmntj3QKxdbWUqn5/GRlMIIlgFgS1g96RwXW+nfSEphrIYxYZedCGXzaHOTN07jCRg1F1vp1+5uU3pxWfuAyCwDwBr0UW6ktIlRG0kpjGAZAFjdCFYxYOBGUgqjDAMAWI0LerAuIyiFkVkGAFYzklUMYB1kltfsupqxTc4+Vr/GKoY0doHlbOK4vOzc8NG0ZNsGtooByxEvXE1mGQBYSWtIPLBVDGAdZJYHri1rMNVPgVdlSYb+mrcx87qNrxmWMfVj441J3pLkyfONA1zFANZBZhkA1mCbLvbxTJLXJ7mfDHoVA1gHmWUAYGXPLG71+ef77gp0SrAMADdwWWnYMldO3YbMM0yJMgwAgItcdIUFmWW4xDKTCVfJDt309w1pObebTrCURQNGZ6QXXeljIvxVf3Mq53+ZZQCA81x0hXNkloHB6HLpwKlkOBiOZcfU2Mfe2Pt/I20XV3HRla0kswwAN1BKeezGRLRdXMVFV7aSYBkA4Ly7d5uLrJznoitbS7AMAHDebNZcZGVvz0VXULMMN7Xur1x9hQvbwbrLIzGbDTY4Hsv7xVRWyhAsA7BVzr+Bj+kNG+iHYHnCbrJuL40uV2W4iU28oQ89aBh6/yBpzg/GKkyLmmWATXNlMIDRkFkG2KSRXhmM9ZF5hnGRWQbYJFcG612t9YUbwHVklgE2qeMrg1lpAWC9ZJYBNsmVwbbOZVf6c+U/tt2YjgfBMsAmuTLYZJwv51DaAdOlDGOknJCHy77hSmeT+I6OmtKL3d0mUDa5D2CQBMsAm9bhlcF8WANYL2UYAADQQrAMsGkuSgIwGsowADbJRUkARkVmeaSWWYpo3cuxjGmZFxisCV2UxJJowDYQLANsUscXJQFgvQTLsKK2tVWtscpSlrwoybVZWnXPo3HdOcN5A4ZNsAywSeu4KMlZ3fPJSVLrw7pnATPA2gmWATZpNkuOj5O9vaSU5v74eLXJfROqewYYOqthAGzabS9Kou6ZAbhsAqeSEqZIZnnD+ljFAhifK+tal6x7HgLnL2BVQ4uFBMsAY7OOumcAliJYBhibddQ9j9QQskzAdlGzzGiohYNzblv3vAZtx6RAFrrl2NsswTIATNTF4EnSAVbXaRlGKeVPlVJ+tpTyM6WUZ0opH1dKeXkp5R2llJ8vpfytUsqLF8/92MXjdy9+vt9l3wDohgtuTMN1k6xcXIVt0VmwXEp5KsmfSHJQa/2tSZ5I8tok35jkzbXWVyT5QJKvXfyXr03ygVrr5yR58+J5APACgRmwaV1P8HtRko8vpbwoyU6S9yf5kiRvX/z8bUm+YvHv1yweZ/HzVxfFN2y5my41eJsbDJ1xu13sb8709d7VWbBca/0XSb4pyYM0QfIHk7wzya/WWj+yeNqzSZ5a/PupJO9b/N+PLJ7/qV31DwCmTpkE3F5nE/xKKZ+cJlv88iS/muT/SPJllzz17Oi97GPBY0d2KeUwyWGS7A5wAf5NcdID4DZM/pueq/ahjPzNdVmG8fuSvLfW+lyt9deTfH+S35PkkxZlGUny0iS/sPj3s0leliSLn//GJL9y8ZfWWo9rrQe11oOXvOQlHXYfYDW+KgaYni6D5QdJXlVK2VnUHr86yc8l+bEkX7l4zuuS/MDi3z+4eJzFz3+0+pgLMGk+XKxmk9vrstUuvC2zjTorw6i1vqOU8vYkP5HkI0l+Mslxkv8zyfeWUv7Sou2ti//y1iTfXUp5d5qM8mu76hvA0LUFQ4IVgM3q9KIktda/kOQvXGh+T5IvvOS5/y7JV3XZHwD6d9us6KY+MEy9prfr17Ou7Te17c74uIIf8AgZzZuzjQCmR7AMACPhAxlsnmB5zZzIAKbB+ZwpsazczXV9BT8AJuixVRnm82R/P7lzp7mfz/vsHgNgFQ2mQmYZYIBGFVzM58nhYXJ62jw+OWkeJ8ls9tjTR/XagK0ns9w12RZg6o6OHgbKZ05Pm3Y6dfFCONashvWTWe7SitkWgLF4JDt8pyXv8uDBZjoD0CGZ5S4NMNvSloW47Q3YYru7S7U7f/Rrk9vd+wab1uVYEyx3qS2rItsCTMndu8nOzqNtOztNO8DICZa7tGS2BWDUZrPk+DjZ20tKae6Pj5WbAZMgWL6N6ybvybYA22I2S+7fT55/vrkXKAMTYYLfTS0zee/s/uioKb3Y3W0CZW8iwIT0cYn0i39zW5ejW+V1b+s2gtsqYz54Dg4O6r179/r54/v7TYB80d5ek1UZqK4mVYx5HA3ZkCbB2Me0GUuwfFk/jetuLHPusu2HY0jvNV24bKwdHBzk3r17S71wmeWbapmk9/zJSZ644aBb5cQx9YENUyYrCrA5t42Z1CzfVMskPetcAIPhokgAtyZYvqmWyXtv7Kc3AI86m1dxcpLU+nBexUQC5lrrI7d1sR7welzcP5fdGA776WqC5ZtqWSrpmb77BQxOL8HXAC+KBDBGguXbsFQSrEZZwAs6z9xs8KJIMlLAlAmWgc2YeFnA4LgoUpJcmdEfbEDvQyUMimAZ2AxlAZvlokjj5EMlDI51ltfsNvWIU106bsxjrG9D2s+33o937jRv/heV0pQyTVjXWc3Wpejm862/KNLolukb6Rr+TNeQ3ofWrda61IuTWQY2Y9vKAs59lV739lK/53ua295eaimb+XrdvIrx2WCtObAcFyUBNuPu3UcvEZ9Mtyzg7Kv0s9d6cpJ8zdc0WfQPf/hh2+Fh829BbGcGn0m+aHf38szyVD9UruB8hnN0+3XE2rb1lDPOF8ksA5vRstziJAPFy+qzf/3XHwbKZ9Rsc9FIas2tR802UbO8ZmqWHzfmMda3Ie1n+3EFbfXZl9mCmu0zbePZ2LpgBLXmm1hd5GImWWZ5WIb0/nRTapYB+rLKV+a+Xr/cNi+fptYcBkWwDLBul32V/jEfk7z4xY+2DfDr9UGwfBqXcMEb+iJYBli3y+qzv/M7k+/4ju2o2b6tljW573/1V6uPBTZOzfKaqVl+3JjHWN+GtJ/tR25r6Zrllprv55M80fZ/gI0a0vvTTalZhgk4/7WjryDZGi113FYaBvogWAZgI5b+8NeyfNr+93yPD4zAxgmWARiWbVqTm95YJ5plqVleMzXLjxvzGBuqPva//bhdNrGOLvTp4hg3vlczpjikzbI1yy53DQCMkoCXTRAsA0yYq571zzdBMG6CZQBg6/hAwbIEyxPmRADclPMH11ECwbYQLAMMmIAE2jke2ATBMsCECSZYlymsfgA3IVgGADq39OXOYWAEy1zrttkEJ0K4OcfP5smgLmcIY1MAzia4gh8AALSQWV6zTX2a9akZgE266n3HRNTts8p4GDuZZVjGfJ7s7yd37jT383nfPQIANkBmGa4znyeHh8npafP45KR5nCSzWX/9AhgR2WbGSmYZrnN09DBQPnN62rQDAJMmswzXefBgtXaALdNX1li2mk0QLA/A1ArhJ2d3tym9uKwdWNl157xNBEDOu8CylGHAde7eTXZ2Hm3b2WnaAYBJEyzDdWaz5Pg42dtLSmnuj48nO7mv1nrpDQCW0fY+Mtb3kjLWjifJwcFBvXfvXt/duLWpfx045jE2VF2OGfuLrinD6J7jmKEa0rFXa12qMzLLAADQQrAMAAAtrIYBAIzS+a/0lZ7QFcHymg2pFgeYLucarnLV+BBUwmqUYQAAQAuZZRghmSG4nauOIVl74DzBMgAwShIHbIJgGQCYlLZvBwTX3IRgGQC2yGWBpCAS2pngBwAALWSWAS647QSvVbJ0JpMBDJvMMgAAtJBZnrBVa9BkuACmoe38vy3n+S5rsC9uQ/Xeq7nJ9up73MosAwBAC5llANgSsqCwOpllYBzm82R/P7lzp7mfz/vuEbCFaq2P3LjExM7XZcw7+uDgoN67d6/vbjyi77qa88ZQszzm8cfNrTrWnk7yliRPnmv7UJLXJ3lmfd1aG6th3E7f54Wp75O+ty8TN58nh4fJ6enDtp2d5Pg4mc1u9Cu7OiZrrUv9YsHymg3pJCtYZqhWHWvvTbJ/Sfv9JC+/fXfWTrB8O32fF6a+T/revkzc/n5ycvJ4+95ecv/+jX5l38GymmVg8HZXbB+TqwKXqQdtU7WJYHQTY8OqD9zIgwertS+h7/OkmmVg8NpOsTc/9QKjN7G62MnYbUljtLWPgGAZGLw3pqlRPu9Di3ZgC53VxZ6cJLU294eHAuYhuHu3qVE+b2enaR8pwTIweM+kmcx3P8nzi/uhTu4DNuDo6NEJZEnz+Oion/7w0GzWTObb20tKae5vMblvCEzwW7Mh1Ria4MdQDek46cK6xvXUt1Obvs8Lt93uQ69ZXrZ/g65ZvnOnyShf8HySJxb/HlR/6cwtjwUT/GCqljk5DP0NG6bqJsfFZcdrV8fXkI/b6/r2wnba3b10xQXzGOiCMgwAYFwuqYs1j4GuCJZ5wcWrEm3iBmwv5wVu7JK62Ce/53vyN42jrbOJGEMZBgAwPrPZqCeNMR4yywAA0EKwDNvMov4AcCVlGLCtzhb1P1ur9GxR/8RXm8AgnV8tQ10ymyKzDD0rpax8WwuL+gMr6PR81JGNnU+ZNMEybKsHLSuStrUDwBZShrFhvjZiMFoW9c/u7ub7wmA4R3XHtr3eVdtIFpi+yCxPlYlbXOeSRf2zs9O0AwyM9bjpi2B5is4mbp2cJLU+nLglYOa8Sxb1z/GxyX3AqKlHZt3KmD+dHRwc1Hv37vXdjUcsfV37Lu3vX/71+t5ecv9+93+flXR1Qt/EWPNmdLl1bfs+tu+Y3xPWZSjH5FiOr6GNmVW329D6z/osEZMtNVhklqfIxC0AgLUwwW+KTNx6wVgyMwBjddV5VtaWKZBZniITtwDYUkOYBGgt52mRWV6zQXyKPpugdXTUlF7s7jaBsolbMCqDOJ/AiDmGtttV+//g4GDp3yOzPFWzWTOZ7/nnm3uBMgAwRR0vlyuzDACwRmcZzbMSjPOlGLLda3a2XO7pafP4bLncZG2JQpnloXExEQCA5RwdPQyUz5yeNu1rIrM8JBv4dHSVMa/pagLF48xQh2FxnrqcrCu3soHlcmWWh2QDn44AACajbVncNS6XK1geEhcTgWudXxZqqDcANmQDy+UKlodkA5+OAIDN8GF6A2az5Pg42dtLSmnuj4/XWr4qWB4SFxMBYMsIJLm1jpfLFSwPyQY+HQEAsDyrYQzNbCY4BgAYCMHyAFhOCABgmJRhAABAC5llJqFtYshUsvZDnfgy1H7BbVw1rpc5p6xyXEzhHOU8wNQJlgGArdP2QUXwz0XKMAAAoIVgGQAAWijDAAAm5XyJhbIKbkuwDADc2CYnKQp86YMyDNg283myv5/cudPcz+d99wi6Z9xPUinlkRt0QWaZzm1bJmDQr3c+Tw4Pk9PT5vHJSfM4ceVIpmuN476P47vvc8oYg9Blttmmt6vSkPGSWWZ4ZIC6c3T0MGA4c3ratMNUGffALcgsMywyn9168GC1dpgC4x64BZllhkUGqFu7u6u1wxQY95NVa33kNmRj6SePEywzLDJA3bp7N9nZebRtZ6dph6ky7oFbECwzLBPIAF2cnd3X7VKzWXJ8nOztJaU098fHSlyYtpGP+0GeS3oyln5uku3QPcEywyID1L3ZLLl/P3n++eZ+JAFDr0w6HT/jHrghE/xYi7XVYJ29gR0dNaUXu7tNoHzDN7YuasN8ct8yJp3SoavOUc41MAxlzIXmBwcH9d69e31349amcEIc8zha1Vj21zbtk07t7zcB8kV7e02GEjoylnPNUGzrOe/iONnW7bCqg4OD3Lt3b6mDTGYZ4ComnQIDJjjunpplgKtMYNIpADcnWGZrTWFmOBtg0ilLmvIqFLDNBMsAVxn5smMA3I6a5QlTxwRrMpsJjoGVmXw3DTLLAADQQrAMAAAtOg2WSymfVEp5eynln5VS3lVK+d2llE8ppfyDUsrPL+4/efHcUkr5llLKu0spP11KeWWXfQMA2DYmkq6u68zyX0vy92qtvynJ70jyriRvSPIjtdZXJPmRxeMk+bIkr1jcDpN8W8d9AwDoTK31kdu6WUllMzqb4FdK+cQkvzfJf5UktdYPJ/lwKeU1Sb548bS3JfnxJH8uyWuSfFdtRtM/WmSlP7PW+v6u+sjt3fagNNmhO1ftG9udPl133hjy+BSIsCnG2nB0mVn+7CTPJfnOUspPllK+vZTyZJLPOAuAF/efvnj+U0ned+7/P7toe0Qp5bCUcq+Ucu+5557rsPsAAGy7LoPlFyV5ZZJvq7V+fpIP5WHJxWUu+wj1WHqh1npcaz2otR685CUvWU9Pe3bxa5p13eAyTyfJ/n5y505zP5/32yGAJXnPuz3bbXVdBsvPJnm21vqOxeO3pwmef7GU8plJsrj/pXPPf9m5///SJL/QYf9g6zyd5C1JcnKS1NrcHx4KmAGgRWfBcq31XyZ5XynlcxdNr07yc0l+MMnrFm2vS/IDi3//YJI/slgV41VJPqheGdbrTUmevNh4epocHfXQGwBuQ5Z9M7q+gt9/l2ReSnlxkvck+Zo0AfrfLqV8bZIHSb5q8dwfTvLlSd6d5HTxXGCNdtt+8ODBJrsBAKPR6dJxtdafWtQX//Za61fUWj9Qa/1XtdZX11pfsbj/lcVza63162qt/2Gt9bfVWu912TfYRq0h8W5rGA3TM5+r2weW5gp+A9e2huIyN7jojWlm2j5iZye5e7eH3jBlgz03zedNnb66fWBJgmXYIs8keX2S7O0lpTT3x8fJbNZzz2BDjo6aOv3z1O0DV+i6ZhnokIkcsKK2+vwl6/avOuZ6z5pPiHObbTAkMssAbI+2+nx1+0ALwTIA2+Pu3aZO/zx1+9BqMPMNeiRYBmB7zGZNnb66fSZukBNsR0rNMgDbZTYTHANLEywDjFAfGaKhTjhaZlsMte/A8AmWAQC4lA+agmW22E1PAF1l9JyQgG3inMdYmOAHAAAtZJYBACZG5n59ZJYBAKCFzDIvsP4iNNSl98u5iD5ZXYWLZJYBAKCFzDKsSEYBhsmxOTz2CVMgswwAAC0EywAA0EKwDAAALQTLAADQwgS/CbP80nR1uW9NyOnOVfttDNvdOYUpuc14HvuxzGpklgEAoIXMMgCdkWWbPvuYqZNZBgCAFoJlAABoIVgGAIAWgmUAAGhhgh+s6Lrlhlad7DK05bja+mMSDzA2fZxfp7as3DLveRefM8bXeRWZZejZ00nem+Sji/un++0OAHCOYLkP83myv5/cudPcz+cb70KtdRC3rTef5y1J9tMcjPtJ3hIBM4zR1M9zU3990EawvGnzeXJ4mJycJLU294eHvQTMDMDRUZ680PRkkjf10RcA4DGC5U07OkpOTx9tOz1t2tk+Dx5c2ry74W4AwE1N/dsGwfKmtQRHre1M2+7lYbHRAMAYDG2SehcEy5vWEhy1tjNtd+/mQxeaPpTkjX30BQB4jGB50+7eTXZ2Hm3b2WnaSSnl2tsm/9ZN/v5Kv++rvzqvT3I/yfOL+9cneWZtr3I7bHLcAMMypuN+LP3kUYLlTZvNkuPjZG8vKaW5Pz5u2tlKzyR5eZInFvcCZQAYDhcl6cNstnRwPMVCedbH+NhebftelorbcE7hootjYhvPMYJlAACWso0fqJRhAABAC8EyAAC0UIbBtYZWn3RVf7bx66GhGsu4MWbG4bbjyX4et6GdT9bN+9qwySwDAEALmWWAC67K5Ew9w8X0yVTCamSWAQCghWAZAABaCJYBAFiv+TzZ30/u3Gnu5/O+e3RjapYBAFif+Tw5PExOT5vHJyfN42TpKxgPSRlzof/BwUG9d+9e393oVB+TicZ8actVxvOYXlebIR+/Y9m+q27Drl7Xuvbl0Lb70F/XJo6hqe6TLgxtWw3dYPfl/n4TIF+0t5dyWfsFmzpv1FqXGnDKMAAAWJ8HD1ZrHzjBMgxErXXl25D10e+xbaOp2MS+Hutx0Bfbi17t7q7WPnCCZQAA1ufu3WRn59G2nZ2mfYQEywAArM9slhwfJ3t7SSnN/fHxKCf3JVbDAABg3Waz0QbHF8ksAwBAC5llGIirlrgxKWc5lp2aLvt2OuzL9Rnq+8Zt93Hb/7/sNW1iPMksAwBAC8EyDNzTyWQuGQqsbmrLSsLYKMOAAXs6yVuSh1dCGvklQwFgbGSWYcDelOTJi42np8nRUQ+9AYDtI1iGAWu91tFILxkKAGMjWIYBaw2JR3rJUAAYG8EyDNgbk3zoYuPOTv7wyUlKKZ3cAG7LOaZ/U9z+fb0ewTIM2DNJXp88dsnQZ3ruFwBsC6thcK2rliEa+6fU6wx2Caav/uq+e9CpVbf71MfhKgY7ZreYfcK2mVrcILMMAAAtBMsAANBCsAwAAC3ULMMGjLFGC4Dpuer9SH395WSWAQCghcwyverzU6xP0P2x7VmVMTMuU1sNgfUZ49iQWQYAgBaCZQAAaCFYZiOeTvLeJB9d3D/db3cAAJaiZpnOPZ3kLUmeXDzeXzwGABg6wfKaXVecvo2TVN6Uh4HymScX7esu5t/G7QsA62BZucspw6Bzuyu2A1tsPk/295M7d5r7+bzvHgFbTmZ54No+yS2Tkd3Ep8Cl/sb+fnJy8ljzg/V3h3Mu2zdDXZZnTLY5u9K5+Tw5PExOT5vHJyfN4ySZzfrrF2s1xqXD2IzbxDxdklmme3fvJjs7jzR9KMkb++kNXZMZ5KaOjh4GymdOT5t2gJ4IlunebJYcHyd7e0kpuZ/k9Ume6blbdOAsM3hyktT6MDMoYGYZD1q+b2prB9iAMuavFA8ODuq9e/f67sYjNjXBbyhlGDfR5dcp2/ia29y0DONW27Cl5CZ7e8n9+zf/vde4zfYd6pjZSj2NH4aj76/badfnubKrcVFrXeoXyywD6yMzyG1cUrKVnZ2mHaAnJvhtmU1m4MacIZhCprKX7b+7e3lmcNfaJyzhbBLf0VHzAWt3twmUTe7r3BTOeXSrbYxsw/6XWQbWR2aQ25rNmpKL559v7gXKQM8Ey2tWa73ytqm/sw2f9LbJaPbthcmc2dtrHgt4gCV4P+Myfcc6yjCA9ZrNBMcATIbMMgAAtBAsAwBAC2UYazakdZbXbQyrWwxhtm4X22kM2x661OE6q5383iFyHqELV42rqRxfMssAANBCsAwAAC2uLMMopfzdJK059FrrH1x7jwBWNJWv+mDI+jzOrvrbyku22ybGxnU1y9909veSvCXJH13LXwUAgBG4Mliutf7Ds3+XUn7t/GMAAJi6VWqWfc8JAMBWua5m+VPOPXyilPLJaUoykiS11l/pqmMAbM6Q6j7VoN/ckPYjTGVZuetqlt+ZJqN89mp/4tzPapLP7qJTAAAwBNcFy19caz3ZSE+AlY3pkzkwfGM8p1gpgzZXjY2Dg4Olf891Nct/Z+nfBAAAE3NdsOwjGQAAW+u6MoynSinf0vbDWuufWHN/AABgMK4Llv9tmkl+AACwda4Llv9VrfVtG+kJAJ0y2Wk67EvGbkzLyl1Xs/zhjfQCAAAG6LrLXb/qsvZSyucm+TO11td30ivo2dA+1QKsyzad3ywrxzpcmVkupfz2UsrfL6X8TCnlL5VSPqOU8n1JfiTJz22miwAA0I/ryjDekuRvJvkvkjyX5gp+70nyObXWN3fcNwAA6NV1E/w+ttb6vy/+/c9LKX8myRtqrR/ttlsAANC/64LljyulfH4eXpzk15L89rIo9Km1/kSXnQMAgD5dFyy/P8k3n3v8L889rkm+pItOTdmYlkoBoD8moMEwXLcaxn+6qY4AAMDQXJdZTinl05N8XZLfkiab/HNJvrXW+ksd940RkRVa5fHeAAAgAElEQVSH4RvzMlrOMauxva7Xto2GfixMxZjG6HVLx31Rkn+yePhdSb5n8e9/vPgZAABM1nWZ5b+S5CtqrT95ru0HSil/J8n/luR3ddYzAADo2XXrLH/ihUA5SVJr/akkn9BNlwAAYBiuyyyXUson11o/cKHxU3J9oA30oI96u3XVnt2272OqgYNEfSyMwXUB75uT/P1Syn9SSvmExe2Lk/xfi58BAMBkXbd03HEp5ReS/MU0q2Ekyc8m+Uu11r/bdedoJ4MGgPeC9RvzqjFDM5Xxee3ScbXWH0ryQxvoCwAADMqVwXIp5X+44se11voX19wfAAAYjOsyyx+6pO3JJF+b5FPTlGcAAMAkXVez/FfO/l1K+YQkX5/ka5J8b5o1mAEAYLKWudz1pyT500lmSd6W5JUXl5IDxr1k221d9do32ce++3GbMTCUfcn6mRQG43ZdzfJfTvKfJzlO8ttqrb+2kV4BAMAAXJdZ/oYk/z7Jf5/k6Nyn45Jmgt8ndti3UZIdom+WPYJp8b4yHM6vj9uG8XldzbKr9AEAsLU6D4ZLKU+UUn6ylPJDi8cvL6W8o5Ty86WUv1VKefGi/WMXj9+9+Pl+130DAICrbCJz/PVJ3nXu8TcmeXOt9RVJPpBmGbos7j9Qa/2cNJfS/sYN9A0AAFp1GiyXUl6a5Pcn+fbF45LkS5K8ffGUtyX5isW/X7N4nMXPX122tQAIAIBBuHbpuFv6q0n+bJJPWDz+1CS/Wmv9yOLxs0meWvz7qSTvS5Ja60dKKR9cPP+XO+4jrMRnuHFq22+rTk7pav/3vewdAJfrLLNcSvkDSX6p1vrO882XPLUu8bPzv/ewlHKvlHLvueeeW0NPB2w+T/b3kzt3mvv5vO8eAQBslS4zy1+U5A+WUr48yccl+cQ0meZPKqW8aJFdfmmSX1g8/9kkL0vybCnlRUl+Y5JfufhLa63HadZ9zsHBwXTTLfN5cniYnJ42j09OmsdJMpv11y+ANZExf8i2GL+pLyu3zWO0s8xyrfXP11pfWmvdT/LaJD9aa50l+bEkX7l42uuS/MDi3z+4eJzFz3+0bvOeOTp6GCifOT1t2gEA2Ig+1lH+c0n+dCnl3Wlqkt+6aH9rkk9dtP/pJG/ooW/D8eDBau0AAKxd1xP8kiS11h9P8uOLf78nyRde8px/l+SrNtGfUdjdbUovLmsHAGAjXKFvqO7eTXZ2Hm3b2WnaAQDYiI1klrmBs0l8R0dN6cXubhMoLzG5b2gTCVYpPe+j79tcGs/wjpfLWFYO6Nu6lt/cpHWd3wXLQzabWfkCAKBHgmVYUdun6DFkKIdu6ksvDdmQs0NjZ9vifWPc1CwDAEALwTIAALQQLAMAQAvB8pqVUm58Axi8+TzZ30/u3Gnu5/O+e3QjzskwDZs4hk3wA2A583lyeJicnjaPT06ax4mVe4DJklnmMbXWpW9D6R+wAUdHDwPlM6enTTuwsqG9v15lLP3sgmAZgOU8eLBaO8AECJYBWM7u7mrtABMgWAZgOXfvJjs7j7bt7DTtABMlWAZgObNZcnyc7O0lpTT3x8cm9wGTVsZcnH1wcFDv3bvXdzcecZulSta1L267XMoq/Vjmb637963bULb7Taza9zH0sY2lvC435nP4JnQ1bmx3uja0c16fY77D43ipXyyzDAAALayzzKRc/OQ7tE/mN7XJT/QyZgD9u+pc7BuTzZJZBgCAFoJlAABoIVgGAIAWgmUAAGhhgh+MQNtkji4mY6x74ogJI3Stj4m8N/mbjoVpW/dSqttoqJPyZZYBAKCFzDK3su5PyUP/fWPXx1JEAKxX27lcdrsbMssAANBCsAwAAC0Ey2Mwnyf7+8mdO839fN53jwAAtoKa5aGbz5PDw+T0tHl8ctI8TpLZrL9+AQBsgTLmQu+Dg4N67969vrvxiNtMkrp0X+zvNwHyRXt7yf37a+/DZf0Y08SvPsfz0LbTZdtiDH1sM7S+D8WYz+GrmPr+35b9OHZdjcN17f8xTPAb0rFca12qM8owhu7Bg9XaAQBYG8HymtVab3y71O7uau237EPfnziBxw36eDWnAgbDe3w3BMtDd/dusrPzaNvOTtMO0KezORUnJ0mtD+dUCJiBCREsD91slhwfNzXKpTT3x8cm9wH9Ozp6OPn4zOlp0w4wESb4ca0hFeNfxwS/h0zwm47Bnqfv3GkyyheVkjz//Fr/1NT3/2D3MY8Y+gS/MRjSsbzsBD9Lx63ZmA+kIQ3gm2p7Ddt0IoKN2d29fLWeK+ZUwND18V541d8c4/vXFOKJ85RhAHAz5lQAW0BmGaDFGDM6G3U2d+LoqFnOcne3CZRHNKdi1X08tYwZcD3BMgA3N5uNKjgGWJUyDAAAaCFYBgCAFoJlAABooWYZgEHY5slzU1s6jPUZ6tjYpuNVZhkAAFrILK/ZTT7lbdOns67JwEzXVft2mWPI2Ji+Tezj245DYHxklgEAoIVgGQAAWgiWAQCghWCZ7TOfJ/v7yZ07zf183nePAICBEiyvWSll5RsbNJ8nh4fJyUlSa3N/eDjJgHm0Y82Hmcnr6ly4ifNsH+dy7x/QrzLmGeIHBwf13r17fXfjEV2dvDaxn6Zw4r12O+3vNwHyRXt7yf37t/rbU9h+m/bY/jr7MHN6+rBtZyc5Pk5ms9bfYzWMcenjWFnX/h/acW5cr9/Q9vFVrLN8O7XWpV6EYHnNxhwstxnaAXGrbXHnTpNRvqiU5Pnnb/57M7ztNAaP7csOP8wwHLc5Vi47/jf5YWlox/mY38PHaMwf9NZtaMfCTSwbLCvDYLvs7q7WzmY9eLBaOwB0TLDMdrl7t/la/7ydnaad/vkwA8DACJbZLrNZU/+6t9eUXuztXVsPywb5MAPAwLjcNdtnNhMcD9XZfjk6akovdnebQNn+AqAnJvitmQl+07LKdt/m7XRTYz7/rOK2Y2OM22ksx8NNJwwO3RjHzBiY4PfQRI4TE/wAAOA2lGFwrU1+qp3CJ9WbGmr2AIDGVefpdS+JOHRdbYshklkGAIAWgmUAAGghWAYAgBZqlkeirf5n1TonM3nZNlOrnWO9pjo+rnpdzsnbaZOXhZ8amWUAAGghs7xmN/lUNtXMxrbxiZzbMH5g3BzDD01tpQyZZQAAaCFYBgCAFoJlAABoIVgGAIAWJvit2RgL1yEZ1tg1UQZgdbc9j1ty8HIyywAAy5rPk/395M6d5n4+77tHdExmecsMfTmXbf7kyvRt0/ge+rnmtlbZl1N4vSzM58nhYXJ62jw+OWkeJ8ls1l+/RmSM5waZZQCAZRwdPQyUz5yeNu1MlmAZAGAZDx6s1s4kCJYBAJaxu7taO5MgWAYAWMbdu8nOzqNtOztNO5Nlgh9swFAnLcCZZcboNk1QvI4ltrbU2SS+o6Om9GJ3twmUNzy5r4/3lLa/uQ3jXbAMALCs2czKF1tGsMwLtuHTIWyCY+lqbdvHNzDrYwwyRkM9N6hZBgCAFoJlAABoIVgGAIAWguWReTrJe5N8dHHvmvQAAN0xwW9Enk7yliRPLh7vJ65Jv9BH8b8JNIzNbY+TbV46ahV9T0a6rg/2F+u0DWNNZnlE3pSHgfILXJMeAKAzMstr1umnqDt3kst+v2vSw8ZNJWMyJFdt0yFka2/rstc3hdfF5g11ibWu9H1ukFkeE9ekBwDYKMHymLgmPQDARgmWx2Q2S46Pk729pJTm/vh46yf3AQB0pYy57u7g4KDeu3ev725M3lRroG5rlWPHNlzNOs9LXW37oZ47hzbW1rWdhva6ttFQx/wYXDd+t/E4GcJrrrUu9Z9llgEAoIXVMGALbdOsfNmwaeh7NjyMwTYeJ5t4zTLLAADQQrAMAAAtBMsAANBCsAwAAC1M8NuQMSxfNdXifx5nX4/fWPbhVf00+XJc2val/di4zTG5jcfJmF6zzDIAG/F0kvcm+ejiPvN5r/0BWIbMMgCdezrJW5I8uXi8nySHh82DW16FdCzLZa2aLRtS3xm/sRwn63TVaz44OFj698gsA9C5N+VhoPyC09Pk6KiH3gAsT7AMQOd2237w4MEmuwGwMsEyAJ1rDYl3W8NogEEQLAPQuTcm+dDFxp2d5O7dHnoDsDwT/NZs00XyY1p6BZZlXDemNOnmmcX9m9KUZDxIsn98fOvJffRvm47XPo7JdW3fMZ1PhrZMocwyABvxTJKXJ3licS9QBsZAZhk2YBuX7IGk/8xi29+/7XF32e91LMM0ySwDAEALwTIAALQQLAMAQAvBMgAAtDDBDwA2wARA1m2blu1L+nu9MssAANBCZpkXdLXEEnAzY1lycIwZrC767BwK0ySzDAAALQTLAADQQrAMAAAt1CzzAnV1ADAN2/ae3uVKGTLLAADQQmaZXl38tLdtn4RhTMa46sUQjGVVE643hX152WsYS9/7IrMMAAAtBMsAANBCsAwAAC0EywBAv+bzZH8/uXOnuZ/P++4RvMAEP3plUgG0c3ywFebz5PAwOT1tHp+cNI+TZDbrr18T5byyOpllAKA/R0cPA+Uzp6dNOwyAzDIvaFsS57afQte9PJxlb5az6jJfQ9qGlijrj23Pqm49Zh48WK19oKa6rFybob2mLs9dMssAQH92d1drhw3rLFgupbyslPJjpZR3lVJ+tpTy9Yv2Tyml/INSys8v7j950V5KKd9SSnl3KeWnSymv7KpvAMBA3L2b7Ow82raz07TDAHSZWf5Ikm+otf7mJK9K8nWllM9L8oYkP1JrfUWSH1k8TpIvS/KKxe0wybd12DcAYAhms+T4ONnbS0pp7o+PTe5jMDqrWa61vj/J+xf//jellHcleSrJa5J88eJpb0vy40n+3KL9u2pTdPKPSimfVEr5zMXvAQCmajYTHDNYG5ngV0rZT/L5Sd6R5DPOAuBa6/tLKZ++eNpTSd537r89u2gTLI/c0CYBwKY5BoAxuOpctc2Tfzuf4FdK+Q1Jvi/Jn6y1/uurnnpJ22N7ppRyWEq5V0q599xzz62rmwAA8JhOM8ullI9JEyjPa63fv2j+xbPyilLKZyb5pUX7s0ledu6/vzTJL1z8nbXW4yTHSXJwcLC9H3Mg6/2k39XSgTf5m7AtxrTcmOP15m6y7Ya2/y/qY+z2NQa7XA2jJHlrknfVWr/53I9+MMnrFv9+XZIfONf+RxarYrwqyQfVKwMA0KcuM8tflOS/TPJPSyk/tWh7Y5L/JcnfLqV8bZIHSb5q8bMfTvLlSd6d5DTJ13TYNwAAuFaXq2H8v7m8DjlJXn3J82uSr+uqPwAANzKfN5fffvCguVjK3btW79giLncNANBmPk8OD5PT0+bxyUnzOBEwbwmXu4YbKqVcehtCH4Y+MYTh2pYxddXrnOLr5RaOjh4GymdOT5v2WxrT+BtLP7sgWAYAaPPgwWrtTI4yjAG5bEmUKXxqu/i6pvCatsltlwey3NS4lgdjOLZpaa5B291tSi8ua+dSU3vfkFkGAGhz926ys/No285O085WECwDALSZzZLj42RvLymluT8+NrlviyjDAAC4ymwmON5igmWYqKvqwoZWD9aHMdbNDc0YtuFtansdQ8Nx2xpt+6s7Q553sa6+KcMAAIAWMsvAjcnW3N7QVsqwT8fBfpqGoR3/V1n3ylZjGsMyywAA0EKwDAAALQTLAADQQrAMAAAtTPCbsKEsezS0SQpsF0uH3c5tj99NbMM+zjHGRje62pf21+2te9+sa59s4viXWQYAgBYyy2vW9mlIdpU2shq06WpZKWMOhq+PZeUu+5viF5llAABoJVgGAIAWgmUAAGghWAYAgBaC5QEppTx2g17N58n+fnLnTnM/n/fdo0tddux0fQz18Te71MdrWeVvjml7j6WffRnavhxSX/o2lO0wtH1iNQzgcvN5cniYnJ42j09OmsdJMpv11y8A2KAy5iWEDg4O6r179/ruxlKG9il1aAt+r2LMY3ZU9vebAPmivb3k/v1N9+ZKxujtDWkbTnX5qjGOiy6MZV+OcX8ts23HfHytu++11qX+szIM4HIPHqzWDgATJFgGLre7u1o7AEyQYBm43N27yc7Oo207O007AGwJwTJwudksOT5uapRLae6Pj03uA2CrCJaBdrNZM5nv+eeb+x4D5aEtN3WVofZzLNtw6P27qbFs/3XYptd6qZEsuzk2fY0nS8cBAKyLZTcnx9JxG9LHp+kx71u4aAoZqb6Pyaluw6m+rjHb6n0yomU3zwxtf13c9l0tiWfpOACATbPs5uQIlgEA1sWym5MjWAYAWBfLbk6OCX4AG3JV/dzUala7MrTaynUZ49iY6r44c+N9cjaJ7+ioKb3Y3W0CZZP7lnaTsdXleDTBb0NM8IPbmfob8yaO16lvw6ka6rl8m8fTUPfJTW3rvlx2gp/M8oZM7cCCTbvqGNrWE/2qpr4Np7pSBnStj3PDJo/Xy/7WwcHB0v9fzTIAALQQLAMAQAvBMgAAtBAsAwBACxP84BK3mWRgMuf1lrgE6YZ6MhxjXDpsaEzmY1Mcr7e3yeP1tn9LZhkAAFrILAOjN/Ul0TahbRsObfvddLmpi/9vaK9rjBx32+E254apLOcoswwAAC0EywAA0EKwDAAALQTLAADQwgS/keqyQL7rZW/GWNzP6m6zny3LRJubjivnHTbF+euhqRx3MssAANBCZhmYtDEviTYEfSwPtsnlpoaytNVQ9/9Nrfv1DO14ZbuWDpRZBgCAFoJlAABoIVgGAIAWapbXbGp1OnCVPsa7mebTt8lxNbVz9k0vQcxmtO0f+2TYZJYBAKCFYBkAAFoowwC20lCWRBuz227Di/9/aCUR615Wbmr7f5O62HZDG29TMrVl5WSWAQCghWAZAABaCJYBAKCFmmVgMrZpWbnr6v7G8HqHXrs45P7dpm/bdJyMhX0ybDLLAADQQmaZjRvaLFmf2m9uaPtyXda9ysM2murYOGMfj1/bPpzC+ByyMZ4bZJYBAOjOfJ7s7yd37jT383nfPVqJYBnYXiM/gQMM3nyeHB4mJydJrc394eGozreCZWA7TeAEDjB4R0fJ6emjbaenTftIlDHXXR0cHNR79+713Y1HDLXeZhV9jgk1y9MxtGPhsf28v98EyBft7SX377f+nqHULI95NYyhjY2b6Hv7uspkd4Y2Pke/T+7caRISF5WSPP/8habNbvta61J/0AQ/tp4le7bDxf380Vz+1drzJyd5YvHcvve/5cG2Wx9Bm3HD2u3uXp6Y2N3dfF9uSBkGsJUerNgOwA3cvZvs7DzatrPTtI+EzPKa3eST922yB5f9vaF9hbSKm2Yuxvyap2roywO9Mclbkjx5ru1Di/aryK7dniW7GDLvQ2s2mzX3R0fJgwdNRvnu3Yft53R1brjs9x4cHCz9/wXLwFZ6ZnH/piS7aTLKbzzXDsCazGaXBsdjIVgGttYzERwDcDU1ywAA0EKwDAAALZRhwBZZ92TSqbvp9lplW5kUtN3GtP/b+rqN5wa2i8wyXOLpJO9NXAYZALaczPIA+FQ+LE/n3JJi5y+DnIx6Nu+QWDqMywxtuUHnZlZhvHSn720rswwXvCmPrr2bZHTXsQcA1kOwDBe0XoDzgWu7AcC2ESzDBa0h8YiuYw8ArIeaZbjgsssgj+k69l3Vdl71e/uuJxuaMdZeWymlO2McD8BDMstwwTNJXp/kfpKUkuztJcfHJvcBwBaSWWYSZLamYWirIYzZ1I6JrsbG1LYTsH4yywAA0EKwDAAALQTLAADQQrAMAAAtTPAbKZOdttvQ9r9l5cbBsoLdGtpxuSn2P1MnswwAAC1kloFRGNOychf72mX/ZO7Gsw0u6+fQxi7wOJllAABoIVgGAIAWgmUAAGghWAYAgBYm+G2ISRzcxBTGjWWlNmtoY2Zq+/8223do+wZYjswyADBM83myv5/cudPcz+d994gtJLMMjF4fy8qtkhW9bf/GmIFlNWNeVq6z8TmfJ4eHyelp8/jkpHmcJLNZN38TLiGzDAAMz9HRw0D5zOlp0w4bJFiG2/AVIUA3HjxYrR06Ilhm+roKaM++Ijw5SWp9+BWhgBng9nZ3V2uHjgiWmbYuA1pfEQJ05+7dZGfn0badnaYdNkiwvGallEtvYzKF1/CCLgPaNX1F2La9R7vNYUlDHfd9HJND3A69m82S4+Nkby8ppbk/Pja5j42zGgbT1mXN2+5uk6m+rB2A25vNBMf0TmaZaeuy5s1XhKNQa730dpv/u86lsq76Gzf5W5voM5s1lH1pXLGtBMtMW5cBra8IAWDylGEwbWeB69FRU3qxu9sEyusKaH1FCACTJrPM9M1myf37yfPPN/eCWwAuY+18LiGzDAyGSz/D46yMsSEur00LmWUAAGvn00JmGejUUDPBQ+1Xl656zbKXD41pO13s6237t43HxQtcXpsWMssAAC6vTQvBMgCAtfNpIVgGALB2Pi3ULAMAJNbO51KCZWDjbjMJ6ar/u9WTkxYsv7ddhjbhEKZIGQYAALSQWYaejWmZqjYylePXtg+7HIPGTT9sd1iNzDIAALQQLAMADMl8nuzvJ3fuNPfzed892mrKMAAAhuKP//Hkb/yN5Kxc5uQkOTxs/m2ljl7ILAMADMF8/migfOb0NDk66qdPCJaBbpRSWm9T+pt9uc3r3KbtBKNydPR4oHzmwYPN9uWMkhBlGAAAg3BVQLy7u7l+nJnPmxKQ09Pm8ZaWhMgsT1itda03Nm9I+8LY2E633e9TGzdDez3O5RPTFhCXkty9u9m+JE2m+yxQPrOFJSGCZQCAIbh7N9nZebStlOSP/bHLM7ldl0i0Zbr7KgnpyaCC5VLKl5ZS/nkp5d2llDf03R8AgI2ZzZLj42RvrwmS9/aS7/7u5K//9cefe1YicXLS1DmflUisM2Buy3T3URLSo8EEy6WUJ5J8a5IvS/J5SZ4upXxev70CANig2Sy5fz95/vnmvq02eBMlEpdlund2+ikJ6dFgguUkX5jk3bXW99RaP5zke5O8puc+AQAMzyZKJC7LdB8fb9XkvmRYq2E8leR95x4/m+R39dQXAIDh2t1tSi8ua1+n2WzrguOLhhQsX7bA52PTdksph0kOk2R3gDUzZhqzCWMYZ2Po45gNYfsOoQ9DYnuwUXfvPrqsW7KVJRI39c53vvOdyz53SGUYzyZ52bnHL03yCxefVGs9rrUe1FoPXvKSl2yscwAAg6FEYmOGlFn+J0leUUp5eZJ/keS1Sf5wv10CABgoJRIbMZhgudb6kVLKf5vk/07yRJLvqLX+bM/dAgBgiw0mWE6SWusPJ/nhvvsBAADJsGqWAQBgUATLAADQQrAMAAAtBMsAANBCsAwAAC0EywAA0EKwDAAALQTLAADQQrAMAAAtBMsAANBCsAwAAC0EywAA0EKwDAAALQTLAADQQrAMAAAtBMsAANBCsAwAAC0EywAA0EKwDAAALQTLAADQQrAMAAAtBMsAANBCsAwAAC0EywAA0EKwDAAALQTLAADQQrAMAAAtSq217z7cWCnluSQnfffjgk9L8st9d4JBMjZoY2xwGeOCNsbG7e3VWl+yzBNHHSwPUSnlXq31oO9+MDzGBm2MDS5jXNDG2NgsZRgAANBCsAwAAC0Ey+t33HcHGCxjgzbGBpcxLmhjbGyQmmUAAGghswwAAC0Ey2tUSvnSUso/L6W8u5Tyhr77w+aUUl5WSvmxUsq7Sik/W0r5+kX7p5RS/kEp5ecX95+8aC+llG9ZjJWfLqW8st9XQNdKKU+UUn6ylPJDi8cvL6W8YzE2/lYp5cWL9o9dPH734uf7ffabbpVSPqmU8vZSyj9bnD9+t/MGpZQ/tXgv+ZlSyjOllI9zzuiPYHlNSilPJPnWJF+W5POSPF1K+bx+e8UGfSTJN9Raf3OSVyX5usX+f0OSH6m1viLJjyweJ804ecXidpjk2zbfZTbs65O869zjb0zy5sXY+ECSr120f22SD9RaPyfJmxfPY7r+WpK/V2v9TUl+R5ox4ryxxUopTyX5E0kOaq2/NckTSV4b54zeCJbX5wuTvLvW+p5a64eTfG+S1/TcJzak1vr+WutPLP79b9K84T2VZgy8bfG0tyX5isW/X5Pku2rjHyX5pFLKZ26422xIKeWlSX5/km9fPC5JviTJ2xdPuTg2zsbM25O8evF8JqaU8olJfm+StyZJrfXDtdZfjfMGyYuSfHwp5UVJdpK8P84ZvREsr89TSd537vGziza2zOIrsM9P8o4kn1FrfX/SBNRJPn3xNONlu/zVJH82yfOLx5+a5FdrrR9ZPD6//18YG4uff3DxfKbns5M8l+Q7FyU6315KeTLOG1ut1vovknxTkgdpguQPJnlnnDN6I1hen8s+xVlqZMuUUn5Dku9L8idrrf/6qqde0ma8TFAp5Q8k+aVa6zvPN1/y1LrEz5iWFyV5ZZJvq7V+fpIP5WHJxWWMjS2wqFF/TZKXJ/msJE+mKcG5yDljQwTL6/Nskpede/zSJL/QU1/oQSnlY9IEyvNa6/cvmn/x7GvSxf0vLdqNl+3xRUn+YCnlfpryrC9Jk2n+pMVXrMmj+/+FsbH4+W9M8iub7DAb82ySZ2ut71g8fnua4Nl5Y7v9viTvrbU+V2v99STfn+T3xDmjN4Ll9fknSV6xmK364jTF+D/Yc5/YkEV92FuTvKvW+s3nfvSDSV63+PfrkvzAufY/spjd/qokHzz72pVpqbX++VrrS2ut+2nOCz9aa50l+bEkX7l42sWxcTZmvnLxfFmiCaq1/ssk7yulfO6i6dVJfi7OG9vuQZJXlVJ2Fu8tZ+PCOaMnLkqyRqWUL0+TMXoiyXfUWu/23CU2pJTyHyf5f5L80zysS31jmrrlv51kN80J8Ktqrb+yOAH+r0m+NMlpkq+ptd7beMfZqFLKFyf5M7XWP1BK+ew0meZPSfKTSb661vrvS0d3OvwAAAIZSURBVCkfl+S709S9/0qS19Za39NXn+lWKeV3ppn4+eIk70nyNWkSWc4bW6yU8j8l+UNpVlr6ySR/NE1tsnNGDwTLAADQQhkGAAC0ECwDAEALwTIAALQQLAMAQAvBMgAAtBAsA4xEKeWjpZSfOnd7w7mfvaSU8uullP/mwv/5r0sp/7SU8tOllJ8ppbymlPKti///c6WUf3vu933l438VYLtZOg5gJEopv1Zr/Q0tP/vjSZ5O8tFa6xcv2l6a5B8meWWt9YOLy7G/pNb63sXP95P8UK31t26g+wCjJLMMMA1PJ/mGJC8tpTy1aPv0JP8mya8lSa31184CZQCWI1gGGI+Pv1CG8YeSpJTysiT/Qa31H6e58tsfWjz//0vyi0neW0r5zlLKf9ZPtwHG60V9dwCApf3bWuvvvKT9tWmC5KS5HO5bk3xzrfWjpZQvTfIfJXl1kjeXUr6g1vo/bqS3ABOgZhlgJNpqlkspP5HkM5L8+qLps5L8llrrz1943kGS76y1/rbF4/2oWQa4kjIMgBErpXxukidrrU/VWvdrrftJ/uckry2lfFYp5ZXnnv47k5z00U+AsVKGATAeH19K+alzj/9ekn+X5O9ceN73pSnH+K4k31RK+azF855L8sc20VGAqVCGAQAALZRhAABAC8EyAAC0ECwDAEALwTIAALQQLAMAQAvBMgAAtBAsAwBAC8EyAAC0+P8BVKc8c4py+NgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8ce340c9e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "plt.imshow(grid, cmap='Greys', origin='lower')\n",
    "\n",
    "nmin = np.min(data[:, 0])\n",
    "emin = np.min(data[:, 1])\n",
    "\n",
    "# draw points\n",
    "all_pts = np.array(to_keep)\n",
    "north_vals = all_pts[:,0]\n",
    "east_vals = all_pts[:,1]\n",
    "plt.scatter(east_vals - emin, north_vals - nmin, c='red')\n",
    "\n",
    "plt.ylabel('NORTH')\n",
    "plt.xlabel('EAST')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Solution](/notebooks/Random-Sampling-Solution.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Epilogue\n",
    "\n",
    "You may have noticed removing points can be quite lengthy. In the implementation provided here we're naively checking to see if the point collides with each polygon when in reality it can only collide with one, the one that's closest to the point. The question then becomes \n",
    "\n",
    "\"How do we efficiently find the closest polygon to the point?\"\n",
    "\n",
    "One such approach is to use a *[k-d tree](https://en.wikipedia.org/wiki/K-d_tree)*, a space-partitioning data structure which allows search queries in $O(log(n))$. The *k-d tree* achieves this by cutting the search space in half on each step of a query.\n",
    "\n",
    "This would bring the total algorithm time down to $O(m * log(n))$ from $O(m*n)$.\n",
    "\n",
    "The scikit-learn library has an efficient implementation [readily available](http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KDTree.html#sklearn.neighbors.KDTree)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}