{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "voronoi.ipynb\n",
    "\n",
    "Discussion: This Jupyter notebook investigates 1D and 2D Voronoi diagrams\n",
    "\n",
    "Licensing: This code is distributed under the GNU LGPL license.\n",
    "    \n",
    "Modified: 23 October 2016\n",
    "\n",
    "Author: John Burkardt, Lukas Bystricky"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: agg\n"
     ]
    }
   ],
   "source": [
    "# Import necessary libraries and set plot option\n",
    "%matplotlib\n",
    "%config InlineBackend.figure_format = 'png'\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import scipy.spatial as spatial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voronoi\n",
    "\n",
    "In this module we will investigate Voronoi diagrams in 1D and 2D.\n",
    "\n",
    "We are given a geometric region R (the whole line or an interval, the plane, \n",
    "a rectangle, circle or other shape).\n",
    "\n",
    "We are also given a set G of \"generators\", points inside the region.\n",
    "\n",
    "The Voronoi diagram simply assigns every point x in R to the nearest\n",
    "generator g in G, dividing R into Voronoi subregions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0.0, 0.30000000000000004, 0.15000000000000002)\n",
      "(1, 0.30000000000000004, 0.5, 0.40000000000000002)\n",
      "(2, 0.5, 0.69999999999999996, 0.59999999999999998)\n",
      "(3, 0.69999999999999996, 1.0, 0.84999999999999998)\n",
      "\n",
      " I            G[I]          Left[I]       Right[I]     Centroid[I]       Length[I]       Energy[I]\n",
      " 0        0.200000        0.000000        0.300000        0.150000        0.300000        0.003000\n",
      " 1        0.400000        0.300000        0.500000        0.400000        0.200000        0.000667\n",
      " 2        0.600000        0.500000        0.700000        0.600000        0.200000        0.000667\n",
      " 3        0.800000        0.700000        1.000000        0.850000        0.300000        0.003000\n"
     ]
    }
   ],
   "source": [
    "## A 1D finite region R.\n",
    "#\n",
    "#  Let R be [0,1], and define G = [ 1/5, 2/5, 3/5, 4/5 ]\n",
    "#\n",
    "g = np.array ( [ 1.0/5.0, 2.0/5.0, 3.0/5.0, 4.0/5.0] )\n",
    "#\n",
    "#  Compute the length and centroid of each subregion.\n",
    "#\n",
    "#  Note that:\n",
    "#    subregion 0 extends from 0 to midway between 1/5 and 2/5.\n",
    "#    subregion 1 extends from midway between 1/5 and 2/5, to midway between 2/5 and 3/5.\n",
    "#\n",
    "length = np.zeros ( 4 )\n",
    "centroid = np.zeros ( 4 )\n",
    "right = np.zeros ( 4 )\n",
    "left = np.zeros ( 4 )\n",
    "for i in range ( 0, 4 ):\n",
    "  if ( i == 0 ):\n",
    "    left[i] = 0.0\n",
    "  else:\n",
    "    left[i] = ( g[i-1] + g[i] ) / 2.0\n",
    "  if ( i < 3 ):\n",
    "    right[i] = ( g[i] + g[i+1] ) / 2.0\n",
    "  else:\n",
    "    right[i] = 1.0\n",
    "  length[i] = right[i] - left[i]\n",
    "  centroid[i] = ( right[i] + left[i] ) / 2.0\n",
    "#\n",
    "#  The \"energy\" of the i-th Voronoi subregion pi is the integral\n",
    "#  ei = Integral(pi) ( x-gi )^2 dx\n",
    "#  Compute the energy for each region.\n",
    "#\n",
    "energy = np.zeros ( 4 )\n",
    "for i in range ( 0, 4 ):\n",
    "  energy[i] = ( ( right[i] - g[i] ) ** 3 - ( left[i] - g[i] ) ** 3 ) / 3.0\n",
    "#\n",
    "#  Print i, g[i], left[i], right[i], centroid[i], length[i], energy[i]\n",
    "#\n",
    "print ( '' )\n",
    "print ( ' I            G[I]          Left[I]       Right[I]     Centroid[I]       Length[I]       Energy[I]' )\n",
    "for i in range ( 0, 4 ):\n",
    "  print ( '%2d  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f' \\\n",
    "  % ( i, g[i], left[i], right[i], centroid[i], length[i], energy[i] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " I            G[I]          Left[I]       Right[I]     Centroid[I]       Length[I]       Energy[I]\n",
      " 0        0.150000        0.000000        0.300000        0.150000        0.300000        0.002250\n",
      " 1        0.400000        0.300000        0.500000        0.400000        0.200000        0.000667\n",
      " 2        0.600000        0.500000        0.700000        0.600000        0.200000        0.000667\n",
      " 3        0.850000        0.700000        1.000000        0.850000        0.300000        0.002250\n"
     ]
    }
   ],
   "source": [
    "# Energy Reduction\n",
    "#\n",
    "#  We will see later that reducing the total energy is an important guide\n",
    "#  for improving the smoothness of a point distribution.\n",
    "#\n",
    "#  Notice the G[0] and G[3] are not the same as the centroids of their subregions.\n",
    "#\n",
    "#  Replace G[0] and G[3] by their centroids, and recompute the energies.\n",
    "#\n",
    "#  Expect the total energy to go down.\n",
    "#\n",
    "g[0] = centroid[0]\n",
    "g[3] = centroid[3]\n",
    "\n",
    "for i in range ( 0, 4 ):\n",
    "  energy[i] = ( ( right[i] - g[i] ) ** 3 - ( left[i] - g[i] ) ** 3 ) / 3.0\n",
    "#\n",
    "#  Print i, g[i], left[i], right[i], centroid[i], length[i], energy[i]\n",
    "#\n",
    "print ( '' )\n",
    "print ( ' I            G[I]          Left[I]       Right[I]     Centroid[I]       Length[I]       Energy[I]' )\n",
    "for i in range ( 0, 4 ):\n",
    "  print ( '%2d  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f' \\\n",
    "  % ( i, g[i], left[i], right[i], centroid[i], length[i], energy[i] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " I            G[I]          Left[I]       Right[I]     Centroid[I]       Length[I]       Energy[I]\n",
      " 0        0.200000        0.000000        0.300000        0.144178        0.320000        0.003303\n",
      " 1        0.400000        0.300000        0.500000        0.404088        0.190000        0.000649\n",
      " 2        0.600000        0.500000        0.700000        0.598942        0.200000        0.000724\n",
      " 3        0.800000        0.700000        1.000000        0.850480        0.290000        0.002969\n"
     ]
    }
   ],
   "source": [
    "#  Sampling estimates\n",
    "#\n",
    "#  Geometry in 1D is easy.  But let's pretend it's hard.  It turns out we can estimate\n",
    "#  the values of length, centroid and energy for each subregion by sampling.\n",
    "#\n",
    "#  If we can generate n uniform random samples of the area, and neari is the number of\n",
    "#  those samples xj closest to generator gi, and l is the length of our region, then:\n",
    "#\n",
    "#  estimated length of subregion i: li = l * ( neari / n )\n",
    "#  estimated centroid:              ci = sum ( xj ) / neari\n",
    "#  estimated energy:                ei = li * sum ( xj - gi )**2 / neari\n",
    "#\n",
    "#  Do these estimates for our 1D problem, using N = 1000.\n",
    "#\n",
    "#  Don't forget to restore the original values of g[0] and g[3]!\n",
    "#\n",
    "g[0] = 0.2\n",
    "g[3] = 0.8\n",
    "\n",
    "n = 1000\n",
    "x = np.random.random ( n )\n",
    "\n",
    "centroid = np.zeros ( 4 )\n",
    "length = np.zeros ( 4 )\n",
    "energy = np.zeros ( 4 )\n",
    "near = np.zeros ( 4 )\n",
    "\n",
    "for j in range ( 0, n ):\n",
    "  i = ( np.abs ( g - x[j] ) ).argmin ( )\n",
    "  near[i] = near[i] + 1\n",
    "  centroid[i] = centroid[i] + x[j]\n",
    "  energy[i] = energy[i] + ( g[i] - x[j] ) ** 2\n",
    "\n",
    "for i in range ( 0, 4 ):\n",
    "  length[i] = 1.0 * float ( near[i] ) / float ( n )\n",
    "  centroid[i] = centroid[i] / float ( near[i] )\n",
    "  energy[i] = length[i] * energy[i] / float ( near[i] )\n",
    "# \n",
    "#  Print i, g[i], left[i], right[i], centroid[i], length[i], energy[i]\n",
    "#\n",
    "print ( '' )\n",
    "print ( ' I            G[I]          Left[I]       Right[I]     Centroid[I]       Length[I]       Energy[I]' )\n",
    "for i in range ( 0, 4 ):\n",
    "  print ( '%2d  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f' \\\n",
    "  % ( i, g[i], left[i], right[i], centroid[i], length[i], energy[i] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.2  0.2]\n",
      " [ 0.4  0.2]\n",
      " [ 0.6  0.2]\n",
      " [ 0.8  0.2]\n",
      " [ 0.2  0.4]\n",
      " [ 0.4  0.4]\n",
      " [ 0.6  0.4]\n",
      " [ 0.8  0.4]\n",
      " [ 0.2  0.6]\n",
      " [ 0.4  0.6]\n",
      " [ 0.6  0.6]\n",
      " [ 0.8  0.6]\n",
      " [ 0.2  0.8]\n",
      " [ 0.4  0.8]\n",
      " [ 0.6  0.8]\n",
      " [ 0.8  0.8]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqAAAAILCAYAAAA35y7oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X90VPWd//HXEDJJ1CyRaQxDsNsmkTIBupUoPQieAIWu\nG6JVJLoRtmpoyaLgIW4lCJYfioWvLYuIft0g5Zx6sGnlx54DVhIMRipdPdXkYDdMVr4Be6J2Ekj4\nYdhOSJrk+0ckkp+Tm0w+k5s8H+d4Qj4zd+Z973mZvHJvZuIoLi5uEQAAAGDIiFAPAAAAgOGFAgoA\nAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADBqZKgH6I0LFy7oww8/1JgxY+R0OkM9\nDgAAADpoaGhQVVWVbrnlFsXExPR4X1sU0A8//FDPPvtsqMcAAABAAGvWrNGcOXN6vI8tCuiYMWMk\nSbt375bH4wnxNMGRk5OjrVu3hnqMIaW8vFyLFi0aUjkZLMhr8JHXgUNeg4+8DpyhlNcrObnS23pi\niwJ65bK7x+PRlClTQjxNcIwaNWrI7MtgM5RyMliQ14FDXoOPvA4c8hp8QzGvvfl1SV6EBAAAAKMo\noAAAADCKAgoAAACjKKAhkpmZGeoRgF4jr7AT8go7Ga55pYCGyHAN3EByu91at26d3G53qEcZcshr\n8JHXgUNeg4+8DpzhmldbvAoe6A23263169eHegygV8gr7IS8Itg4AwoAAACjKKAAAAAwigIKAAAA\noyigAAAAMIoCCgAAAKMooAAAADCKAoohw+/368SJE/L7/aEeBQiIvMJOyCuCjQKKIaO8vFyTJk1S\neXl5qEcBAiKvsBPyimCjgAIAAMAoCigAAACMooACAADAKAooAAAAjKKAAgAAwCgKKAAAAIyigAIA\nAMCokaEeAAgWj8ejsrIyJSQkhHoUICDyCjshrwg2CiiGjKioKE2cODHUYwC9Ql5hJ+QVwcYleAAA\nABhFAQUAAIBRFFAAAAAYRQEFAACAURRQAAAAGEUBBQAAgFEUUAwZPp9P69evl8/nC/UoQEDkFXZC\nXhFslgtoY2Oj8vLylJGRoTvuuEOPPPKISkpKerVtSUmJ/u3f/k333HOP0tLStHjxYu3fv1/Nzc2W\nBwc68vl82rBhA18gYQvkFXZCXhFslgvo5s2btW/fPs2dO1fLly9XWFiYVq1apbKysh63++Mf/6iV\nK1fqwoULWrhwoR555BHFx8frxRdf1Msvv9znHQAAAIC9WCqg5eXlKi4u1o9//GMtWbJE8+bN05Yt\nWxQXF6e8vLwet33rrbc0cuRIbdu2TQsWLFB6erqefvppffvb31ZBQUG/dgIAAAD2YamAHj16VGFh\nYUpPT29bczqdSktLk9fr1dmzZ7vd1ul0yul06rrrrmu3Pnr0aEVERFgcGwAAAHZlqYCeOnVK48aN\nU1RUVLv1CRMmSJIqKiq63faee+5RS0uLtmzZosrKSlVXV+vAgQM6duyYFi5c2IfRAQAAYEeWCmht\nba1cLlendZfLpZaWFtXW1na7bVJSkrZs2aL33ntPDz30kDIzM7V9+3YtX75c99xzj/XJbaq6Wpox\nQ0pMbP145kyoJwK6R15hJ+QVdjLc8zrSyp0vX76s8PDwTutOp7Pt9u5UVlbqySefVFxcnLKzs+V0\nOnXkyBG98MILGj16tKZPn25xdHu6917pD39o/ffp09L8+dKxY6GdCegOeYWdkFfYyXDPq6UCGhER\nocbGxk7rDQ0Nbbd35+WXX1ZYWJief/75tvulpqbq8ccf17Zt2zRt2jSNGNHzCdmcnByNGjWq3Vpm\nZqYyMzOt7EZIdXwHC97RIngiIyOVnJysyMjIUI8yZJDXgUNeg4+8DhzyGnx2z2t+fr7y8/PbrV28\neLHX21sqoC6XSzU1NZ3Wr1x67+ry/BVlZWWaNm1ap5J622236eWXX1ZVVZXGjh3b4/Nv3bpVU6ZM\nsTLyoON2t/6kc/XnCI7k5GSdOHEi1GMMKeR14JDX4COvA4e8Bp/d89rVCcDS0lKlpKT0antLBTQx\nMVHHjx+X3+9v90Ikr9crh8OhpKSkbrf929/+pqampi7XJXV521C0f3/raXafrzVs+/eHeiKge+QV\ndkJeYSfDPa+WCmhqaqpef/11HTx4UPfdd5+k1r+MVFhYKI/Ho9jYWEnSuXPndOnSJcXHxyssLEyS\ndNNNN6mkpER1dXWKjo6WJDU3N6u4uFhRUVEBz34OFTfcMLx+xwP2Rl5hJ+QVdjLc82qpgHo8HqWm\npmrnzp06f/684uPjVVBQoOrqaq1cubLtfjt27NDhw4eVn5+vuLg4SdKiRYv05JNPaunSpUpPT1dE\nRISOHDmiiooKLV68uK2oAgAAYGizVEAlafXq1dq1a5eKiopUV1enhIQEbdq0SZMnT267j8PhkMPh\naLfd1KlTtXnzZu3evVuvvvqqmpqadOONN+rxxx/XvHnz+r8nAAAAsAVHcXFxS6iHCOTkyZPKzs5W\nSUmJ7V+EBAAAMBRdeRFSXl6exo8f3+N9Lb0RPQAAANBfFFAAAAAYRQHFkOH1ejVx4kR5vd5QjwIE\nRF5hJ+QVwUYBxZBRX18vr9er+vr6UI8CBEReYSfkFcFGAQUAAIBRFFAAAAAYRQEFAACAURRQAAAA\nGEUBBQAAgFEUUAAAABhFAcWQ4Xa7tW7dOrnd7lCPAgREXmEn5BXBNjLUAwDB4na7tX79+lCPAfQK\neYWdkFcEG2dAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAMWT4/X6dOHFC\nfr8/1KMAAZFX2Al5RbBRQDFklJeXa9KkSSovLw/1KEBA5BV2Ql4RbBRQAAAAGEUBBQAAgFEUUAAA\nABhFAQUAAIBRFFAAAAAYRQEFAACAURRQAAAAGDUy1AMAweLxeFRWVqaEhIRQjwIERF5hJ+QVwUYB\nxZARFRWliRMnhnoMoFfIK+yEvCLYuAQPAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAA\nAAAwigKKIcPn82n9+vXy+XyhHgUIiLzCTsgrgo0CiiHD5/Npw4YNfIGELZBX2Al5RbBRQAEAAGAU\nBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQDFkREZGKjk5WZGRkaEeBQiIvMJOyCuC\nbWSoBwCCJTk5WSdOnAj1GECvkFfYCXlFsHEGFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAA\nRlFAAQAAYBQFFAAAAEZRQDFkeL1eTZw4UV6vN9SjAAGRV9gJeUWwWX4j+sbGRu3atUtFRUWqq6tT\nQkKCFi9erJSUlB63y8nJ0UcffdT1ECNH6vDhw1ZHAdqpr6+X1+tVfX19qEcBAiKvsBPyimCzXEA3\nb96sd999VwsWLFB8fLwKCgq0atUqbd26VZMmTep2u0WLFmnevHnt1urr6/Xv//7vuvXWW61PDgAA\nAFuyVEDLy8tVXFyspUuXKiMjQ5I0d+5cZWVlKS8vT9u3b+92267OkL711luSpO9973tWxgAAAICN\nWfod0KNHjyosLEzp6elta06nU2lpafJ6vTp79qylJz9y5IiioqI0ffp0S9sBV6urq9NjKx9T+gOt\nuUx/IF2PrXxMdXV1IZ4M6Iy8wk7IKwaKpQJ66tQpjRs3TlFRUe3WJ0yYIEmqqKjo9WNdvHhRJSUl\nmjFjhiIiIqyMAbSpq6vTtO9P00u+l+RL9UmSfKk+vVT1kqZ9fxpfJDGokFfYCXnFQLJUQGtra+Vy\nuTqtu1wutbS0qLa2tteP9fbbb6u5uVlz5syxMgLQzppn1qg8qVzNSc3t1psTm1WeVK6nNj4VosmA\nzsgr7IS8YiBZ+h3Qy5cvKzw8vNO60+lsu723jhw5olGjRgV89fzVysvL233u8Xg6nY29ms/nk8/n\n6/b2yMhIJScn9/icgV7153a75Xa7u73d7/d3mrsj9qNVX/Zj74G9ak5tlv4i6fMvF6taPzRf06w9\nB/bowfsf7PExAVO6zOuXH8krBpt2ef3y62pPeR3O34s6Gq77EWifrmapgEZERKixsbHTekNDQ9vt\nveHz+eT1ejV//nyNGNH7k7CLFi1q9/nMmTO1ZMkSZWZmdnn/vLw8bdiwodvHS05O1okTJ3p8zoyM\njB7f92zdunVav359t7efPn06YMkuKyvTxIkTu72d/fhKl/vxcYc7Hfjqnz75LP2QAwy4jnn93Vf/\nJK8YdCzkddh/L7rKcNiP/Px87dixQ++8806Pz9EdSwXU5XKppqam0/qVS+9dXZ7vSlFRkRwOh+VX\nv+/evVsej6ft80A/HWRnZ+uuu+7q9vbIyMiAz7lnz56APx30JCEhQSUlJQHv0xP24ysd9yP9gfS2\n301SlVrL53xJX2tdch91641fvxHwcQETuszrXZLGtC6RVwwm7fJaI2m/eszrcP5e1NFw2I/MzEzd\nfffdnc6AdjxZ2B1LBTQxMVHHjx+X3+9vV/y8Xq8cDoeSkpJ69ThHjhzR2LFj25XJ3vB4PJoyZUqv\n7x/o1HFvBDp9HUhUVJSlmbvCfnyl434suGuBXqp6Sc2JV/2O0tckjZVGVIxQxg8y+j03ECxd5nWM\nyCsGpYHI61D9XmQV+2HxRUipqalqamrSwYMH29YaGxtVWFgoj8ej2NhYSdK5c+dUWVmppqamTo9R\nUVGhyspKXnyEoHj2p8/K8/88GlHRPsojKkbIU+HRxqc2hmgyoDPyCjshrxhIlgqox+NRamqqdu7c\nqby8PL3xxhvKyclRdXW1srOz2+63Y8cOPfTQQ11ern/rrbf6dPkd6Ep0dLTeO/yelo1dJvfvW38S\ndP/erWVjl+m9w+8pOjo6xBMCXyGvsBPyioFk+U9xrl69utPfgt+0aZMmT57cdh+HwyGHw9Fp25aW\nFhUXF2v8+PEaN25c/yYHvhQdHa1t/2ebHrz/QaWkpOiN197gMiYGLfIKOyGvGCiWC2h4eLiys7Pb\nnfHsKDc3V7m5uZ3WHQ6HXn/9datPCQAAgCHE0iV4AAAAoL8ooAAAADCKAgoAAACjKKAAAAAwigIK\nAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMo\noAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAw\nigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAA\nAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAA\nAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoC\nCgAAAKMooAAAADCKAgoAAACjKKAAAAAwaqTVDRobG7Vr1y4VFRWprq5OCQkJWrx4sVJSUnq1fUlJ\niV577TWdPHlSLS0tGjdunDIzMzVz5kyrowAAAMCGLJ8B3bx5s/bt26e5c+dq+fLlCgsL06pVq1RW\nVhZw20OHDmnlypUKDw/Xj3/8Y/3rv/6r/uEf/kFnzpzp0/AAAACwH0tnQMvLy1VcXKylS5cqIyND\nkjR37lxlZWUpLy9P27dv73bbqqoqvfDCC5o/f74effTR/k0NAAAA27J0BvTo0aMKCwtTenp625rT\n6VRaWpq8Xq/Onj3b7bYHDhxQc3OzHn74YUmS3+/v48gAAACwM0tnQE+dOqVx48YpKiqq3fqECRMk\nSRUVFYqNje1y29LSUn3961/X+++/r//4j/9QTU2NoqOj9YMf/EAPP/ywHA5HH3cBAAAAdmKpgNbW\n1srlcnVad7lcamlpUW1tbbfbfvbZZwoLC9Nzzz2nzMxMJSQk6N1339Xu3bvV3NysH/3oR9anBwAA\ngO1YugR/+fJlhYeHd1p3Op1tt3fH7/fr0qVLevjhh/Xggw/q9ttv1+rVqzV16lTt27dv2FySr66W\nZsyQEhNbP/L6Kwxm5BV2Ql5hJ8M9r5YKaEREhBobGzutNzQ0tN3e07aSNHv27Hbrs2fPVkNDgyoq\nKqyMYlv33iv94Q/S6dOtH+fPD/VEQPfIK+yEvMJOhnteLV2Cd7lcqqmp6bR+5dJ7V5fnr/ja176m\nzz//XNdff3279ZiYGLW0tKiuri7g8+fk5GjUqFHt1jIzM5WZmdmb8QcFn6/nz4HBhLzCTsgr7MTu\nec3Pz1d+fn67tYsXL/Z6e0sFNDExUcePH5ff72/3QiSv1yuHw6GkpKRut73pppv0+eefq6amRmPG\njGlbr6mpkcPhUExMTMDn37p1q6ZMmWJl5EHH7W79aefqz4HBirzCTsgr7MTuee3qBGBpaWmv/zCR\npUvwqampampq0sGDB9vWGhsbVVhYKI/H0/YK+HPnzqmyslJNTU1t95s1a5ZaWlr05ptvtq21tLSo\noKBA0dHRGj9+vJVRbGv/fmn6dCkhofXj/v2hngjoHnmFnZBX2Mlwz6ulM6Aej0epqanauXOnzp8/\nr/j4eBUUFKi6ulorV65su9+OHTt0+PBh5efnKy4uTpI0Y8YMTZkyRa+99pouXLigxMREHTt2TCdO\nnNDjjz+ukSMt/1VQW7rhBunYsVBPAfQOeYWdkFfYyXDPq+XWt3r16k5/C37Tpk2aPHly230cDkeX\n7+u5ceNG/fKXv9Q777yjwsJC3XjjjVqzZk2nFyYBAABg6LJcQMPDw5Wdna3s7Oxu75Obm6vc3NxO\n65GRkXr00Uf5U5wAAADDmKXfAQUAAAD6iwIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACj\nKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAA\nMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoA\nAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyig\nAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCK\nAgoAAACjKKAAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigIKAAAA\noyigAAAAMIoCCgAAAKMooAAAADBqpNUNGhsbtWvXLhUVFamurk4JCQlavHixUlJSetyuoKBAzz33\nXKd1h8OhvXv36vrrr7c6CgAAAGzIcgHdvHmz3n33XS1YsEDx8fEqKCjQqlWrtHXrVk2aNKnHbR0O\nhx5++GGNGTOm3fp1111ndQwAAADYlKUCWl5eruLiYi1dulQZGRmSpLlz5yorK0t5eXnavn17wMeY\nOnWqxo8f37dpAQAAYHuWfgf06NGjCgsLU3p6etua0+lUWlqavF6vzp4926vH8fv9am5utjYp0I26\nujo9tvIxpT/Qmsv0B9L12MrHVFdXF+LJgM7IK+yEvGKgWDoDeurUKY0bN05RUVHt1idMmCBJqqio\nUGxsbLfbt7S0KCcnR36/XyNHjtStt96qRx55RPHx8X0YHWj94jjt+9NUnlSu5tRm6WPJl+rTS1Uv\n6e3vv633Dr+n6OjoUI8JSCKvsBfyioFk6QxobW2tXC5Xp3WXy6WWlhbV1tZ2u21kZKTuuOMOrVix\nQs8884wyMzNVWlqq5cuX9/rMKdDRmmfWtH5xTGp/Rr05sVnlSeV6auNTIZoM6Iy8wk7IKwaSpTOg\nly9fVnh4eKd1p9PZdnt3Zs6cqZkzZ7Z9Pn36dN1yyy1asWKFdu/erZycnIDPX15e3u5zj8fT6Wzs\n1Xw+n3w+X7e3R0ZGKjk5ucfn9Hq9qq+v7/Z2t9stt9vd7e1+v7/T3B2xH636sh97D+xt/cn8L5I+\n/3KxqvVD8zXN2nNgjx68/8EeHxMwpV1eP/ty8cvcklcMNu3y+uXX1Z7yOpy/F3U0XPcj0D5dzVIB\njYiIUGNjY6f1hoaGttutmDx5sjwej0pLS3t1/0WLFrX7fObMmVqyZIkyMzO7vH9eXp42bNjQ7eMl\nJyfrxIkTPT5nRkaGvF5vt7evW7dO69ev7/b206dPB3yLqrKyMk2cOLHb29mPr3S5Hx93uNOBr/7p\nky/g3IBRHfP6u6/+SV4x6FjI67D/XnSV4bAf+fn52rFjh955550en6M7lgqoy+VSTU1Np/Url967\nujwfSGxsrD799NNe3Xf37t3yeDxtnwf66SA7O1t33XVXt7dHRkYGfM49e/YE/OmgJwkJCSopKQl4\nn56wH1/puB/pD6TLl/rlT4Cfq/WL412SvnynL/dRt9749RsBHxcwoV1eP5P0pqR5kr78NXjyisGk\nXV6r1PrDfQ95Hc7fizoaDvuRmZmpu+++u9MZ0I4nC7tjqYAmJibq+PHj8vv97Yqf1+uVw+FQUlKS\nlYeT1HoKOSYmplf39Xg8mjJlSq8fO9Cp494IdPo6kKioKEszd4X9+ErH/Vhw1wK9VPWSmhOv+h2l\nMZLGSiMqRijjBxn9nhsIli7zGi/yikFpIPI6VL8XWcV+WHwRUmpqqpqamnTw4MG2tcbGRhUWFsrj\n8bS9Av7cuXOqrKxUU1NT2/0uXrzY6fHef/99nTx5UlOnTu3T8MCzP31Wnv/n0YiK9lEeUTFCngqP\nNj61MUSTAZ2RV9gJecVAslRAPR6PUlNTtXPnTuXl5emNN95QTk6OqqurlZ2d3Xa/HTt26KGHHmp3\nuX7ZsmXasGGDfvOb3+jgwYPasmWLfvrTnyouLk4LFy4M3h5hWImOjtZ7h9/TsrHL5P5960+C7t+7\ntWzsMt4iBIMOeYWdkFcMJMt/inP16tWd/hb8pk2bNHny5Lb7OBwOORyOdtvNmjVL77//vkpKSlRf\nXy+Xy6U777xTP/zhD3t9CR7oSnR0tLb9n2168P4HlZKSojdee4PLmBi0yCvshLxioFguoOHh4crO\nzm53xrOj3Nxc5ebmtlvLyspSVlaW9QkBAAAwpFi6BA8AAAD0FwUUAAAARlFAAQAAYBQFFAAAAEZR\nQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABg\nFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAA\nAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUAB\nAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQF\nFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABG\nUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGCU5QLa2NiovLw8ZWRk6I477tAjjzyikpIS\ny0/8i1/8QrNnz9aaNWssbwsAAAD7slxAN2/erH379mnu3Llavny5wsLCtGrVKpWVlfX6MT7++GMV\nFhYqIiLC6tMDAADA5iwV0PLychUXF+vHP/6xlixZonnz5mnLli2Ki4tTXl5erx9n+/bt+sd//EfF\nxMRYHhgAAAD2ZqmAHj16VGFhYUpPT29bczqdSktLk9fr1dmzZwM+RmFhof785z9r8eLF1qcFAACA\n7VkqoKdOndK4ceMUFRXVbn3ChAmSpIqKih639/v9euWVV7Ro0SJdf/31FkcFAADAUGCpgNbW1srl\ncnVad7lcamlpUW1tbY/b/+pXv1JERIQWLFhgbcohpLpamjFDSkxs/XjmTKgnArpHXmEn5BV2Mtzz\naqmAXr58WeHh4Z3WnU5n2+3d+fTTT7V//34tXbpUI0eOtDjm0HHvvdIf/iCdPt36cf78UE8EdI+8\nwk7IK+xkuOfVUhOMiIhQY2Njp/WGhoa227vz4osvatKkSZoxY4bFEb+Sk5OjUaNGtVvLzMxUZmZm\nnx/TNJ+v58+BwYS8wk7IK+zE7nnNz89Xfn5+u7WLFy/2entLBdTlcqmmpqbT+pVL711dnpek0tJS\nffDBB3r66adVVVXVtt7U1KTLly+rqqpKf/d3f6drrrmmx+ffunWrpkyZYmXkQcftbv1p5+rPgcGK\nvMJOyCvsxO557eoEYGlpqVJSUnq1vaUCmpiYqOPHj8vv97d7IZLX65XD4VBSUlKX2505c0YOh0Nr\n165tt+5wOFRTU6OFCxfqkUce0b333mtlHFvav7/1NLvP1xq2/ftDPRHQPfIKOyGvsJPhnldLBTQ1\nNVWvv/66Dh48qPvuu09S619GKiwslMfjUWxsrCTp3LlzunTpkuLj4xUWFqaUlBQ9/fTTnR7vynuI\n/su//Iu++c1vBmF3Br8bbpCOHQv1FEDvkFfYCXmFnQz3vFoqoB6PR6mpqdq5c6fOnz+v+Ph4FRQU\nqLq6WitXrmy7344dO3T48GHl5+crLi5OsbGxbeX0ai+++KJGjx6t2267rf97AgAAAFuw/HL01atX\na9euXSoqKlJdXZ0SEhK0adMmTZ48ue0+DodDDocj4GP15j4AAAAYWiwX0PDwcGVnZys7O7vb++Tm\n5io3NzfgY/3617+2+vQAAACwOUvvAwoAAAD0FwUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQA\nAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFA\nAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAU\nBRQAAACaLsRIAAAZtUlEQVRGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABG\nUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAA\nYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQA\nAABGUUABAABgFAUUAAAARlFAAQAAYBQFFAAAAEaNtLpBY2Ojdu3apaKiItXV1SkhIUGLFy9WSkpK\nj9v96U9/0m9/+1tVVFTowoULuvbaa/WNb3xD999/v7773e/2eQcAAABgL5bPgG7evFn79u3T3Llz\ntXz5coWFhWnVqlUqKyvrcbtPP/1UYWFhuuuuu7RixQr98z//sy5duqQnn3xSRUVFfd4BAAAA2Iul\nM6Dl5eUqLi7W0qVLlZGRIUmaO3eusrKylJeXp+3bt3e77bx58zRv3rx2a3fffbcyMzP1xhtvaM6c\nOX0YHwAAAHZj6Qzo0aNHFRYWpvT09LY1p9OptLQ0eb1enT171tKTO51OxcTEKCwszNJ2AAAAsC9L\nBfTUqVMaN26coqKi2q1PmDBBklRRURHwMf7617/q4sWLqqys1CuvvKJPP/1U999/v5UxgHbq6ur0\n2MrHlP5A6w9G6Q+k67GVj6muri7EkwGdkVfYCXnFQLF0Cb62tlYul6vTusvlUktLi2prawM+xoYN\nG/TBBx9IkqKiorR27VpNnTrVyhhAm7q6Ok37/jSVJ5WrObVZ+ljypfr0UtVLevv7b+u9w+8pOjo6\n1GMCksgr7IW8YiBZOgN6+fJlhYeHd1p3Op1ttweyZMkS/fznP9fKlSv193//93rmmWdUUlJiZQyg\nzZpn1rR+cUxqbrfenNis8qRyPbXxqRBNBnRGXmEn5BUDydIZ0IiICDU2NnZab2hoaLs9kMTExLZ/\nz5kzR0uWLNG2bdv06quvBty2vLy83ecej6fTrwNczefzyefzdXt7ZGSkkpOTe3xOr9er+vr6bm93\nu91yu93d3u73+zvN3RH70aov+7H3wN7Wn8z/Iqnqy8Wa1g/N1zRrz4E9evD+B3t8TMCULvP65Ufy\nisGmXV6//LraU16H8/eijobrfgTap6tZKqAul0s1NTWd1q9ceu/q8nyPTz5ypG677Tbl5+fr0qVL\nuu6663q8/6JFi9p9PnPmTC1ZskSZmZld3j8vL08bNmzo9vGSk5N14sSJHp8zIyNDXq+329vXrVun\n9evXd3v76dOnA75HallZmSZOnNjt7ezHV7rcj4873Gn/V//0yRdwbsCojnk98NU/ySsGHQt5Hfbf\ni64yHPYjPz9fO3bs0DvvvNPjc3THUgFNTEzU8ePH5ff727Vyr9crh8OhpKQkywNcuWzvcDgC3nf3\n7t3yeDxtnwf66SA7O1t33XVXt7dHRkYGfM49e/YE/OmgJwkJCQF/xSAhIaHH29mPr3Tcj/QH0uVL\n/fInwM8l/U7SXZLGtC65j7r1xq/fCPi4gAld5nWepPjWJfKKwaRdXqvUWj57yOtw/l7U0XDYj8zM\nTN19992dzoB2PFnYHUsFNDU1Va+//roOHjyo++67T1LrX0YqLCyUx+NRbGysJOncuXO6dOmS4uPj\n295i6cKFC4qJiWn3eJcuXdK7776rb37zm7r22msDPr/H49GUKVN6PW+gU8e9Eej0dSBRUVGWZu4K\n+/GVjvux4K4FeqnqJTUnXvU7SmMkjZVGVIxQxg8y+j03ECxd5jVe5BWD0kDkdah+L7KK/bBYQD0e\nj1JTU7Vz506dP39e8fHxKigoUHV1tVauXNl2vx07dujw4cPKz89XXFycJCk3N1exsbHyeDyKiYlR\ndXW1CgsLdf78eeXm5vZpeODZnz6rt7//tspbytV8zVdfJEdUjJCnwqON/3djCKcD2iOvsBPyioFk\n+U9xrl69Wvfee6+Kior04osvqrm5WZs2bdLkyZPb7uNwODpdUk9LS1NdXZ327dunbdu26Xe/+508\nHo9eeukl3Xzzzf3fEwxL0dHReu/we1o2dpncv2/9SdD9e7eWjV3GW4Rg0CGvsBPyioHkKC4ubgn1\nEIGcPHlS2dnZKikp4fIUulVaWqqUlBRyAlsgr7AT8oreuJKTvLw8jR8/vsf7Wj4DCgAAAPQHBRQA\nAABGUUAxZFx5093evPUEEGrkFXZCXhFsll4FDwxmvXnTXWCwIK+wE/KKYOMMKAAAAIyigAIAAMAo\nCigAAACMooACAADAKAooAAAAjKKAAgAAwCgKKAAAAIyigGLI8Hq9mjhxorxeb6hHAQIir7AT8opg\no4BiyKivr5fX61V9fX2oRwECIq+wE/KKYKOAAgAAwCgKKAAAAIyigAIAAMAoCigAAACMooACAADA\nKAooAAAAjKKAYshwu91at26d3G53qEcBAiKvsBPyimAbGeoBgGBxu91av359qMcAeoW8wk7IK4KN\nM6AAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAYMvx+v06cOCG/3x/qUYCA\nyCvshLwi2CigGDLKy8s1adIklZeXh3oUICDyCjshrwg2CigAAACMooACAADAKAooAAAAjKKAAgAA\nwCgKKAAAAIyigAIAAMAoCigAAACMGhnqAYBg8Xg8KisrU0JCQqhHAQIir7AT8opgo4BiyIiKitLE\niRNDPQbQK+QVdkJeEWxcggcAAIBRFFAAAAAYRQEFAACAURRQAAAAGEUBBQAAgFEUUAAAABhFAcWQ\n4fP5tH79evl8vlCPAgREXmEn5BXBRgHFkOHz+bRhwwa+QMIWyCvshLwi2CigAAAAMIoCCgAAAKMo\noAAAADCKAgoAAACjKKAAAAAwigIKAAAAoyigGDIiIyOVnJysyMjIUI8CBEReYSfkFcE20uoGjY2N\n2rVrl4qKilRXV6eEhAQtXrxYKSkpPW5XWlqqoqIilZWV6ezZsxo9erRuvvlmZWVlafTo0X3eAeCK\n5ORknThxItRjAL1CXmEn5BXBZvkM6ObNm7Vv3z7NnTtXy5cvV1hYmFatWqWysrIet9uxY4c++ugj\n3X777Vq+fLlmz56td955R0uWLNH58+f7vAMAAACwF0tnQMvLy1VcXKylS5cqIyNDkjR37lxlZWUp\nLy9P27dv73bbRx99VJMnT263duutt2rFihX6z//8T2VlZfVhfAAAANiNpTOgR48eVVhYmNLT09vW\nnE6n0tLS5PV6dfbs2W637Vg+Jenb3/62oqOjVVlZaWUMAAAA2JilAnrq1CmNGzdOUVFR7dYnTJgg\nSaqoqLD05H6/X36/X6NGjbK0HQAAAOzLUgGtra2Vy+XqtO5yudTS0qLa2lpLT7537141NTVp1qxZ\nlrazs+pqacYMKTGx9eOZM6GeCOgeeYWdkFfYyXDPq6UCevnyZYWHh3dadzqdbbf31kcffaRXX31V\nM2fO1He+8x0rY9javfdKf/iDdPp068f580M9EdA98go7Ia+wk+GeV0svQoqIiFBjY2On9YaGhrbb\ne6OyslJr165VQkKCfvKTn/T6+XNycjpdrs/MzFRmZmavHyPUfL6ePwcGE/IKOyGvsBO75zU/P1/5\n+fnt1i5evNjr7S0VUJfLpZqamk7rVy69d3V5vqMzZ87oiSeeUHR0tDZt2tTp90l7snXrVk2ZMqX3\nAw9CbnfrTztXf47g8Hq9ysjI0J49e5ScnBzqcYYE8jpwyGvwkdeBQ16Dz+557eoEYGlpacD3hb/C\nUgFNTEzU8ePH5ff72xVHr9crh8OhpKSkHrf/4osv9MQTT6ipqUlbt24dlm9Av39/62l2n681bPv3\nh3qioaO+vl5er1f19fWhHmXIIK8Dh7wGH3kdOOQ1+IZ7Xi0V0NTUVL3++us6ePCg7rvvPkmtfxmp\nsLBQHo9HsbGxkqRz587p0qVLio+PV1hYmKTW8Obm5qq2tlbPP/+8xo4dG+RdsYcbbpCOHQv1FEDv\nkFfYCXmFnQz3vFoqoB6PR6mpqdq5c6fOnz+v+Ph4FRQUqLq6WitXrmy7344dO3T48GHl5+crLi5O\nkrRx40Z9/PHHSktL0yeffKJPPvmk7f5RUVGaMWNGkHYJAAAAg5nlvwW/evXqTn8LftOmTe3eaN7h\ncMjhcLTb7tSpU3I4HDp06JAOHTrU7ra4uDgKKAAAwDBhuYCGh4crOztb2dnZ3d4nNzdXubm57dY6\nvlIKAAAAw5Ol9wEFAAAA+osCCgAAAKMooBgy3G631q1bJ7fd3kwNwxJ5hZ2QVwSb5d8BBQYrt9ut\n9evXh3oMoFfIK+yEvCLYOAMKAAAAoyigAAAAMIoCCgAAAKMooAAAADCKAgoAAACjKKAAAAAwigKK\nIcPv9+vEiRPy+/2hHgUIiLzCTsgrgo0CiiGjvLxckyZNUnl5eahHAQIir7AT8opgo4ACAADAKAoo\nAAAAjKKAAgAAwCgKKAAAAIyigAIAAMAoCigAAACMooACAADAqJGhHgAIFo/Ho7KyMiUkJIR6FCAg\n8go7Ia8INgoohoyoqChNnDgx1GMAvUJeYSfkFcHGJXgAAAAYRQEFAACAURRQAAAAGEUBBQAAgFEU\nUAAAABhFAQUAAIBRFNAQyc/PD/UIQ47P59O9994rn88X6lGGHPIafOR14JDX4COvA2e45pUCGiLD\nNXADyefzaf/+/XyBHADkNfjI68Ahr8FHXgfOcM0rBRQAAABGUUABAABgFAUUAAAARtnib8E3NDRI\nksrLy0M8SfBcvHhRpaWloR5jSLmSj6GUk8GCvAYfeR045DX4yOvAGUp5vZKPK72tJ47i4uKWgR6o\nv4qKivTss8+GegwAAAAEsGbNGs2ZM6fH+9iigF64cEEffvihxowZI6fTGepxAAAA0EFDQ4Oqqqp0\nyy23KCYmpsf72qKAAgAAYOjgRUgAAAAwigIKAAAAoyigAAAAMIoCCgAAAKMooAAAADDKFm9EbxeN\njY3atWuXioqKVFdXp4SEBC1evFgpKSk9bldaWqqioiKVlZXp7NmzGj16tG6++WZlZWVp9OjRhqYf\nvPp6XP/0pz/pt7/9rSoqKnThwgVde+21+sY3vqH7779f3/3udw1NP3j19bh29Itf/EJvvvmmpk2b\nxvv1qu/HtaCgQM8991yndYfDob179+r6668fqJFtob95LSkp0WuvvaaTJ0+qpaVF48aNU2ZmpmbO\nnDmwgw9yfT2uOTk5+uijj7q8beTIkTp8+PBAjGsb/clrSUmJfv3rX+v06dO6fPmy3G635s2bp7vv\nvlsjRgyd84ZhDz300PpQDzFU/OxnP1NBQYHuvPNOzZkzR6dPn9ZvfvMbTZkyRTfccEO3223YsEGV\nlZVKTU3VnDlzFBsbq4KCAh06dEhz585VVFSUwb0YfPp6XD/88EP95S9/0fTp0zVr1iyNHz9e//M/\n/6M9e/YoPj5eCQkJBvdi8Onrcb3axx9/rO3bt8vpdMrtdut73/veAE89+PX1uFZUVOi//uu/lJWV\npX/6p3/S7bff3vbf+PHjFRYWZnAvBp/+5PXQoUPasGGDbrzxRt1999269dZbFRUVpZaWFk2cONHQ\nHgxOfT2usbGxuvXWW9vl9NZbb9V7772n7373u8P+a0Ffj+sf//hHrVq1SpGRkZo/f76mTZumuro6\n7d27V//7v/+rqVOnGtyLgcX7gAZJeXm5Hn30US1dulQZGRmSWt+QNSsrS9dff722b9/e7bb//d//\nrcmTJ7db+9Of/qQVK1Zo0aJFysrKGtDZB7P+HNeuNDQ0KDMzUzfeeKOef/75gRjZFoJ1XJctW6Zv\nfOMbKikpUUJCwrA/A9qf41pQUKCf//znevnllzV+/HhTI9tCf45rVVWVHn74YaWnp+vRRx81NbIt\nBPvr61tvvaVNmzbpqaee0uzZswdiZFvoz3F99tln9fvf/1779u3Tdddd17a+YsUKnTp1SgcPHhzw\n+U0ZOudyQ+zo0aMKCwtTenp625rT6VRaWpq8Xq/Onj3b7bYdy6ckffvb31Z0dLQqKysHZF676M9x\n7YrT6VRMTMywP5sUjONaWFioP//5z1q8ePFAjmorwcqr3+9Xc3PzQI1pO/05rgcOHFBzc7Mefvhh\nSa3HFq2C/fX1yJEjioqK0vTp04M9qq3057g6nU45nc525VOSRo8erYiIiAGbORQooEFy6tQpjRs3\nrtPl8gkTJkhqvbxmhd/vl9/v16hRo4I2ox0F47j+9a9/1cWLF1VZWalXXnlFn376qe6///4Bmdcu\n+ntc/X6/XnnlFS1atGjY/27i1fp7XFtaWpSTk6N58+bpjjvu0Jo1a/T5558P2Lx20Z/jWlpaqq9/\n/et6//33dd9992nevHn6wQ9+oF27dqmlZXhfAAzm962LFy+qpKREM2bMGHJFyar+HNd77rlHLS0t\n2rJliyorK1VdXa0DBw7o2LFjWrhw4YDObRovQgqS2tpauVyuTusul0stLS2qra219Hh79+5VU1OT\nZs2aFawRbSkYx3XDhg364IMPJElRUVFau3btkPo9mr7o73H91a9+pYiICC1YsGCgRrSl/hzXyMhI\n3XHHHbr55pt1zTXX6OTJk3r99de1fPly5eXlKTY2diBHH9T6c1w/++wzhYWF6bnnnlNmZqYSEhL0\n7rvvavfu3WpubtaPfvSjgRx9UAvm9623335bzc3NmjNnTjBHtKX+HNekpCRt2bJFa9as0e9+9ztJ\nUlhYmB577DHdeeedAzZzKFBAg+Ty5csKDw/vtO50Ottu762PPvpIr776qmbOnKnvfOc7QZvRjoJx\nXJcsWaL77rtPZ8+e1YEDB/TMM8/oZz/7meVXew8l/Tmun376qfbv36+1a9dq5Ei+hFytP8d15syZ\n7V6RPX36dN1yyy1asWKFdu/erZycnKDPaxf9Oa5XLrkvWbKk7crH7bffri+++EL79u3TwoULh+0L\nPYP5fevIkSMaNWrUsP66ekV/jmtlZaWefPJJxcXFKTs7W06nU0eOHNELL7yg0aNHD6lfb+ASfJBE\nRESosbGx03pDQ0Pb7b1RWVmptWvXKiEhQT/5yU+COqMdBeO4JiYmKiUlRXfccYdeeOEFxcfHa9u2\nbUGf1U76c1xffPFFTZo0STNmzBiw+ewqWF8Hrpg8ebI8Ho9KS0uDMp9d9ee4Xrmt44tiZs+erYaG\nBsu/HjWUBCuvPp9PXq9Xs2fPHlJvE9RX/TmuL7/8ssLCwvT8889r7ty5Sk1N1dNPP63Jkydr27Zt\nQ+p3w0lKkLhcri5Pq19Z6+p0fEdnzpzRE088oejoaG3atGnY/lR+tWAc16uNHDlSt912mz7//HNd\nunQpKDPaUV+Pa2lpqT744APNnz9fVVVVbf81NTXp8uXLqqqq0l//+tcBnX0wC3Zepda3u/niiy/6\nPZud9ee4fu1rX5OkTr+rHBMTo5aWFtXV1QVxUnsJVl6LiorkcDiG/VsvXdGf41pWVqabb765U0m9\n7bbbVFtbq6qqquAOG0IU0CBJTEzUZ5991ukVll6vVw6HQ0lJST1u/8UXX+iJJ55QU1OTnnvuOd6A\n/kv9Pa5duXL5w+FwBGVGO+rrcT1z5owcDofWrl2rBx54oO2/2tpalZaWauHChTp06JCJXRiUBiKv\nPp9PMTExwRrRlvpzXG+66SZJUk1NTbv1mpoaORyOYX1sg5XXI0eOaOzYsfJ4PAMxpu3057j+7W9/\nU1NTU5frkrq8za4ooEGSmpqqpqamdu/R1djYqMLCQnk8nrYXEJw7d06VlZXtQlRfX6/c3FzV1tZq\n8+bNGjt2rPH5B6v+HNcLFy50erxLly7p3Xff1Te/+U1de+21A78Dg1Rfj2tKSoqefvppPfPMM+3+\nGzVqlL71rW/pmWee0W233RaSfRoM+pPXixcvdnq8999/XydPnhz2L5rrz3GdNWuWWlpa9Oabb7at\ntbS0qKCgQNHR0cP6PVf7c1yvqKioUGVlJS8+ukp/jutNN92kkpKSdmfmm5ubVVxcrKioqCHVD3gF\nQZB4PB6lpqZq586dOn/+vOLj41VQUKDq6mqtXLmy7X47duzQ4cOHlZ+fr7i4OEnSxo0b9fHHHyst\nLU2ffPKJPvnkk7b7R0VFDevftevPcc3NzVVsbKw8Ho9iYmJUXV2twsJCnT9/Xrm5uaHapUGhr8c1\nNja2y1djv/jiixo9evSwLp9S//K6bNkyJSUl6Vvf+pauvfZanTx5UgUFBYqLixtyb79iVX+O64wZ\nMzRlyhS99tprunDhghITE3Xs2DGdOHFCjz/++LB+IV1/jusVb731FpffO+jPcV20aJGefPJJLV26\nVOnp6YqIiNCRI0dUUVGhxYsXD6n3sB6+/+cNgNWrV3f626+bNm1q90bzDoej06XfU6dOyeFw6NCh\nQ50uX8bFxQ3rAir1/bimpaXp7bff1r59+3Tp0iVFR0dr8uTJWrhwYdtlueGsr8e1K8P51xk66utx\nnTVrlt5//32VlJSovr5eLpdLd955p374wx8O68vEV/Qnrxs3btQvf/lLvfPOOyosLNSNN96oNWvW\nDOu/1nNFf45rS0uLiouLNX78eI0bN87k2INeX4/r1KlTtXnzZu3evVuvvvqqmpqadOONN+rxxx/X\nvHnzTO/GgOJPcQIAAMAofgcUAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABgFAUU\nAAAARlFAAQAAYBQFFAAAAEZRQAEAAGAUBRQAAABGUUABAABg1P8HNyOw9OYX1aEAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6dd4b80c90>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Exact 2D Voronoi diagrams\n",
    "#\n",
    "#  Given a set of points in 2D, we can compute the exact 2D \n",
    "#  Voronoi diagram.  A plot of the diagram shows the boundaries\n",
    "#  defining each Voronoi polygon.  Some of the lines seem to\n",
    "#  go off to infinity, and they really do.\n",
    "#\n",
    "#  We saw in the 2D geometry notebook that we could define a set of\n",
    "#  points, and then compute the Voronoi diagram using the \n",
    "#  vor=scipy.spatial.Voronoi() command, followed by the \n",
    "#  spatial.voronoi_plot_2d(vor) command to actually see the diagram.\n",
    "#\n",
    "#  If we use a very regular set of points, we may be able to guess what\n",
    "#  the Voronoi diagram looks like.  So make a tensor product of the points\n",
    "#  we called \"g\" for the 1d problem.\n",
    "#\n",
    "#  We do this using the numpy.meshgrid() command:\n",
    "#    X, Y = numpy.meshgrid ( g, g )\n",
    "#\n",
    "g = np.array ( [ 1.0/5.0, 2.0/5.0, 3.0/5.0, 4.0/5.0] )\n",
    "\n",
    "X, Y = np.meshgrid ( g, g )\n",
    "#\n",
    "#  Unfortunately, the Voronoi command doesn't want two matrices of X and Y coordinates.\n",
    "#  It wants a single Nx2 array, with each row containing an (X,Y) pair.\n",
    "#  Here's one way to make Voronoi happy:\n",
    "#\n",
    "XY = np.zeros ( [ 16, 2 ] )\n",
    "XY[:,0] = np.ndarray.flatten ( X )\n",
    "XY[:,1] = np.ndarray.flatten ( Y )\n",
    "#\n",
    "#  Now compute the Voronoi diagram.\n",
    "#\n",
    "vor = spatial.Voronoi ( XY )\n",
    "#\n",
    "#  Now plot the Voronoi diagram. Note that the plotting command arbitrarily chooses \n",
    "#  the plot range, so it shows the nonunit square [0.1,0.9]x[0.1,0.9], which may give \n",
    "#  the false impression that the Voronoi polygons all have the same size, and\n",
    "#  the generators g are the centroids of their region.  However, the\n",
    "#  layer of Voronoi subregions along the boundary are not the same shape\n",
    "#  as those away from the boundary, and their generators are not at the\n",
    "#  centers of their subregions!  We will demonstrate this in a minute.\n",
    "#\n",
    "spatial.voronoi_plot_2d ( vor )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 862.  612.  599.  886.  595.  421.  422.  599.  602.  404.  424.  583.\n",
      "  878.  578.  628.  907.]\n",
      "\n",
      " I            G[I,0]     G[I,1]          Centroid[I,0]  Centroid[I,1]       Area[I]       Energy[I]\n",
      " 0        0.200000        0.200000        0.157865        0.154298        0.086200        0.001601\n",
      " 1        0.400000        0.200000        0.395086        0.153764        0.061200        0.000765\n",
      " 2        0.600000        0.200000        0.601660        0.148880        0.059900        0.000817\n",
      " 3        0.800000        0.200000        0.849418        0.147855        0.088600        0.001764\n",
      " 4        0.200000        0.400000        0.149250        0.401418        0.059500        0.000767\n",
      " 5        0.400000        0.400000        0.396536        0.397369        0.042100        0.000273\n",
      " 6        0.600000        0.400000        0.601648        0.402212        0.042200        0.000290\n",
      " 7        0.800000        0.400000        0.853284        0.404502        0.059900        0.000823\n",
      " 8        0.200000        0.600000        0.151824        0.595748        0.060200        0.000787\n",
      " 9        0.400000        0.600000        0.397330        0.601619        0.040400        0.000270\n",
      "10        0.600000        0.600000        0.600903        0.605820        0.042400        0.000292\n",
      "11        0.800000        0.600000        0.850759        0.603529        0.058300        0.000813\n",
      "12        0.200000        0.800000        0.150795        0.846389        0.087800        0.001723\n",
      "13        0.400000        0.800000        0.401234        0.852154        0.057800        0.000795\n",
      "14        0.600000        0.800000        0.600214        0.850952        0.062800        0.000847\n",
      "15        0.800000        0.800000        0.849170        0.849427        0.090700        0.001789\n"
     ]
    }
   ],
   "source": [
    "## Approximate 2D Voronoi diagrams by sampling\n",
    "#\n",
    "#  Let's consider our 2D problem, and try to analyze it using sampling, to get\n",
    "#  the area, centroid, and energy of each Voronoi subregion.\n",
    "#\n",
    "#  We are restricting our attention to the unit square, so we can easily compute\n",
    "#  a 10000x2 array of random values by calling np.random.rand ( 10000, 2 )\n",
    "#\n",
    "g = np.array ( [ 1.0/5.0, 2.0/5.0, 3.0/5.0, 4.0/5.0] )\n",
    "\n",
    "GX, GY = np.meshgrid ( g, g )\n",
    "G = np.zeros ( [ 16, 2 ] )\n",
    "G[:,0] = np.ndarray.flatten ( GX )\n",
    "G[:,1] = np.ndarray.flatten ( GY )\n",
    "\n",
    "centroid = np.zeros ( [ 16, 2 ] )\n",
    "area = np.zeros ( 16 )\n",
    "energy = np.zeros ( 16 )\n",
    "near = np.zeros ( 16 )\n",
    "\n",
    "n = 10000\n",
    "XY = np.random.rand ( n, 2 )\n",
    "\n",
    "def nearest ( G, xy ):\n",
    "  idx = -1\n",
    "  dmin = np.Inf\n",
    "  for i in range ( 0, 16 ):\n",
    "    d = np.linalg.norm ( G[i,:] - xy[:] )    \n",
    "    if ( d < dmin ):\n",
    "      idx = i\n",
    "      dmin = d\n",
    "  return idx  \n",
    "\n",
    "for j in range ( 0, n ):\n",
    "  i = nearest ( G, XY[j,:] )\n",
    "  near[i] = near[i] + 1\n",
    "  centroid[i,:] = centroid[i,:] + XY[j]\n",
    "  energy[i] = energy[i] + ( G[i,0] - XY[j,0] ) ** 2 + ( G[i,1] - XY[j,1] ) ** 2\n",
    "\n",
    "print ( near )\n",
    "    \n",
    "for i in range ( 0, 16 ):\n",
    "  area[i] = 1.0 * float ( near[i] ) / float ( n )\n",
    "  centroid[i,:] = centroid[i,:] / float ( near[i] )\n",
    "  energy[i] = area[i] * energy[i] / float ( near[i] )\n",
    "# \n",
    "#  Print i, g[i], centroid[i,:], area[i], energy[i]\n",
    "#\n",
    "print ( '' )\n",
    "print ( ' I            G[I,0]     G[I,1]          Centroid[I,0]  Centroid[I,1]       Area[I]       Energy[I]' )\n",
    "for i in range ( 0, 16 ):\n",
    "  print ( '%2d  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f  %14.6f' \\\n",
    "  % ( i, G[i,0], G[i,1], centroid[i,0], centroid[i,1], area[i], energy[i] ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comments##\n",
    "\n",
    "There are four generators in corners, such as (0.2,0.2).  \n",
    "Our results suggest that these generators have energies around 0.0017,\n",
    "and areas around 0.09.\n",
    "\n",
    "There are 8 generators away from corners, but near boundaries.\n",
    "Our results suggest all these generators have energies around 0.0007,\n",
    "and areas around 0.06.\n",
    "\n",
    "Finally, there are four generators tucked in the middle of the region,\n",
    "away from the boundaries.  They have energies around 0.00028 (the least),\n",
    "and areas around 0.04.\n",
    "\n",
    "Notice that the generators near boundaries are not close to the centroids\n",
    "we compute, but the four interior generators are.\n",
    "\n",
    "Our goal, later on, will be to automatically adjust the points G so that\n",
    "each subregion has about the same area, and is very close to the centroid...\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}