diff --git a/analytics.py b/analytics.py new file mode 100644 index 0000000..270361a --- /dev/null +++ b/analytics.py @@ -0,0 +1,156 @@ +import math + +from point import Point + +def mean_center(points): + """ + Given a set of points, compute the mean center + + Parameters + ---------- + points : list + A list of points in the form (x,y) + + Returns + ------- + x : float + Mean x coordinate + + y : float + Mean y coordinate + """ + x = None + y = None + sum_x=[] + sum_y=[] + for x_tmp,y_tmp in points: + sum_x.append(x_tmp) + sum_y.append(y_tmp) + + x=float(sum(sum_x)/len(sum_x)) + y=float(sum(sum_y)/len(sum_y)) + + return x, y + +def euclidean_distance(a, b): + """ + Compute the Euclidean distance between two points + + Parameters + ---------- + a : tuple + A point in the form (x,y) + + b : tuple + A point in the form (x,y) + + Returns + ------- + + distance : float + The Euclidean distance between the two points + """ + distance = math.sqrt((a[0] - b[0])**2 + (a[1] - b[1])**2) + return distance + +def average_nearest_neighbor_distance(points,mark=None): + """ + Given a set of points, compute the average nearest neighbor. + + Parameters + ---------- + points : list + A list of points + mark : str + + Returns + ------- + mean_d : float + Average nearest neighbor distance + + References + ---------- + Clark and Evan (1954 Distance to Nearest Neighbor as a + Measure of Spatial Relationships in Populations. Ecology. 35(4) + p. 445-453. + """ + mean_d = 0 + if mark==None: + points_tmp=points + else: + points_tmp=[ point for point in points if point.mark==mark] + + length=len(points_tmp) + nearest_distances=[] + for i in range(length): + distance=[] + for j in range(length): + if i==j: + continue + else: + distance.append(euclidean_distance((points_tmp[i].x,points_tmp[i].y),(points_tmp[j].x,points_tmp[j].y))) + nearest_distances.append(min(distance)) + + mean_d=float(sum(nearest_distances)/len(nearest_distances)) + return mean_d + + +def minimum_bounding_rectangle(points): + """ + Given a set of points, compute the minimum bounding rectangle. + + Parameters + ---------- + points : list + A list of points in the form (x,y) + + Returns + ------- + : list + Corners of the MBR in the form [xmin, ymin, xmax, ymax] + """ + + mbr = [0,0,0,0] + + x_list=[] + y_list=[] + for x,y in points: + x_list.append(x) + y_list.append(y) + mbr=[min(x_list),min(y_list),max(x_list),max(y_list)] + + return mbr + + +def mbr_area(mbr): + """ + Compute the area of a minimum bounding rectangle + """ + area = 0 + area=(mbr[2]-mbr[0])*(mbr[3]-mbr[1]) + return area + + +def expected_distance(area, n): + """ + Compute the expected mean distance given + some study area. + + This makes lots of assumptions and is not + necessarily how you would want to compute + this. This is just an example of the full + analysis pipe, e.g. compute the mean distance + and the expected mean distance. + + Parameters + ---------- + area : float + The area of the study area + + n : int + The number of points + """ + + expected = 0 + expected =float((math.sqrt(area/n))/2) + return expected diff --git a/point.py b/point.py new file mode 100644 index 0000000..829745b --- /dev/null +++ b/point.py @@ -0,0 +1,28 @@ + + + +class Point(): + + def __init__(self,x,y,mark=None): + self.x=x + self.y=y + self.mark=mark + + def check_coincident(self, peer_p): + + return (self.x == peer_p.x and self.y == peer_p.y and self.mark == peer_p.mark) + + def shift_point(self, x_shift, y_shift): + + self.x += x_shift + self.y += y_shift + + def __eq__(self, other): + return self.x == other.x and self.y == other.y and self.mark == other.mark + + def __str__(self): + return "x=%f,y=%f,mark=%s"%(self.x,self.y,self.mark) + + def __add__(self, other): + return Point(self.x+other.x,self.y+other.y,self.mark) + diff --git a/point_pattern.ipynb b/point_pattern.ipynb new file mode 100644 index 0000000..53d0182 --- /dev/null +++ b/point_pattern.ipynb @@ -0,0 +1,281 @@ +{ + "metadata": { + "name": "point_pattern" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from point_pattern import PointPattern\n", + "from point import Point\n", + "from utils import *\n", + "\n", + "import pysal as ps\n", + "\n", + "shapefile = ps.open(ps.examples.get_path('new_haven_merged.shp'))\n", + "dbf = ps.open(ps.examples.get_path('new_haven_merged.dbf'))\n", + "\n", + "point_list=[]\n", + "\n", + "for geometry, attributes in zip(shapefile, dbf):\n", + " point_list.append(Point(geometry[0],geometry[1],attributes[0]))\n", + "\n", + "ptp=PointPattern(point_list)\n", + "\n", + "# Run a few tests to explore the dataset.\n", + "nn = ptp.average_nearest_neighbor_distance_KDTree()\n", + "print('This dataset has a nearest neighbor distance of {}'.format(nn))\n", + "# Use your monte carlo simulation code to see if the result is statistically significant\n", + "\n", + "smallest, largest = ptp.critical_points(ptp.permutations())\n", + "\n", + "if is_observed_distance_significant(smallest, largest, nn):\n", + " print('The mark is significant.')\n", + "else:\n", + " print('The mark is not significant.')\n", + " \n", + "nn = ptp.average_nearest_neighbor_distance_KDTree(' Thu, Sept. 18th 2014')\n", + "print('This Thu, Sept. 18th 2014 has a nearest neighbor distance of {}'.format(nn))\n", + "# Use your monte carlo simulation code to see if the result is statistically significant\n", + "if is_observed_distance_significant(smallest, largest, nn):\n", + " print('The mark is significant.')\n", + "else:\n", + " print('The mark is not significant.')\n", + "\n", + " \n", + "nn = ptp.average_nearest_neighbor_distance_KDTree(' Wed, Sept. 3rd 2014')\n", + "print('This Wed, Sept. 3rd 2014 has a nearest neighbor distance of {}'.format(nn))\n", + "# Use your monte carlo simulation code to see if the result is statistically significant\n", + "if is_observed_distance_significant(smallest, largest, nn):\n", + " print('The mark is significant.')\n", + "else:\n", + " print('The mark is not significant.')" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "This dataset has a nearest neighbor distance of 0.000240361352928\n", + "The mark is not significant." + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "This Thu, Sept. 18th 2014 has a nearest neighbor distance of 0.00259760178027\n", + "The mark is significant.\n", + "This Wed, Sept. 3rd 2014 has a nearest neighbor distance of 0.00189031697659" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "The mark is significant.\n" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "compute_g = ptp.compute_g(10)\n", + "print(compute_g)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "0.111111111111\n" + ] + } + ], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "new_point_list=[ point for point in point_list if point.mark==' Thu, Sept. 18th 2014']\n", + "new_ptp=PointPattern(new_point_list)\n", + "compute_g = new_ptp.compute_g(10)\n", + "print(compute_g)\n", + "\n", + "new_point_list=[ point for point in point_list if point.mark==' Wed, Sept. 3rd 2014']\n", + "new_ptp=PointPattern(new_point_list)\n", + "compute_g = new_ptp.compute_g(10)\n", + "print(compute_g)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "0.111111111111\n", + "0.111111111111\n" + ] + } + ], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "%pylab inline\n", + "import matplotlib.pyplot as plt\n", + "from point_pattern import PointPattern\n", + "from point import Point\n", + "from utils import *\n", + "\n", + "import pysal as ps\n", + "\n", + "shapefile = ps.open(ps.examples.get_path('new_haven_merged.shp'))\n", + "dbf = ps.open(ps.examples.get_path('new_haven_merged.dbf'))\n", + "\n", + "point_list=[]\n", + "\n", + "for geometry, attributes in zip(shapefile, dbf):\n", + " point_list.append(Point(geometry[0],geometry[1],attributes[0]))\n", + "\n", + "ptp=PointPattern(point_list)\n", + "\n", + "def plot_fun(x,y):\n", + " \n", + " plt.plot(x, y, 'g+')\n", + "\n", + "x=[point.x for point in ptp.points]\n", + "y=[point.y for point in ptp.points]\n", + "\n", + "plot_fun(x,y)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].\n", + "For more information, type 'help(pylab)'.\n" + ] + }, + { + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAD9CAYAAACrxZCnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX14VOWd8P+ZlySTF8gQAkSIZCioaYpksL7QLphB60JR\nd0PXZ9e2rubx0qvah7W4P9vari0T6mKqtI0u7bJPrx87uOz+6FNaZi1NYdnKRCJOheqE8tAIic5A\njBMMyUxIwrzP74/T++bMZCaZQHjT87kuLiYzZ865z5lz7u/9fdclk8kkGhoaGhoaOaK/3APQ0NDQ\n0Li60ASHhoaGhsaE0ASHhoaGhsaE0ASHhoaGhsaE0ASHhoaGhsaE0ASHhoaGhsaEyElwxONxFi9e\nzL333gvAz3/+cz71qU9hMBh46623Mn4nFApx2223YbVaqamp4Vvf+lbK5//0T//EJz/5SRYuXMg3\nv/nNCzwNDQ0NDY1LhTGXjV588UVqamo4c+YMADfeeCM7d+7kK1/5StbvmEwm9u3bR1FREbFYjKVL\nl9LW1sbSpUvZt28fr7zyCocPHyYvL48PP/xwcs5GQ0NDQ+OiM67G0d3dTUtLC4888ggiV7C6uprr\nr79+3J0XFRUBEIlEiMfjlJWVAfDP//zPfOtb3yIvLw+AGTNmnPcJaGhoaGhcWsbVOJ588kleeOEF\nBgcHJ7zzRCLBTTfdRFdXF48//jg1NTUAHD9+nNdee41vf/vbmEwmNm7cyM0335zyXZ1ON+HjaWho\naGjAxS4IMqbGsWvXLmbOnMnixYvPayB6vR6Px0N3dzevvfYaLpcLgFgsxsDAAG63mxdeeIG//uu/\nzvj9ZDJ51f5bt27dZR/Dx3X8V/PYtfFf/n9X+/gvBWMKjgMHDvDKK68wb948vvjFL/Lqq6/y4IMP\nTvggpaWl3H333Rw6dAiAyspKvvCFLwBwyy23oNfrOX369HkMX0NDQ0PjUjOm4NiwYQMnT57kvffe\nY/v27dxxxx28/PLLKdtkk3B9fX0EAgEAzp49y969e1m8eDEA9fX1vPrqqwAcO3aMSCTC9OnTL/hk\nNDQ0NDQuPhPK4xB+h507d3Lttdfidru5++67+fznPw9AT08Pd999t3x9xx13YLVaue2227j33nu5\n8847AXj44Yd59913ufHGG/niF784Shh9FLDZbJd7CBfE1Tz+q3nsoI3/cnO1j/9SoEteKqPYBNHp\ndJfMXqehoaHxUeFSzJ1a5riGhoaGxoTQBIeGhoaGxoTQBIeGhoaGxoTQBIeGhoaGxoTQBIeGhoaG\nxoTQBIeGhoaGxoTQBIeGhoaGxoTQBIeGhoaGxoTQBIeGhoaGxoTQBIeGhoaGxoTQBIeGhoaGxoTQ\nBIeGhoaGxoTQBIeGhoaGxoTQBIdGTqxpWXO5h6ChoXGFoAkOjZzYdWzX5R6ChobGFYImODQ0NDQ0\nJoTWyEkjK2ta1khNwxf0UVVaBcA919/DplWbLufQNDQ0snAp5k5NcGjkhKXZgnet93IPQ0NDYxy0\nDoAaFxWX1zXqf/FaQ0NDIxua4PgYIwRFLoLjnuvvuXQD09DQuKIxXu4BaFweXF4XW97ekvP2mk9D\nQ0NDoPk4PmYIrcLZ4aS9t51ZxbPoHe6V/wPUVdVhMVtosDZgs9gu74A1NDQmxKWYOzWN42OGzWLD\n4/dgNpnRoaN3uJfaWbXUV9fjDXixmC3YbXaa3c2a0NDQ0MiIJjg+Rri8Lh7a+RAnB0+SRFmR6NDR\n3tvOlPwp3PmJO+W2zg4na5esvVxD1dDQuILRBMfHCJvFhu9JH3aXHYAN+zfw7WXfxhvw4qh3aBFV\nGhoaOaEJjo8J6oipxtZGamfVEkvEcHe72dO1B4vZgrvbjX/Ij9lkptXXis1hA6C+ul7TPjQ0NCSa\nc/xjiN1lx2ax8Z1Xv8P+h/djd9mx2+wp29gcNlwNrssyPg0NjfNHSwDUuGjYLDb2P7z/cg9DQ0Pj\nKiQnwRGPx1m8eDH33nsvAD//+c/51Kc+hcFg4K233sr4nVAoxG233YbVaqWmpoZvfetbo7b5wQ9+\ngF6vp7+//wJO4aODMCVd7BLm6dFSmaKn6qvrL+oYNDQ0rl5yEhwvvvgiNTU16HQ6AG688UZ27tzJ\n7bffnvU7JpOJffv24fF4OHz4MPv27aOtrU1+fvLkSfbu3UtVVdUFnsJHB+GHyFbCfLKc17kIDs2n\noaGhkY1xBUd3dzctLS088sgj0m5WXV3N9ddfP+7Oi4qKAIhEIsTjccrKyuRnf//3f8/zzz9/vuO+\nqml2N2f9zOV1EYqFsn52JTDW+DU0ND76jBtV9eSTT/LCCy8wODg44Z0nEgluuukmurq6ePzxx6mp\nqQHgP//zP6msrGTRokVjft9ut8vXNpsNm8024TFcLlZvX83O+3dm/EydIyG0jJbjLRzsOUiBoYBw\nPEzFxgpMRlNKCfOW4y2jnNiXCpfXJTUTLcdDQ+PKweVy4XK5LukxxxQcu3btYubMmSxevPi8BqbX\n6/F4PASDQVasWIHL5eLWW29lw4YN7N27V26XLQJALTiuNlqOt4z5ecGzBYSfCaeYiYryimj1tQLw\n2M2Pyc9E3sXBnoPytc1iu6SZ3WrBoaGhceWQvqhubGy86MccU3AcOHCAV155hZaWFkKhEIODgzz4\n4IO8/PLLEzpIaWkpd999N4cOHaK8vByv10ttbS2gmMI+/elP8+abbzJz5szzP5MrjGgimvJ3s7sZ\nZ4cTQAoHm8OWkiPh8rqoq6qTn4vJWmgZze7my6ZxuLvdMq9Dy/HQ0Ph4k3MeR2trKxs3buRXv/qV\nfG/58uVs3LiRT3/606O27+vrw2g0YjabOXv2LCtWrGDdunXceeedKdvNmzeP3//+9yn+D7g68zhW\nb1/NPu8+AILhIKUFpQAstyxPMVvZHDZafa0k1507v2VblvG7939HIpkgnoxTYCgAYJppGgVG5fWl\n7sKXnjS4rm4doJiqPI95LupxNe1GQ+P8uOKKHIqoqp07d/LEE0/Q19fH3XffzeLFi/nNb35DT08P\njz76KL/+9a/p6emhoaGBRCJBIpHgb//2b0cJDfU+PwrMmToHs8kMKIJDvJ4zdQ6gmKci8YjcXteo\nnHu+IZ/wM2FAmTTv2HoHoWdGO8gLni3AUe+4ZJNqujlMaDsX20mvCQ4NjSuc5BXKFTy0nNDZdVk/\n+9EbP0piz35+2Eku/X+Xyr/3vbcvmUwmk4ZGQ3LdvnXy70vJun3r5OsfvfGjS3asyWKsMauvtYbG\n1c6lmDu1WlUXiUJjYdbPrBXWUe+ld997/eTr2F12zCYzzg4nFrOFeDIut7Pb7Jd0Va4+1sXwaaSb\nxdTHnYzzTI8EU2s1B3sOXvD+NTQ+TmiC4yLxPxf/z6yfZTL1ePyelPeTJHF5XUwzTaOjrwOPX/Ep\ntPpaKS0opamtCcicvHcxuNjHyWYWmywCoUDK3w6PQx53PNRCptndrAUDaIzLR93cqtWqukioHddC\nIDS7m7E5bOcmLYcNm8MmJ6N4Io672y2/13aiDec7TnqHe7GYLRQYCphdMptgOMiSyiWjtJQrJUHw\nSkFcb5vDRntve8r13vnHnXzu5c9hetZEOB7G9KwJ07Mmlm1ZNuo6qv/e9OamUe9paKST7f6YzPvm\nYpcmGgtNcFwCxAQfCAWUnhhBH4AMtRUr2NqKWipKKsjT5wFQXlSODh1T8qfQ3tuO2WQmGA7yUO1D\n2G32Ueaqj8pkNlkrtbVL1uJqcOFqcFFVWoXdZsdithAIBRiMDLJ07lLuX3g/OnSEngkReibE/of3\nZ7yOze5m1rSsoW+kj2Z3syawNXJGXWlhsu6VZndz1tJElwLNVHWRSFdVRRnzgz0HMegMxJNxNuzf\nQFNbE7fMvoX9D+9n06pNbFq1iWZ3M0/ueRL/U34szRYsZgsmo4mVC1bS1NaExWyZ8PGvRLKNcbLG\n3exuliYqX9BHU1sTHr+HipIKQDH7eQNeaRZMF8Jqn0uBQYmIS5KUviezyawIJ9V3r9l4Dd9c+s2s\n/hSNjy7p94xoxZxrpYWJ3CciJ+xyoQmOi0RTW9MoZ+/8svnML5uPN+Cl1dfKZ6/9LDaLDbPJLLe1\nWWx09ncCYG5SNAyhobx+8nWGIkM4O5w4O5w0WBuwVlgzOpXd3e4rarLK9FBc7Ak1EApIX4m7282S\nyiUsqVyC3WanYmMFj938GHDuuoms/HTnvLvbreTd/KndbigWwhf0sbtzN4FQAG/AK7f3D/vlRLFs\nyzKpwVxJv4XGxUHtpxNCw26zs+3wtoz3llh4CMa6T8Rna1rWsOvYLvxDfsLxMJZmC3Bp8rrUaILj\nIhGKhVIcvOrXLq+Lz738OXlj2V12XF4X3oAXm8Umb4BNqzZhabZgs9jY3bmbx25+jMbWRuqr63F5\nXWx5awsvrXop43FEZveVQqYV/aWi2d1MKBbC2eGkvbedTW9u4vTZ06xvXS+FwcptKzHqjdz1ibuo\nq6qjb6SP9wff56e//ym9w73Ek3G5v3A8TKGxkNdPvo434OWd0+9w9MOj1MyoSTmuFq318aXtRJvU\nQLoGuuT9vmL+CvmMCmGSCy6vC4/fw5FTR7CYLfiCPgoMBVjMlstSvUETHJOIWlVt9bVid9lTVqPq\nbYrzitnavhWL2YLL68JitshtXV4XC8oW0Oxuxj/kx+V10Tvci8vrwqAzYDaZsVlsUqvZcXQHm1Zt\nwhvwyptRHB8ufV2rsRDn7/K6ZGkVmLwxZjIXuLwuqsurae9tp66qDo/fQ74hny8u/KLU/p5e+jTu\nbjdfW/I17C47u47twmK2sPqTq9navpWhyBAAOnQY9Ubmls4lFAtRUVLB+4Pvc7DnoBQUrb5Wmdy5\ncttK9nTtmfTz1LjyEPed+N8X9NFgbcDj90h/ZLqwEPerN+Bla/tW+X6m+2TtkrVSQNgcNuXevkxd\nOjXBcZHxBryjmiLZbYqNvKmtSa6CxSRq3Xwux8PzmIfO/k4pRCxmC62+VpwdTlp9reh1elxeFx19\nHWxatQlrhXVU2OnlRO1jSDf/iP8nO+w2U1ivaI3b4GyQwjkSj7C7c7csYe/scNLR1wGcC3k+3HuY\nQCiQUr7BoDdQkl8io9oA/qrmr9hxdAflReXSrCj8WELY//bd347ZcXGyzFmaWezy4fA4pP+xa6CL\nqtIqXF4XwXBQCoi+kb6MZiuAdXXrRlkmsuU2gWKeulxoguM8yfSAqicth8chJy21GunwOOTN0Dvc\ny8oFKwmEAhj1RroGuqQZSky4wmwlfBotx1s4cPIAOnQkkgkpcOb+aC4vr35ZHl8kCY7FZOQkZJuo\nxDmox2Cz2OTD0NjaSF1VHQ3OBhqsDZM+2clESZedxtZGWo638NYHSrdKYXY6NXwKk9EEKMUaf/DG\nDzhw8gB6nZ5gOAjAYHiQRDIBKKVhxHvbDm8jnozjC/rkbzAcHAZgSv4UXvniK/z5v/05oWdCGXu6\np+PwOEaZ8s7nmkz0e5qgmTyE6RmU6/q77t/hXetl5baVWX//scxWY+U2Xe7ioprgOE/Gc2SFYiE5\naQlsFlvKzeUNeHHUO1i5bSX+IT+FxkIcHodctVqaLZhNZhqsDVKgrLpuFY56B7f+9FYO9hyktKCU\nYDiIXqenwdkgnWS5RF5NpK/GmpY1o3JThCDIJjgykUmQTDQJLxeENiNL09vsUpg4PA76Rvq4efbN\nmIwm9nTtweFxMBQZonZWLd6Al0g8Qr4hn/KicroGugBS6oyBolXk6fOYM3UO3YPdhONKvbHPXvtZ\nXnS/KAVOLueZyaR5KSZ0TXBcGOlagdpUDYpAEFF86ahNy42tjUqYvevcPTvW73K5k1A1wTGJqG+i\n3uFeAOqq6jCbzARCAbnSTqezv5Olc5fiH/LjXeuVzvJ0+6V6VfL8Xc9zx9Y7pEBR2z8BGqwNk3pu\nu47tShEcj+16jI41HaO2E+apdB+GiBxTR52ICTUXzmeCS9/e4/fg7HDKiJRWXyt1VXUYdAYc9Q7l\nuje4WLltJR6/B5PRJH0b86fNp3JqJR19HVSXV2Oz2Hjpdy9xNnqW/rP9hONhaZ6qKKng3//w75QV\nlslxNzhHa1WZfGKZxp3LtbmY5Vo0ziGc1NYKa8qiRCBeNzgb2Nq+VS5YMpGuaafsRyVA4NJViMgV\nTXBMgGyOLBFWJx5UcaNkMlW5u92p/TVcdvpG+tjavpVZxbNGOdTTJwW1eWdKwRQAGTGkzinw+D0p\nN5swS4m+IIFQQGZTB0IBGqwN465ixD5cXhfHTh9L0ai8AS/WCivODic2i01OynBuAhMT45qWNbSd\naMuokV2MB8RsMkth3N7bzi2zb+H9wfeZUTxDjtXhceDxe2hwNrCnaw+3zL6FssIyPH4PpQWlrFyw\nkqHIEB19HXQPduPscBJNRAnFQ4TiIUryS4jFY0zNm4rFbCGWiPHhyIdyAhAlY8bCG/DKf+M5StWk\nr1AzmUXUglcTNBNHrWGrk3nTFzTTvz+d/lC//FsESZSZyjj9zdMp209Ug76S0ATHBFBPgGpHVqZI\nCfV31AgHrJjAuwe7CYaDFBgK6B3u5ae//ykzimcwzzwv5fsu77kmT6Ku1WB4EJfXxUh0hO7Bbrmd\nzWIbZYYSf4t/dpedjr4OXA2uUcJNjK+zv1Nmp/qCPr6x9xvYXXaun369DGOtq6qTD9DaJWvH9WuA\nor2YTWa5Gstm/52sCU6cm92lZI476h00OBtoO9Emt7GYLSypXEKDVXk/Eo/Q0dchNcd/PvTPJJIJ\nCo2F9A73csvsW5heOJ1IPEIkHmEkOkIimSAUD7G+dT2g+ETc3W5GoiO097bT4GzAYrbI8avPw+V1\n4ah3pIxnIj6q8TSy9AlrIr6wjzvCvCmuWbYAFHWwh8vrYvnW5Sk9d8T72ZJeryaBrgmOi0S2VYXo\n0aGOgGpsbaRyaiV9I33cNf+ulAlEILQXm8Umk4uefe1ZfEEfpQWldA10Yd2s7FMcQ3A+NW2cHU5c\nDS5pnrI0WxiKDBFLxDh2+hhwLj/iO69+B4PekKJBODucTMmfkuJnEJpUKBaivbddaldjPUyZVHm1\nYM5lwlQ/jFWlVdhddqwVVva+uxe7zc761vUpE2nXQBezp8zmurLr8AV91M6qldfVWmHlZ//3Z5wI\nnpBCBcCoNxKLx0iSlI10YokYiWSCdwfeBRRTWbrvSYw/F59UOpl8VFfaBPNRwOFx8Mo7r1CxsUL+\n5u297Ww8sJHh6LB8HnMxLalr0anJNF9cyQJdExw5IlYdFrOFxtZGVsxfgc2hPPBb27dKE4P4G5Qb\nTsRxixIBwjwE5yIjHB6HYupIxOQEkm5a8Pg98hitvlYeqn2IeDJO7axa2nvbKc4r5r2B9xiMDNLs\nbsYX9DHnh3MYjgwzFBkinoxj3Wzl3YF3icYVEwsgi/y1HG9h1XWrsFls7Di6Qwo1oRkJv4DoagjK\nqtg/5Kc4v5g9XXuwWWxSA2k53sLJwZNy22Z3M7s7d2MymuTD1+xuJhgOykizbCsrtaAQPhT1tRnP\njKP+3BvwYrfZlRyZM36sm62yjAjAwpkLeaj2IV555xXqq+s5cPIAnf2dDEeH8QV9tPe2A3DD9Bso\nyS8BYCgyRHV5NR19HTK8NxgO8ulrPk1RXlGK0E+fDKQPJM0nlYsACIQCUotJ952NtYIVPje4fPk+\nV4NTXm2aHggNsK5unQx/f+zmx1JCvTNRZipL2Q/Anq49V2R+1UTRBEeOiAdRmqdUN0u6WUGdES5u\nDnXijtrp7fK6iMajeANeguGgjMwIhAI0r2xOUY8d9Q5ZME0IGLPJTF1VHR19HZiMJhZfsxhXgyvl\nOOYmZaWsPq7dZcfhceBd68XSbOHNR99UTB+719LR10E4Hk4RcAtnLuTIqSNyH4b1hpTxiQdInPOp\n4VMp12/tkrXSoSjqbwkz2XgrK/UkI3woEyE96VCY+RIkpOmwe7CbyqmVstzLSHQEa4WVO+bdwdNL\nn2bt7rXyerf6Wqmvrmfzoc3SVCUCAfToSaBEUx07fYzrp18vf89M5oeW4y0Ze6tkO0d17/r23nYC\noQAWsyUlIzn9GILxQkIvJWIivRInzmZ3M9YKq/J8/GmhBvDS714i35Avk3EzJfiqOf3N00DmMP3x\nuBKvixpNcEyQXKOAckFMAjOKZ9DT2yP7iddX10s/hjrnwe6y0+5vl9oHIFe5ak0A4Pjp47KOTTAc\nxN3txtJsGbOmjfBRODwOJTKqwSW1nfcH30/Zdppp2qhVo8vrYvuR7VKTAEVomYwmPlP5GWoranF5\nXfiCPumwVmfLjxe2KrbPFPYI2Vdw6qCFprYmQrEQ3YPdDIQGUq9Z/3H6R/qJJ+NEE1Hqt9djMprY\ncXSHzNYHpMaz7fA2ABmuW2AokCG5hcZCTEYTX7rxSzJwQiw81GHBB3sOTmgFmp49LLSYiZSvuFK4\nUrUOYQJUL4q2Hd5G5dRK7DY7TW1NPL306ZyjAtXC3hf0jbI4ZOJKvC5qNMExDpnitEWZkEx2SaHa\nZosYUmeRp08CouxIIBRga/tWHqp9CG/AS+2sWlp9rdgsNmoraqmz1LF2yVpKNpQwt3Qu7b3tcqL2\nBX2UbCihKK+Iu+bfBcDW9q3YLDaWVC6RY+0b6SOaiGJz2KRZ68PhD0kkEzJBztxkpjivmBnFM6SJ\nR/DM7c+kOAlF5FBFSQUVJRXnMuErrCycuTBFWLm73aMevHTBkX7dRQ8TX9CXImDVk+14D5vNcq4X\nCigTvpjoS/JLKCsso+dMDwWGAgBMRhPV5dUsKFvAgrIFWCusePweWQama6BLmi8CoQD11fVsO7yN\n7sFuZk+ZTddAlyxIuXDmQsqLyuU4XnS/SG1FLaUFpTlFQWW6NkKIAhl9RblE8FzKCSrd3CsWQ1eK\nySaTEBAavgjLTl/c5DJutbZt3Wy9bGVCJhNNcIxDpps6k0kg2/bp244V8irs3HKlo4raElFIakry\nS3jz0TflalP4VNLjwUsLStn9wG7g3MMhSriDUuakvrpehoGKlbN4TwgIEYqrfvjV560WhJZmpRCb\nMEepUQuwsSa0TNdRvcoW5zcRs4c6dr7s+2U8cdsTrG9dT4O1QTr0z0TOAEoujghiEOMRAlJtsmxw\nNkgtsGugS+Z8lOSXyElCHZe/pmUNznec7PPuIxgOZqxwmi4k1BOW+Cd8FWNFpY2nxV3KCVstuEWE\noFiEXeqxpNPsbpZCQq0VvP3B28ybNo+hyJCyeHPYZKBErkl4q7evpraiVv5mHwU0wZED6TZywURW\nStkeXvG+WD0L1LbTbCrxfTX35XTs66dfP+Y41ELJP+TH4/cQiUf4P//3/3A2dlZuZ1hvIN+QT1Vp\nFY56h7TPqxFqed9InzJ2hxIFpn7Q1CtpdV6MuL7qKJVspqfxrk0u3DrnVkBp0ytWhV/42Reoq6rj\neP9xHr3p0ZRILiEA1cUpXV7XKIHfYG1I0ZjSx1peVI5BZyAYDlKcVyyF2XiBAemIsOrxyPSbi/eu\n2XgNHzz1wbj7mCyE/0+dEHmpfCxjCVBhphXbuRpcNLubcXe7MZvMDEeHpTlYtDPI9Zj7vPuoragF\nGFW3bjK4HCY/TXDkgNpGPpG493TzwliCI9PqRT35ZMosLi8qlyaxh2ofYmbxTJmhrTbzpJcySEd9\nM3f0dXBfzX109ney+4HdckXv8XsIPK08WGIcIixXTWd/J96AF71ODyghqMFwkB++8cOUlbs4D/Fa\nPaGkhyOrxy0ErMvr4ht7v0FRXpEU5sKMOFbtK3G9ReVhYXsWAm44otSb+vpnvz7KFJeuAYl9rd6+\nmp3375RCUPilbpl9CzaHTZY1+dmRn/Fe4D2MeiPxZJwCQ4GS4+FvT/F9iN9NRPAIZ2wmoZpJIKRH\nUqm1lvTr4B/2p/x9Mcg0JpEQeT5hyBcyjkzXK90JbjKaqNiomFyFKbO0oJTyovKMGuRYx7O77IRi\noRTtfLKvtSY4rnAmurq1WWw53WCZyLTiHstENp75TO1oV7+fLrTMJjObVm2iwdmQsr9YIjaqqmcs\nERt106pNYKZnTQSeDmBz2Oge7M5J4GaKUlHvX621FOUVpZjC1JP5ePzi6C9kRBXAa77XZFKjyAxf\nOHPhmPsQx9nn3Sf/Tq+PBcokI8xz6gl0SeUSRqIjfG3J1+T3081TIgPfZrHJ5MV01Kae9H0AMpJM\nbepSf0+8vliTT/p9nH4dsm13KVAfU9xHfSN90h/V3tuOzaIsKlYuWJnzfldvX83urt0kk0nC8TAF\nhgKa3c0styy/ovMzckUTHBNkIje21FDSHOXqOPpc3r/Qh0n9cKQ7uWF0iKdYfTe7m6VddpppGnab\nnTUtaygtKMXhcTAcHZYCJlO0VlFekXwtenWPF0VyISvQTCtr8f7f/PxvePyWx+V1DcfDVJdXy3M2\n6o1EE1FAcYQGQgEWlC0Yc6zn85nQdgoMBRnDcNXbpddAymSWElpOpgQ09T7UwuLe/7iXoeiQvBai\nLMa/HPqXi2q2ShdmNsu5vJ+LtQLPpO3A2EJKfR9vPrRZtkGwVlizdolM1+b+2PdHZhXPkg2XROfJ\nyTrPsc7rUqAJjgmSyw8z0QfkQuPrxzKFpb+XaUUvnNour0vanhtbG2WxQrPJzM77dwLKQ1VeVI7d\nZsfSbMG7dvT+xLmX5JdIJznA0//9tHTgp7fMfGzXY9y/8H7pPxBjTx+/WsipnZVL5y6Vk2ima3D6\n7OlRJU7E/7pGHd9e9m0ZhrxywUoc9Y6sgkiYp4SmEQwHMTcpTs/lluXsvH+nDDeG0eZCZ4eTeeZ5\nOQsgtWkvE9lyCcT2IltZNhP6kxASvqWHah9ia/tWvnLzVy5qlJM4tsfvkYLzfDXysUiPehpLM08f\nXzoiyTPdPzdWPobNYsM/5MdaYaXB2qDUpJtEoSGOket5XQw0wXERyPSjXsgDMp6fZCyTmHrlJcwd\nY41bncyYrRZXpuOk102y2+z8xQ1/wZFTRwAlUkX0vhCIUu3iQTObzKOa2ajLuYtjpCdTpgcvjCV4\n3va/DUCnjX8oAAAgAElEQVTL8RYcHodMVPzea98jkUxg0BlwdjhTmlBlOl8hSEEJWxb+H0G6ZqU+\npyOnjrDjr3eM2qeaTBOf+rcURRlFJruoYSYmK/X3/EN+ud9MNn1vwKsItguYfNT35pqWNdxXc1/G\n67bj6A55T6jP6Uog01iWzl2adftMGrRY2IgoLG/ASzKZHFV09GpHExyXiGw3TS7v5xJWORbCXGaz\n2MZd0adHfWTy6YjvqDuQqesmiRVueVE5C2culL6EYDhIIBSQGoNoeaveR/p41OXcM9VmStfuRI6H\nzWLjiZYnONp3lGQySYKENMmUbChhODrMDdNvQK/TK3WmEjGK84opNZXSc6ZHjtNaYZ2UhldqTen0\n2dNSGI/XkGe8xYDIZK8qrRqV1yKujX/IP6okidC+fnDgB5Myoan9LLuO7UrJWRGs3r6a/3znP9Hr\n9MST8VFa2oUefzyTVK7nqd6X+nnpG+lLWQjZXXacHU7qq+t5cs+T1N9QL7XQeDKOQWegb6SPyqmV\nF7V/xuUQSJrguMiMd9NO9H01udpv0+3dY60sxQ2uLkWeLQRZPekHQgE5GYp6PNsOb5OJU4AMZ+wf\n6WftkrV8fe/X5eQRDAdpO9HGoZ5DtPvbc5pI1CHMDo9D5gaIa/PSqpfk5/pGPT9c8UOe+q+nGPr2\nEDaHjbc/eJvyonICoQAmo4nB8CAzimbQc6ZHNtRq723H2eFUSq5nidZablmecXxqZ7TNci4IwbDe\ncMFJYGrhUFVaJc2B4n5Q10frHe6VbUzFwkDUBxO+jvnT5uMNeC+oG6N6YdNyvCVlrM3uZgZCA1Jo\nGHQGrBXWSetkl4vpJtfzymlfDuVz8Ts+uedJec8u+ski3g28S0VJBWtuXTNmO+fzWQzqGnUpVXev\nWMERj8e5+eabqays5Fe/+hU///nPsdvtdHR0cPDgQW666aZR3wmFQtTV1REOh4lEIvzlX/4lzz33\nHABf//rX2bVrF/n5+cyfP59//dd/pbS0dNQ+PgpcyEMobPZj5TiMdYPnEpqZ7bW6FLla6KhJd6qb\nTWaOnDoiTVKi77Io/AfICc7cZCaWiKHX6RmJjqBDRzwZ5+bZN7P/xH7Kvl/G1IKpionrWWV/6fWz\n1HZnEYprMVtGJWct27KMJEn+fs/fkySJcb1RZscPRgbRc65VbCAUYP60+XQ+0Yl1sxVrhVV2ahSk\nd0Pcef/OjCXO1YJD7RNJJBM5r7bHmljE7yKS6tJX1+pkTGGGFFUJ1tWtk9uKFbUIfshVu02/R/P0\neTS1NRGOh/EFfRzuPcxPf/9Tbp1zK290v4HJaJJhyJF4hIUzF2YVGunX+EpB3PMizFzcj2qO9h1F\nr9PTf7ZfFhxVFxFVX9+mtqar0oSVk+B48cUXqamp4cwZJaP2xhtvZOfOnXzlK1/J+h2TycS+ffso\nKioiFouxdOlS2traWLp0KX/+53/O97//ffR6PU8//TTPPfccTU1Nk3NGHxHUk4C6aOJE7NDp2o6o\nCpvJDKZ+nS5wspWGUPsbFrykRCANRYYIx8NsPrSZAkOBFBR5+jyun349ZYVlvPreq6xdspbG1kae\nuO0JADYe2MjNs28eVfhQ7YBPLxApxqw+X4fHIU1M4vM7P3EnbSfbmFk8k97hXqYWTE2pU5UggShO\n6Av6mFU8iwUvLaBroEsWfLRutmI2mamvrh/VDRFSzWiZJl21cNA16kb5RLIh+kCMF/yQaUEgro34\nDdTbZ/NfiXPLRXDIfbnsUttbMX8F/9X1Xxj1RsLxMIV5hZSaSqkoqaC+up7G1kYqSiowGU0ygTWT\nkMh0jQXjjW2yndBqxD3f7G7myT1PpmjjwhQKsKBsQUrL2FPDp1KeXXEO6pDwq4lxBUd3dzctLS38\nwz/8Az/84Q8BqK6uzmnnRUVKOGYkEiEej1NWppQZvuuuu+Q2t912G7/4xS8mPPCPApMZO59tP+k+\nirHUZvW+Mpm6xkIUgGtwNmA2mfE85pGTkrvbzf4T+zl2+hgj0RGiiags79DuV5y7d33irlFFB7Oh\nFnDiGEIbE4llwkEvTDZqc04gFECHTuZtALKiLSjC7zOVn6Ekv0RqHGqBJcpTZPK3yGOPEYIN5ybs\n8a6tugik2F59DdTvqetwqffr7nanaEzqcYr9i1YAIpIuV9THGYmO4O52kyRJNBFFj56ugS7KCstk\nkcj0AAcYW0hkQt1UadmWZex/eP+oMU0W2fa1dslanB1OGS2la9Sh1+lJJpMkSfLO6Xd45/Q7cnuj\n3pjSQVPcI8JXIo6V7XhqoZT+d3qzqEvBuILjySef5IUXXmBwcHDCO08kEtx00010dXXx+OOPU1NT\nM2qbLVu28MUvfjHj9+12u3xts9mw2WwTHsOVTK6CIxfn3lifpR8nXaNQO5fV+8vVqbimZQ1vffAW\n9dvrCYaD+If8WJotLK5YTKmplI6+Dv5Hzf+Qnfe2H9lO4OkAukadErraZMZ5v5MHdz44asJdXLFY\nvq6vrk/px6Eeo8loSllhiwlw4cyFUotxdjjpP9vPHfPuSGnNCkq1XyG4FpQt4L3Ae7T3tsteJ9Wb\nqqUvxBf08fR/P8361vUU5hUyo2gG7b3t1G+vx2wyK7WOLKkh2M3u5lEJd0BK1d3030dEPglfk+jT\nrv5N081FmQIfllQuSdm/+ne1WWwyN2fzoc0Ew0H5O6abVzKNEZBZ8sdOH8NkNBEMB8nT5xFLxDDq\njKy6bpW8z4rzijOW63B5Xew4uiOl4+TU56ZSVlg2KkdIHX58sOfgqH1dKuqr61Oerfh3FfOncb2R\nz33ic+zz7iMaj5JEaer1jb3f4Fu//RYzimZwcvCkPA9nhxOzyZzxXhCohUO6j8PlcuFyuS7GKWZl\nTMGxa9cuZs6cyeLFi89rYHq9Ho/HQzAYZMWKFbhcrpTJ/x//8R/Jz8/nS1/6UsbvqwXHx5lcJ/Bs\npqbaWbXSB5DeyxpSzRZjZaBnO6bI7QBlNSgedJfXxdrdawnFQrLZVduJNuLJuFwx6Rv1JEnysPNh\n9Dp9igN/Tcsauvq7Ulapoorwym0rCcVCUtiVFpSi1+kZCA1QWlAqJ+8jp47g8rr47bu/lZoQKM7b\nD0c+lJqHWtsRD/S1U6/l4cUP4+52yyKRQEo/EXHdQJlI1CZFdT0mtYnD4/fk5BxXa0+tvlYZgJBJ\nGxTmomxO4XSTnppMuTnjmUXVwkT43Wpm1LC1fSvTTNMYDA+iQ0csGZMOfLPJzL//1b9js9hY9JNF\nDEYGCcVC9A73yhL2n6n8DDvv3ykLPwozZXr0nLjmiWSCy8VYtcKWVC4hFAvxRvcbXFNyDUBKzpO4\nvpZmi7wnz5f0RXVjY2P2jSeJMQXHgQMHeOWVV2hpaSEUCjE4OMiDDz7Iyy+/PKGDlJaWcvfdd3Po\n0CF5gg6Hg5aWFn7729+e9+CvBsZb7QtyMQfleiz1qlTYntfVrQMYFa8/nllivDGJ6BwxuZUWlHLk\n1BHWtKzhyKkjBEIBeod7MegMeANe4sm4DIGNxCPkG/IJx8OcGjnFSHREVhIFxYQhyn+0elupraiV\nq3CxmjcZTYRiISkU66rq8Aa81FfX46h3yIn2D6f+kHLeX73lq2x5ewufmPYJWn2tGHSGcw7z8CB6\nnZ58Qz5wbsWu7oboC/qY84M5zCiewVBkiL6RvlEBCOm/jcihCIaDY5onnmh5gsNfPTzqN2j1tcpa\nZOocDPF/tm5+QrCM9ZteaNVWcW3d3W6O9x9Hh06a/9zdbplIJ449GBmUAqqprSmjk1wdVJHOs689\nCyhhr3nfy8OgM3DL7FtGma0uBpme4bqqOvn81ZTXYDaZeW/gPaLxqGzvnMlPeLVWyx1TcGzYsIEN\nGzYA0NraysaNG0cJjWQys32tr68Po9GI2Wzm7Nmz7N27l3XrlMlr9+7dvPDCC7S2tmIymTJ+/6NC\nuuBIv3EmI+Mzk8lLvRK1WWxZnaHpY5voccVqU0zmI9ERTEYT5UXl2CxKBvo00zQi8YjMxhZl2Vt9\nrTxy0yP8+OCPKS8qxxf08bb/bd7ofkNW1wWldtSBkwcYCA3gC/oozivGG/DS7m/nt+/+li8v+jKO\neseoeHth9nJ3uwnFQkoXRO+55LmTgydle9uFMxdSUVLBnq49TC2YKh35wgQksNvs7Di6g7YTbXgD\nXh6tflSGs4qeC30jfVlX+S6vi+ry6jFLoR/58NzYa2cpQlRE0qWX2hefbW3fKredyCJE3DvC/yJy\nc853EXP/wvtlVJu5STG/eNd6cXldKQJiMDwoJ+BYIiaDMI73H2fzoc0MR5WWxxUbFUe62lwlOleC\nUg9tSeUSqcFNpt8wG+L6pmtB4u/bLbfj7HDSO9xLkiTTTNMoNBammKLE/+ktg3Phcvg00plQHodO\np5gXdu7cyRNPPEFfXx933303ixcv5je/+Q09PT08+uij/PrXv6anp4eGhgYSiQSJRIK//du/5c47\n7wTg7/7u74hEItJJ/pnPfIaf/OQnk3xqHw+Ek01M0pDaKGk8M8X5CI70B8ZmsVE7q1Y2MxIF/Nzd\nbmYVz6JvpI94Ms49/3EPI9ERkiSlo/pffv8vAMwsnkkgFGBxxWJ2Hd/Fjw/+WB5PdNhzd7uZZppG\nIpmQ348lYymZ0epzftH9IgWGAv773f+WxzfqjSy3LJfOe2/AS1d/l1yp6hp1mE1mivKKpHkqvTWw\nzWJjzg/nyLBjodVVTq3E4XGwu3M35UXlKb+Hf8iPf8gvcypsDlvGHIYdR3ekmOwyTVIC0TZWXdwx\nXZtMzzAH2HZ4G0vnLqXB2pDiUIdzuTkTmXzV1QjShVZ1ebUc273/371UbFQijQZCA3zu5c9JLc+g\nM0gtqnZWLdYKKzuO7pDlX9Skr9K7B7tTzvliCY5si8D0/Kg1LWvwBrzEEjFAMaflGfJkW2LxXRi7\nP8+VTM6Co66ujrq6OgBWr17N6tWrR20ze/Zsfv3rXwOwaNEi3nrrrYz7On78+PmM9aohV3PU+azw\ns2oWqptXNDtKP8ZkZECrjylIn+TU74lEwFgihqPeQVNbEx19HVJzeOqzT/Gvb/8rbz76JpZmS0rY\n6pTnpjAUGSJPn0c0EZVJevFknNpZtbw78C5wzpSUPq46Sx2lplIZbbV2ydqU/hqgmHj2P7yfNS1r\n2HVsF3qdHl/QJx3891x/j/TfCGH0tv9t/EN+EskEja2NlBaUUpxXLCNsbI7MyZY2h01qG+m/46Kf\nLOJo31E5kRrXG0kkE/zy6C85/NXDGX/3dEGS7d4Qv4O4f9zdbrnSVTvhM30vG2Pd4+I4yy3LqbOc\nC+Ueigwxf9p86qvr2bB/A5HvRGhwNqQ4/dWRccPRYbkwUkeltfe2s3LbSt58/00MOgNdA10ynyK9\nrM1kku35UyPuI3GfijGN1bL5akTLHL8I5GqOuhDBoa6lNJZtXX2cbCU7JmOFln7Oze5mHB6HXGUD\nNDgbGAgNMM88T7bJVV+be66/J2U8xXnFDEWGKMorIhgOYq2wcuTUEU6fPZ2ygm5sbaSxtRGjzsje\nB/disyiJVe5ut0zsM+gM/PCNH0ozijiO6J+xadUmFpQtwOFxyNwNi9nCkVNHZHl1l9dFbUUtO+/f\nKVfOvcO9OO93cufLd6aEBosKw+oS3aI0SHpoLcDhrx5m9fbVvO1/m+7BbmLfHV2yXkzK6bWmRJa7\n8DeNNfG3HG+hs79TTua+oA+L2TLhKrXpznG7zc6yLctS6rKJxEhvwJuSKOoNeGXbYo/fI82aIoBD\ndFOcPWW21CxFHS510ufa3WvxPOahZEOJFIiiakH6GC8W6Qu0Tas28eODP5a+tsHwIP6n/Nl3cJWi\nCY6rFI/fk6LCixXZ2t1rZWl0GL8WUi6CY6xtMmk24oFVFyM0GU3sfmA35iYzDVal1erCmQtZ07JG\nli/ftGqTfOhdXqUOU+9wL+VF5YRiIblStzRbCMVCBEIBwvEwK+avYOWClaPO01phlZN1obGQjtMd\nSumQP7V6BWS0lxiz5zHPqKgpQcvxFo5+eBSHx0HvcC86FNPtjqM7MOjOlRFRF18U5yLqd4lS7mrE\nKvVE8ITMLRHaznhahHhP5K2kIxYT7m43rb5Wqb0d7j0sk8/ENVJnoGcim8YqrqU6NFZEj1nM56oj\niwCE+up6fvnHX8qsdxGeK65Re287VaVV5Onz5P7UIbjegFdqrjaHTWomHr+HFfNXTGqlWLVP7nys\nB1+6MXPE6NWOJjguMhe64slmEhAPq0CEewq1Xkxize5mKUTEShjGFyjpYxhPcIz1ntCM7C67LHQI\nSuvbBmdDioMwPVchFAvRNdDF7JLZODwOHtv1GCcHT6aEYe59dy9/6P0D1gorTW1NHDh5QPYNB+g5\n08Md8+6g43SH3P/W9q3UV9dTV1U3yodxz/X3yDBeQK7wxcQ4HFS6BIpJXvhjdI06CgwFFBgKUq6D\nMF/VVdWx+4HdLPrJopTjbVq1iftq7pN5IFPzp3L4q4dzTsQbazsxBrPJjMfvkUEHImx5mmkat865\nVfrEsu3fZrFJjTXT/SCc3A3OBhndBorz1xvwSu1rODrM1vatsqeLyD1xeBx09ndKv4Av6KN2Vi2h\nWGiUJru7czcmo4lwPCxza0QNrvSclckgU8XoTDS7m9n0pmKOavW1MrtkNkdOHZk0E/GVhCY4LjIX\nKjjSTQKC9t52ueIcK6QvvQz5icAJttRvkbbjS9GBLZFI8G7gXTYf2gycy3AWIZregDeln7d4LxAK\n4Av6qKuqw2Q0cfTDo5waPjUqdr84r1j2Rt/9wO5RWbbRRJQ9XXsApBnkhuk3yM/ViXNmk5kFZQtY\nULaAprYmdj+w+1wEjLOB3Z27qS6vTvEJPFT7EP/xh/8g8p0IgMyKFyzbsgyD3iD7lB/tOyonX1HD\nKhQLEY6HGQwPMrVgKqu3r6bUVDrm76EWsOrt0oMhvrH3G5wIniAQCqQ4o6PxKMuqlo2anNWIhUv6\n/j1+D4FQgC1vb+Hk4EnufPlOEskE/3b43zDoDLw38B7L5y3H4XHwH3/4D9kkS3AmfIau/i7qquow\nm8yU5JfwwKIHpNYhii1++Rdfxu6yy+s9q3gWQ5Ehaf48ETxBPBHH4XHQ3tsuEz8nq3hirsI7/Z4D\n6BnqoWeo54ILWl6JaILjI4K6adCK+SuyNuTxBr2j3h/LxCH2KTgf4XLT7Jv4Qs0XAGhqa2JxxWLe\n6H6DbYe3EQwH5Qr0+unXM3/afCxmC/4hP+297TLEVKyIF7y0gLmlc/lj3x+JJ+IY9UapXXz+3z8P\nSWSZ9NKCUunjEATDQfme2WROae60pmWNPPf66no6+jpYvX01X1vyNelE7h3ulaYVkW0uImjSy2mI\nvI/XT75OkiSBUIAGZwPJZFJOxiIYwO5S6og9sOiBc5FNjnPXOZd8IHUCqPr3WnXdKnme9++4n6HI\nEEa9kWA4SEVJhZKQmEWrFCYhkTsifDczi2ey6rpVPLz4YRpbG5lnnsd7gfeYZ57HUGSIE4MnpOYx\nt3QuZYVlFOUVyRa9c6bO4c5P3Mm2w9toamuSZkdQFgLbDm/DZrFx1/y7ZMUBoS2K0N1DPYdYULaA\nL934JZkcOlmTtLi+6Tkz6uukfi1CZG2Oc+1+P4oCQ6AJjqsIdcnrbDZtyK5O11fX85rvtZyOlauD\nfyzUD5/QHKLxKKWmUp5e+jRrl6ylZEMJsURMTrYi76Kjr4OHah9KaYm6cttKTgRPKHWQdHoqiiu4\na/5duLwueod7uWH6DbI/RXFesRQQojmUrlHHurp1yoR+vzNlonV5Xew6tgubRSk8J6qa+of81FbU\nyuvh8Dhw1DswPWti0axFSvx+g4tlW5aNKqWh7k9i0Bk4cuoIiWSCJEkaWxtZ37qe8qJyvnrLV2ls\nbaTAUMBQZIgGZ4PMcxGv1Sa88X4fIeSEJpfei2Pp3KV09nemJE6qNQv1NQnFQjILX/gs0ptlbdi/\nAf+QHx06uga6mD9tPtXl1VjMFhqsDbi8rlFjODl4Ur5XVVpFkiS1s5QEz2A4SEl+ScaoMbvNLs8t\n35DPwZ6DimB02XOqw5YL6eHmarNbuvC2WWyjulKCEiL8UTRRCTTBcRWhvmkn0pv7E82fwBv0Aopd\nXt+oV/ZRamFL/ZZxv5+tNelYiISvipIKmTmr7pWxu3M31gorsUSMcDxM30gfm1ZtouV4C6FYCKPe\nKIWGKFR4qOcQ0USU4rxiQrEQUwqm0HaiDf+Qn3A8TGd/p4xgyjfkMxwdRoeOxtZGfvr7nwJKZFl6\nRrLZZOaRVx6he7B7VA2raCKKy+vil0d/ycM3PaxMFg4b4XhYOpeb3c3sf3h/SjmQZVuW4X7fLW32\ngFIET1VYcWbxTB67+TH5nXZ/u5xo1dc+l986W85GIBSgdlatvI5nImd49b1Xpemoqa0Jk8Ekr7/Y\nV/rEaTFb5LUVeSsWswW7y06SJMPRYa6dei2fmPYJAqEAf/zwjwDSzGQymvibT/0Nr3lf473gewBM\nyZ/CSHREOs/FuI16I96AV0ZJiUgrIZiFQBLakBCak5WFnUmrHq+HjRAQ1s1W2nvb6XyiM+v2HwU0\nwXGVki3jNJO54d2178rX+ka90vtiXfwijewcHX0d0lkpekGcHDwpczganA3SPPG/f/+/+cnBn6DX\n6TEZTQxHh3F2ONl2eBvXlFyDu9tNLKlMwsNRxTn9zul30Ov0JJIJpuRPoaywDEe9g/rt9SwoW8DB\nnoN8t+67ADz/+vOA0rY1noyzcttKwvEwz7U9R6GxkFgiJu3/ahLJBL/r/h2heIiX3C9h1Bt58/03\nAaSd/et7v85T//UU15VdJyeY/Q/vx+awSZt90+eUtgEev4dth7dRkl9CdXl1ygpfmK3E3xsPbJTj\nGKtzo/g7PdoKUtsWqyPdjOuNlBeV43/KP6ouVbpfTTiHReSexWxhJDrCL/74C8KxsBSOPWd66DnT\nI8vKmIwm/EN+KWhefe/VlGS9kegIAAWGAhLJBM/f9bwiyP5USLKipEKGCQMy81+0/RUCx9Ks9GA5\nnyzsXBjPRKi+ZoFQQGooH2U0wXGVMpEop3TUJcRzYSLaTfr37DYlARCUMhOJZII8fR7D0WEZnQTI\nFXA8GWc4OowePWaTme7BboYiQ1JogGIDjyai1M6qpWZGDd6AVyYV3rH1DpIkZQSU+uEWbWUDoQCB\npxV/g4gUamxtlBVdkyQx6AyUF5XTO9zL2WfOYlxv5HbL7TjqHRjXpz42f3btn7H/xH46Tndgd9lp\nOd5CJB7hD6f+QCKZQIcOZ4eTt/1vMxIdIZ6MEwwHafW18rmXP8es4lms/qSSUCv6dafkhFhs0oks\nypqoGe83z/a5SJYb6/uiKyCcK+onthfX3t3tJhwPs3TuUg71HGI4OkwymaSzv5OugS7Z+VEIFCGg\n48m4Eh4cj5IgIXOTwvEwVaVVKUUeRXVdm8XG83c9LyPVAqHABRcJHItsQlogBG6zuxm7yy7rmJ1P\n9OLVhCY4PmZYSi3SVDAWY62svvPqd7IWkxP23o6+DnqHe7E5bJwaPgVAobGQAQaYPWU2gVCAs7Gz\nsux0OgkSdA92E4qFqCipkKt7UMqeH/3wKMdOH+Ngz0GqSqsoyS+hd7iXQmMhI7GRlH0JZ7kv6EOH\n0jdBJJuJctZ6nT4l8ieejMtjFjxbQDwZZ9vhbTg7nKM0E+HwnTt1LmaTmY6+DkaiIzL6K0mSw72H\nCUVDfPnGL0vfgWhupJ6MXm5/mSOnjuANeBmODsuEQZPRJCOGJioo1OazP5z6gzy/D4Y+wNxkZrll\nedZ9yu6AquALUATK0Q+Psuq6VXj8Horyiogn4rKbY5IkJ4InlAq5KnNdPBlnVvEseW1Fkp9Rb6TB\n2iATJtX+FrX2o/a/iPyTbIEgk0Eu+3R5XbLUvxhT+tg/auiS2aoUXmZ0Ol3WAooaE8fQaMioaejR\nj2u2SjdlmJ41EXome+VSOBcRJBLhxArR3e3mkZseYdvhbbKp0/liMpj45tJvygfUutnK4d7DowTR\nDdNv4NjpYynv63V6CgwFnI2dldFRQsCIia/QWMjZ2NmUyrljUWAokBPlnKlz8AV9GHQGWfRRmFbE\n/ovzihmODktzTLu/nd90/oanlz5NY2sjVaVVKTkEQpCLCCkYPbFds/EaPnjqg4zjU/+O5iZz1i6E\n6YsGUVm53d9OqUnRHtJ9KTp0ssOiEAw/WvEj6bDeeGAjT332KWXfDS5ZUuXm2TdTlFdEq6+Vh2of\nwuP3pFQbEFF1oPhAxFhEoII6m3wsLkYNq0y/Qf738vnstZ+9rBFVl2Lu1DSOjwlq4ZDeCGaySHeq\nqleDgiRJdh3bRSKZIJqIyklU/b/ZZGYgNMBDtQ+x/ch2ivKKUvpliIkpFA/x/bbv4+xw0mBtkCVJ\n4sl4Soe/UCyEXqdnVvEs4sm4UnjvTwLtzq13UlZYxkBoQK6MDXoDiWSCW+fcKosXqtuvpiP8LNXl\n1dI8JbYXpreCsJIUqA4Rjifj3DL7FkaiI7I8SzgexuFxUFpQSjQeVcrDp5UYESvcTGYr/3BqeYtM\nHfLGI5s5xtKsOKXVf9ssNhn11WBt4As/+wLb79vO8q3LsVZYpcN+ODos/RL538uX0WUHew5SYCig\n2FgsKyGI+0jd693ldfHIK49IJ7/QPNSFA8cifZLPRZCs3r46az948X27y57SgEpk5YvM/49SfSo1\nmuDQGBebxcayLcuk3yAcD2N6VrGPj9UDobq8Gne3WxY5FJNpKBZiSv4UzkbPSt+FmOTVTZX2vrsX\no94o+2IYdMqE7n/Kj75Rz9SCqVjMFhKJBE/ueTLl2Grt4kTwBHqdnikFU1Iq6a7691UkSMhiiQIh\nQEaiI0qEz1ovxvXGrFqHMEm197bLMiTpiHMaDJ/rpBmKhTg1fEpmSQv7voisWjhzYcpqWpiLhMNb\nOF5WomwAACAASURBVN6zsXr7al4/+brcdk/XHrwBL/4hPzfOvHHM72ajb6SP8ufLWThzIb6gj22H\nt2HQG4jEI7zc/jJJkizfuhxA/g+KsPcFfUwvnE7l1Er6z/bLkvbV5dUpCa1q05RaiHQNdMmKxtYK\nK30jfRNqO5stByMb+7z7ctrXfTX3yYizxtZGnrjtiZTz+CiiCY6PIXr0E9reZrGlCIdspir1StXZ\n4WT3A7tTHjDhzGxe2aysVJstVJdXE4qFOH76OKeGTzG9aDrWCisev4dCYyE90R7pLC8wFMgoLIEo\nryIcq+tb16doMABzS+dy7dRrORM5QyKZkIlsZ2Nn0aHj9qrbZbOr77V+D51Op0RPPfomU5+bypqW\nNeTp84jHUwWHMG2pSTeTzSudR89QD5VTK+ka6Br1+YngCfIN+ZhNZoLhIAadAWeHk6Vzl46aEEXl\nWHVSms1hk1FfgvQsZlExWZh5REXgXGj3t8tufL6gj13HdtF/tp+FMxdy4OQB8g35LChbQHtvOzOL\nZ3Jq+BTfrfuuNLU1WJUABFESxNnh5MORD+kdOuez8ga8o0KCQcl/EfeOw+NgVvEsbBabzJk4cupI\nxjE3u5vp7O/kvpr7UkxudVV1ozTg9O/lavZSV14QeSXnG0RyNaL5ODQmTC4+Dutmq+x5Icuu/ynr\nWJg71rSsYSgyhMfvkRVpSwtKsVZYGYmO0NnfydzSubT3tvNQ7UP8/OjPZQinGtGtr66qTvbeKC0o\nZTA8SEl+CQ/WPpiyIgRSBIsg3ZdRVVolHeowWigsvXYpbSfbKDAUEIlHUiKGBKKooBB6t8y+hfcH\n3+f02dPYLDYZriwm95HoiMz0htTJTJT6ANjduZuOvg68a72jEtDEuYvxivFVTq3k4cUPs+3wtvPK\nMxC/e8mGEiqnVtLZ30k8GcdkMKUILnHdphdOp+8bfbLMvNAcfnP8N5waOZVy3RPJhMwrUk/2K+av\noO1EG6FYKOXaqn8TES0nTEPiPhMNpNTO9LqqOgKhwCh/ifDBTTNNk5pGMByUEWHLLctlBQExtnV1\n62QirhBGBc8WEH4mdXFzqbkUc6cmODQmTC528wZngywVIUwv249s553T70hnsCi8B8gOdvmGfFlR\n9JlXn2FB2QLeHXgXvU4vfQPCTzC7ZDYzimdIP4Av6JMTNSiaVZ4hT5nk/jRRTmuaRiCc2SksfBUT\nRTjAh749JBPABGIShXPCSvTktpgtPLDoAWwWGzuO7lC6FqpKoACjymis3LZSdjQMx8MyZ0CEfeoa\ndRk1IfXxxbhE7kMuq+zV21fjfMdJviGfSDyS87VRC5FIPJJSfFJcb2G+qiipSCm2uOClBTyw6AE5\n6c8qnkVJfgkngidkNFaBoYAllUtSwl7VgkOgNvOJbPx0DUcsdgRjBRCIBZEoKqoW3Om/yaVGc45r\nXJFkExpqe/TW9q0y8kZ0vxPJX8LMMhIdIRKPyIxfb8BLeVE5m97cRCwRYzg6jH/IT1FeEU8vfVo+\npMLnEAgFZEltX9Ano5BEH/KFMxey7fA2KqdWyjF+bcnXePa1Z1k6d2lKoUJglNDINgGryTfkE41H\nGY4OY3MozYoE6+rW0e5vl5nzwXCQqtIqhiPDlOSXyJapoJhiambUyJwFYfpIL6OxpHKJLCkPStTa\n6u2rlYnsT5NjMplM0Z7qquroHuym84lOWRNLPamOhVqbAfhM5WdoO9FGPBmnqrQK/5CfJZVLpKYT\nioUIPRNS6n1tr5d9xeGcb+b7bd8nSVL6Ns5Ezsh7QPQUET4NdaJf/9l+hiJDRBPn+njHEjHsNjse\nv4c5P5jD6bOnpTkz/3v5lOSXyIWI0DrFtRUObWHySm9HkAvpiZXWzdaPdI0qgSY4NCYNtY9DXW1W\ntGAVD6WoOaRGROhYK6xsenOTnBhFeOfuzt3EEjFCMWVSqimv4fjAcew2uyxwJ/Z7609vpf9sPwBD\nkSFlUmsqTXFMq4WGCLtNrktielYp163X6XMSHGL1nW/IZ5ppGsf7z3W33LB/g2y3LMplmIwmOvo6\nZM0nkaeRb8jnzUffxOVVemeIPJP23nasm62cCJ4gFAvJKsCgrNgXvLSANbeukdE/hc8WcvYZZRtd\no072CVnw0gK5Eg+GgzknqKnLaZQ/X47H75ECSfQPef3k64CS0R9LxGSf8GA4iHWzlUAokNKHxKA3\nYNQbqa+up723nWg8Ks1HoryJ6NkizJ2i8KQonb+kcokMIEifvNNNoqBoL0vnLpUlUIRQHooMyYne\n0pzag6XVm7qwUCP9dmn+kskqe3Klo5mqNCYd4TgEpIkBzgkBk9EkJyShpWw+tJne4d6UQnegTOp6\nnV4m+OUb8ik0FvLAogfYdnibNCVYN1uxVlhlwx2jzohBbyAcD2PQGSjJL6GipIKONR2YnjVRXV7N\nHz/8I9HEuQTEuqo6mcw3FiLvYyIIM9G80nk8aH0Qu83OnB/MoT/UT2lBqexFHggFuH769dLPkW6q\nEpMoZJ70C54t4JqSawDFmV1oLGRm8UyunXot+x/eL8Naz8fHUfpcKeF4mEg8khKAIHwUSZKUF5bz\nv279X9gsNr78yy9zXdl1spBiviGfZDIpTYlqrag4r1jWLROmHhGwIBIgxQJE3U4gk9ATwkotOCo2\nVvDYzY/RcryFgz0HpQZoMpow6AzUzKjhYM/BUX6PiXIlFDbUTFUaVx1CEIjJYnrhdFkpNb1VLJx7\nQEUdJGHOAHjpdy8xt3QunseUSKIGZwMmowlvwCtLspublBVemalMxvw3u5sJPB2Q5hJrhWI+EPuN\nJ+MMRYYozCskGj6XgNjqa6UoryijAx7O+VYmKjRAERw6dJwOnZbmqVMjp6gqreKBRQ8orXbrHdhd\n9pRuc+mlvNVkmqBWLVgltQ9do46Rf0g9F5vFlmK6mwiNyxvlMXWNOlkHTODyujjce1j+pteVXafk\nqNQrOSqhWAhXg4uifywiEo+kmAaHo8PodXqmF06XgtK62ZqS7Kg+B/X/4tg2i1KpdigyxMzimTS2\nNsreL73Dvbi8Loryipg9ZTaP3vQoze5mguFgSj2sC+0eeLmFxqVCExwak0q6uarZ3TyqZAQwKizS\nWmGVdnoROllWWMaHIx/i8rpYu3utTJADpF3deb8Tm8XGrT+9lbLvl5FvyCcYDmJptkgfiMloSikv\nHkvE8A54SZAYpV2EopmjxYw6I877nXx+2+cx6A0U5xWnRAZlitJSk0ApxBiOhWWZ9++99j0qp1bi\nDXgZiY6c63nxp9pUQqCmX9+xzCFd/V0ptbTE65ryGg5/9TCQu/0+HfWkqEMnExFF35T23nYMOgNr\nd6/FWmHFZDTR6muVeSTCfDm1YCorF6yk7UQbXQNd565RMsHps6dpcCrJnOI8c6nNJQRHIBTggUUP\nSCc4IAWx+J7JaJJta8U5TWa72Y8DBrvdbr/cg8hEY2MjV+jQNHLA5VVyK3rO9PCq91Vafa0Y9Ube\n+uAtxcGd1l/CYrZI27Y34MVsMjO3dC5tJ9qwmC0c6jnEn137Z3QNdLFs7jJcDYoTWWgZ7595n7Ox\ns9xXcx8dfR10/303Lq+LT834FLu+tEsey2wy819d/8WyqmV4g1506CjOK+a2ObdxYvAEt1fdLqOz\n1MJARGcdOXWEqQVTef//eR+Hx0EwHKTAUMCXF32Zzv7OMUuoROIR4sk4rb5Wnn/9eenkbe9tJ5FU\nanOVFZbJBLIGawNLKpdgd9nlZG8xW6SDXHQUVOcPfHLGJ0kkE/zFDX9Bq6+V79Z9l9urbufxWx6X\n211Ie1VRvr1roIvDpw6TSCYoMBRQalJ8SJFERBZPvOf6e2j1tRJLxOge7Obk4EkZyHDrnFtlaRGD\nziBLvAifksvrktFpQuMYK09CZKi7vC5ZGt/j92A2mWle2YzD45CLDCGsD5w8wPtn3k/pOmk2mTMe\np9ndfFHa0l4MLsXcqWkcGpPOmpY17Dq2i1AsRDQRlbHwN868UU52mUpSpwsS0QfE5XXRNdDF0rlL\niSaiHDh5QCallWwokd3sqkqrcHY4GYoMYXfZ2dO1R1ndu+z0jfSx4+gOqsurpTAQpqfh6DBtJ9uA\nc07zaCJKXVUdHr+H+up6nB1OadoAJVRTRFDFErFRfTzGQ+3k1qGjyFhE5xOdWDdb2XF0B/fV3Cc/\nz9QPRd2EKr16qwjnXd+6HrtNiaLacXTHpGQyi9DhE0+ewO6y09TWJE1fos+JL+ije7Cb10++TiwR\nk9Vz8w35BEIBovEomw9tTqmSq74mH5z5gE/O+GTGIIr0a6CO4hPl20V9q/rqehlkoXbOi+sgghHU\n4c/ZEFWUNRQ0waEx6WxatUnGx4v6RjDaHJDtYbVWWGXfB0BmGYvy7NFElBPBE4ASWRSOhWUPB1/Q\nx5T8/7+9cw9vssrz+DdpQwMtNFwbsaVBsJSWSiognRFtUBQGcLaOrOIMQgb10dlh3ToOipfZBtFS\nhwqtyzzrzuyMAdlZ1su0i7XAojQVxYoDpo6XihQSWkoLhSYltEmb5OwfmXN4kya9pjf4fZ6nj22S\n8+a8L77n957f5fsbjYraCsRFx4nvWFe6TrR5BXy9Qi63X8Y45Tj88/x/BuAr6hoZORIe5sEri14B\nABw5cwT1jnqRFcUNi7QlbWDRX1duqzWz10Cv1ePOnXciYUyCyLDistzvfPMOAIgMK2knwEAD21mT\nLZ5+WlxVDIvNEjbdpJLjJaIq2+11o7rJ1zs8Uh4pDC7fPa58ZyXqf10vAtZatVZkT5Vby0Vc43zL\neXH8C60XUG4txxdnvwh6zpxA5Vwex+LXxOa0iQcUdYw6aKyIG9ru0B9CicMVMhxEv9PTm03aszta\nEY2oSJ9AYKQ8Eh6PB9PGTsOqm1Yh91AufvWDX4lxJosJ31/8HjeOuxFOt1PIugO+p8/rYq6DJdsC\n5UtKPD73cfzh6B9w5qkzYvzG8o1odbcibWIa8j7OQ5unDa3uVpjrzVBGKrE8abmvoZTXLVJ1g+lX\ncaMxTjkOo6NGi7RVnn3Ea1wYY1ietBzTx033Sz+OjYrFO9+8A61aK7oM8sXQaDbCaDaKdq+BtSjS\nwPrSG5eKOEQ4Oes4iyW7lsDDPOL8+TzkMjnM9Wa0tLdAp9EJt5W02JAbA4PJgEOnD3UoKOQZcK8v\nf71XC3WgYm5Xu4nOJEikhX3ceF+tPTZ6AhkOol9ZnrQ85NNiV2RnZMNcb4Yxywjt675gaUVtBRZM\nWQDAZ0h4IJ3n5SvkCiEBMf216SIrx+a0wWq3ilanFpsFighF0O+1uXxPxk6309ea9tf14mnzj1/8\nEXHRcXh0zqNigY6KiIJOo8OS6Uvw5P4nMTV2KpImJEEdo8aOyh2IkEVALpNj4qiJuOX6WzBbPdtX\ngzBxFrYv3Q6TxYTzl8+joKIAgG8342hz4OPTHyM2KlYEyqULIJdb1xl1Is4jRaVU4UD1Aez+ardI\nKFC+pMQoxSj8NO2nPd59cPcj4IvV8NRhR5sDL9z+Agw6g9Ae27Bgg3iKX560HEDH1GEem/F4Pdh0\nxyY88PYDsLvscHlcIk27KwLVmPluo95Rj+yMbKwrXSd2XZ0R6v9FaQ2LtE6IoDoOYoizrnQdvjr3\nFcz1ZthddoyMHImk8UlYMGUB9p3Yh9b2VtQ56kR1NHeb6DQ6vP7X17FhwQa/p8bAGgxpU6FQjFWO\nhd1lR7QiWkhmJMYmwtHmwIXWC5gcMxl1jjrkZOZgY/lGlK0pg06jg75Yj5LjJWhqbRK9UHi18+X2\ny341A+988w7ONJ9BmaUMl9suw83cQj6Fy3wsnrYYGfEZQuJCp9H59cvgrpvA/ik8HmTMMobF1aIp\n0IiK8JcPvYwDDx2ATqODOl+NJdOXiF1PqJoIvuAHvq7KUyE2KhbWJ62duoWCvReoiWbSm6Ap0PT5\nnPlceQV7X+s8BgKq4yCuefiTMXcbSAXluG7TmtlroFFpUPp9qV9655TYKX6Fhot2LsJfHvgLFu1c\nJCRHkickw+604/ox1+O0/XSHrCgZZJgSOwXqGDX2rdqHyBcj8cLtL/iydf7uSjvdfFoUramj1X7p\nyDEjfMH7xpZGeJgH2RnZIpgrhcteZMRnoKK2AiqlCla7FXHRcZgSOwUpE1P8dha8F0eo2hjp0ziX\n5+CLq7ne3GtXC9+98fTmEREjxPesSFnht5MJnFdggyj+Gl+EF2oWYrZ6tjiHUBjNHY2BxWYRab98\n9xEOF12wOiOim4bD4/Fg7ty5iI+Px3vvvYe3334bBoMBVVVV+Pzzz3HzzTd3GON0OpGZmQmXy4W2\ntjb8wz/8AzZv3gwAuHjxIh544AFYrVZoNBq89dZbUKmujVJ9ovMGOaHIzsgWxVy8kG9/9X5EK6JF\nf4bjF44jLjoO5dZyUV+wZNcSsRgzMJgsJuGXj42KRUt7CxgYzl46GzSVloFBHaPG/ur9MJgM8DKv\n0DrSaXxuIqkbQ6fRQZWnEn3MY0bE4HLbld7qfOGcNnZa0EVow4INQvBQr9XDaDZilGJUhxTRQGmX\nYO/x/1bUVmD6uOmYMGqC2JH01nDoNDpRJ8HPRzzpd/EEHviUHnj+RSuLQrqopLuMYAkBXIr/cM1h\nABAZcCvfWSniU31JDuD91gkf3TIchYWFSElJwaVLvm16WloaioqK8Nhjj4Uco1QqUVZWhlGjRsHt\ndmPBggX45JNPcOuttyIvLw933XUXnn76abzyyivIy8tDXl5eeM6IGPJ01iCnM+ou1Qm3TsnxEkTI\nInC5/TKqGqvwYvmL8DKviE3oNDpkJWeh9PtSGM1GnLl0Bl7mxYvlLwK4EszlzakUcgUUcgXcXjfk\nMjnix8RDo9Jg1qRZONN8BlERUSioKAADw5//9md4mRce5hFpw9JiNbvLjj8c/QMA+IkeAsCntZ9i\n6fSlojhy+5HtaGxpFFlaWbuzMGHUBNFoifcrCbUod7VYmyy+/uYTRk3oVjV6d5g+7oruldVuDap7\n1dtjhxrHxQ+BK/EM/nmuU2UwGXDzdb6HWO6q6q6YY1/mdi0i7+oDtbW1KC0txSOPPCL8ZsnJyUhK\nSury4KNGjQIAtLW1wePxYOzYsQCAPXv2YM2aNQCANWvWoLi4OOQxCILT7GrGutJ1+Pj0x2j3tMPD\nPJBDDpvThkttl3Cp7RI+Pv0xPj79MUq/L4XFZsGZ5jPQa/Vwe92IjYrFv2b+K+QyORLGJKBsTRly\nMnOQGJuIlIkp+GHCD7F18VbIZL72tVnJWdi+dDuKVhb5YiUrixEVEYXnbntOqMPygratn27105H6\nr/v+S2h0SWnztKH4u2LsqNyB7IxsnHjiBGwbbMhMzPT1ithgw3/++D/R6m7F6399Hfur94uMHu3r\nWhFA53TnKd+g86kKS+drMPkC64HH6w7ZGdkw6U0wZhmRmZjpa8OrN/ntYrozr+7Ar6e0c2MgBRUF\nviQBsxHl1nJUNVZBZ9T5iVqGAzIcV+hyx/Hkk09iy5YtaG7u+T+C1+vFzTffjOrqavziF79ASkoK\nAKChoQFxcb6bKi4uDg0NwYOT0upHnU4HnU7X4zkQQ4N7d9/r1yBHlaeC2+vGXTfc1anbSuoXd3lc\n+OrcV1ApVah31AtBRIPOgC2Ht8DutPtcQ+2X8Xnd5zhtP42Gyw2w2CzITMxEVWMVDDoDcg/lou5S\nnTg+z7haPG0xbE4bxirHduhTwZ9sI+WRMOh8vnR1jBrGLCMmbZmEmuYavFj+otjRSNumRkVEIVIe\nicvtl7vV613q0uHGqC9S3fx43EcfKJzYl+P251gulqnT6LC/er+o3I5WRItz4v9+PFEgMTbRl5ig\n1Yt6mFDHvloMgclkgslkGtDv7NRwlJSUYNKkSUhPT+/VxORyOcxmM+x2OxYvXgyTydRh8ZfJZEJ6\nOhCSHLl6kBoH3iBn5vaZHYxGoLpoYUWhn2uL99nmdREe5sELB18QInmX2y8jQhaBBVMWICs5C+sP\nrIcxy4glu5ZApVRh+mvTRSwje182HG0OqGPUsLvsyIjPgMVmwdIbl3aIAfBFZvzI8QAg6isAXx9u\n/p3AFTcYLwRkYF1mufRWPyoUoYrdwk1/zttcbxYGghd/Aj5BRIPJAGWkUjR+4mNMFpNIIujsfK8m\nwxH4UL1x48bQHw4TnRqOw4cPY8+ePSgtLYXT6URzczNWr16NnTt39uhLYmNjsWzZMhw9ehQ6nQ5x\ncXGor6+HWq3G2bNnMWnSpD6dBDG84HpE0t4VnEBpBx4w5UVv0p7jdQ7friFmRAwiZBGQyWRQyBVw\neVwot5YLiXRNgQYNjgaolCqhjeRhHnxz/ht4mRfXj75efDffffDYBa9654uMIkIBk8WEdm87Pq39\nFMqXlGBgQoNKCi8EfHDWg9Br9X67kECyM7KDuo14Rk9nXJd/Hc7++iyA0FIk0nMI52If7kI4vruU\nSu3zVNjiqmK/a8GNt/Rcw23IiOB0GuPIzc1FTU0NTp06hd27d+OOO+7oYDRCPUk1NjbCZvOlw7W2\ntuLAgQPQan2yDT/+8Y+xY4cvHXHHjh3IyqJ/7GsF3lvBYDLAwzxCEbezgi+evRStiEZOZg5yMnMw\nImIElBFKREVEwdHmgId5wBgT/TeiFdFYPXs1ZJDBkm3B/Pj5OPvrs7g98Xa8cPsLiJBFoO03bXD/\nqxvWJ62YHTcbWclZSJ6QjDWz1/gK2v6eOfX0gadhMBmQvD0Z1U3VuOvNu1B3qQ7NrmY/QxaKHZU7\noNPounRTSTvtSc+9q8W5/vIV/393CueGatUznzuPy6hj1MjJzIFeq4dCrkDBEp/ScmJsojDSXKuK\nqzB3dm7c9cVTibvz/x4RnB7VcXCXUlFREZ544gk0NjZi2bJlSE9Px969e1FXV4dHH30U77//Purq\n6qDX6+H1euH1evHQQw/hzjvvBABs2LAB999/P/74xz+KdFzi2qCwohCONod4ui6oKIDL7cL2I9sx\na9Isn6x4kCwdAEhXp4tdAZepiIuOg6PNAS/zotXdKrr5XW6/LGol5BvlYGAY9fIotLpbUdtcCy/z\nijlwHSou4ZEYmygkuW/70204bT8t+nwDvgZETObbZeRk5uDF8heFTlNURJQwJrFRsfAwD9bMXtO/\nF/XvCME+iYQ8MLSL1YAr8+bX32g2otnVjCZnE6oaq8T1NJgMaGlvQbOrGXHRcbDarWI3sue7PUJV\nuLMMtM7SgYnuQ5XjxIDAF4aK2grsr96PxNhEWO2+J32ewcRbwHYWuC2oKIC53ox3v30XjjYH1sxe\nA5PFhEh5JMaNHIdjZ4/5dhPeNsghhxdezJs8T3R3M5qNMGYZ8ZuDv8GhtYegfV2LrOQsv6K8mBEx\niIuOQ3VTtV+XOwbm13woUh6J1Imp+LLhS4weMRrNbb4EEj7ml/N+2WXtQKAeEi8k5Ncj0OW0rnQd\nti/djuvyr/PbaUjhXfOkxZLDAWlv8uTtyai7VIdIeSRa2lvg8riQGJsIwPcAwWVbuLBhTxMIAqvr\nryaocpy4agjMFDLoDIh4MUJ09+tulg1XzXW0OZAYmwiNSoPlSctx4uIJZMRn4GTTSay7ZR2MZiPO\nXT4nenmr89VCt0qn0eHQ2kMAfHUWKqUK5nqzkFl3tDng9rihjFSKXUYwMUO5TI4GRwMYGC63X/Z7\n3cM83So4C9RDClz8TBaTX6V3yfESbF+6XcQ0AF83Pu4G49c2sFHWUCEw+YHHr7hMybSx02CymHD8\nwnEwMGQmZsJis6CxpdHvWpksJuz6cleX9SShGIrXZjhBhoMYNMYpxwEIHdwsqChAuaVcZF5JjQ+v\nHm5sacT0cdNRcrwE9Y56XGi9gPzD+aIifJRiFLSva9HsakZBRQGykrP8nvKrm6qx4YMNHRRanZ6O\nnQB5UF0ZqcT86+d3eKKPfDESC6YsgElvwrrSdeG4RHjjizcwdezUHsclBnthDNZ722QxobiqWMil\nAFfkUQw6n0Ahz6TiumLcKEQrolFQUSDG6jQ6rLtlndiVGUyGHu04Bvv6DHe6LAAkiHDDb9rnb3++\nw3vSxaa4qhhllrKg2UYalQZ2px27vtwFc70ZVrsVWclZkEGGkp+W4J/m/RPkMjnW3bIOKqUKre5W\nXy/yv+sX8aK12KhYxI+Jx93T7hYS4IDPSIwZMQZymRyz42ZjcsxkxI+JhwwyON1OKCOVfjUmHG74\neiNvwcdKg7inm0/j6NmjUOWpoM5X+zK+CjTQFGiEcVJHqztc28FeGIMF+qUZU8H4rPYzVDVWIf9w\nPgBfwafVbsX4keORMjEFKqXK77z6WqVO9B7acRADjrRwqzuBXP6Uyt+XHoNLrcdGxcJoNoKB4Y4d\ndwDwpcPmfZwn0n8B307lm/PfiIXN7rLD6XYiZkQMlJFKuDwuTBw1EedbzmPsyLG41HYJeq3et5Df\nasBLH72EuZPnYt+qfR3OK2VCSp8ylqRjuQIu4HOn8dhHm6etg4yG1G012ItoMLcj3+Hx3YPFZkFx\nVTHUMWohPc9TnF0el+gnv+iGRThx8QQcbQ6MUozC53Wfh3RHURruwELBcWLQkQYqCyoK8NJHL8Hm\ntPnFFSJkEVDIFdi7aq/Qqqp31MPlcUEh9/XV0Kq1IghusVnw5pdvYqpqqp9ibrQiGjEjYvCD+B+g\nydkkaj3GKsei2dUMhVwBp8cJGWQYEzUG7Z52vHzny3j2w2cx//r5IuuqLw197t19LzI1mSHHJm5L\nFFLvLo8LMsgwImIE5k2eh5rmmrDqL4WTgooCFFQUiI6GmYmZsDltfgaCx5GkmWa8twjgH//6cPWH\nfuq50hjQ1VTAF24GYu0kVxUxpMjOyEbj043IvzsfsVGxiJBFAADix8TD6XHi8ZLHhewIr9lo97Zj\nTNQYLL1xKRRyhZCcGD1iNE48cQKLpy1GpDwSsVGxmDt5rq/LntLXB10GX4p5s6sZHuYRsQ0GWGMe\nxAAAHLxJREFUBrvLjhZ3C57c/yScbie+qP8CsVGxWJ60vIM2U08os5SJ3UQwN9zP038O5wtOOF9w\nQiFX4PbE2+F8wYlDaw+JxkhDkeyMbOi1epj0JqFhZX7cjH2r9sGYZUROZg60aq2ozdCoNEHVfQFf\nD5RgPTs4VHsxuJCrihh0gi0MXEb9m/PfIGZEDKx2KwDA6XZCo9JAq9aKmIXdZYfb6xbuDp7WOTJy\nJDQFGtQ018DLvGh2NaPcWo4RESOw57s9GDdyHLzwiiyoSHkkbk24VUiyR8giED8mHsYso+gup1Fp\nwta7G/C54V75+BX894r/7tCrgiN1w4Tzu8NFsD4boRR4Z02aJXpnAL4U5DWz18Bg8gkxWmwWLNm1\nBBdaL/gZB652zOmrwi/RN8hwEINOqJtfr9VjU/kmNDmbRG2E1W7Fri93oexUGRQRCqhj1KLpERfy\n40+yJr1vQVv5zkqcbzkPL/OK9E6VUoWCJQXQF+txsfWi6OwHAOoYNSobKoVkiLneDMBntHrrS793\n973YV71PqPqWW8sR+aLv9vMwT8jitA9PfjhkK705weYeGMjmn+M9VYArdSu8Ql+n0WFd6TpUNVYh\nMzHTVwz691gW4HNFcoOzo3IHLDbLFTcWGZABhWIcxJCH7yxe+ugl0d8auNJWlu8QuFou720RGxUL\nZaQSDZcbRMxiwqgJqG6qRmJsIpInJOPExROIGREDR5sDre2tcLqdSIhNQGVDJQBghHwEPMyDEREj\n0OpuFUVofWkMNPKlkYiKjIJWrRXSGZmJmSJmMpyL0zqbu1TtlveKB+CX1jz9temIHxN/peYnxLGo\nB3hoqACQIACRq597KNfv9e1LtwsXh8FkEMWE/An1w5MfIkIeAafbCbvLDrlMjnpHPZSRSljtVtQ7\n6sHA0OZpE/29AWB01GgAPsPj8rjwYOqDMGYZMfLlkWEJTEdF+vS1pKKI5dZylFvLhdtquBIqHiFt\nwqTT6MROQ6/VQ6fR4YaCGzBGOQYWmwXVTdWwOW2wOW1obGn0ax9srjfDYrOI62WymPyOQwwMZDiI\nYYFOo8P86+cL14XUp56ZmCliHFLXRk1zjXDzWGwWX5vXv2fn8IVJpVShsqESIyJGCMNx2n4aAMTO\nZUflDpjrzXC6naIVbV/cI2mT0nDWcRarblolYgI5mTlXhcsl2PyDubL4zkOn8TWTstgtmIIp8DAP\nohXRUClVWDBlASaMmiDG2Zw2YSD0xXq/bCxiYCHDQQwbuEwIALHgJMYmwmKziB4NgC8AuyJlBV7/\n6+soqChA8oRk0fWutrkW1U3VogqcB92lkiE8tsGJVkTD5rQhUu67Xfq6WB1ae0i4dLh217W2AHZm\nJHndzfRx00XBJoc/HPBsrGAV6kT/Q4aDGJZIdat4bYVOo0NlfSVKjpeg5HgJGi43ICoiCgCQMCZB\n1ADoi/UiK0uv1cNis4iMqZ2VO3H3tLtxoPoAEmITYLVb0epuhc1pu7KrMer8aji48GB35hyYfeR0\nO4X0yrWAtIDztj/dhs/OfAYv84rEBwDwMi80Kg1sThs2lm/0c3NlJmYK9yQPlpPhGHiojoMYlvDA\nqEFnQE5mjnBbFK0sQnZGNjQqjZA416g0kMvkombCYrP4ubB4syb+34z4DPww4YcwZhkRFRGFqSqf\nVlRcdBxcHpfI3uILWsnxEgBd1xYEzpn/fuGZC2G/PkMV6S7j0NpDaPtNGz5Y/QEUcgUyEzMxVjkW\n7d52saPgvVF4XITHhQJ3IsTAQjsOYljDYxrSNFmuoqop8E/LlUqd6Iv1sNqtMNf7AuqzJs3yuUBi\nfamhvE5AEaFAdVM1jGajWKx4J7rAlNPe1BUM95hGODBZTJDL5KJWg/cQ12l0frEQ/juXY+H9O3qi\nikuEBzIcxLAnlL98edJyfHXuK/E3j4FUNlRCpVSJzB5pfOFk9kkAgNFshMVmwa0Jt2J/9X7otXpU\n1FbAXG9GVnIWGlsahcQGFx4MzALqas7ElV3avMnz/AxEsKpxbvi7kqIn+h8yHMSwJ9QivH3pduGe\nki5Cpd+XAgBqm2tRbi0XbhGVUiV2FTsqd4h2svWOemFctK9rxbG3L93u2+3szoJeq8fG8o2YMGoC\n9MX6LtNDr3XDEdj1jzefAq4YCOnvV0PG2dUEGQ7iqiawaRBvEmR32REXHYfZcbM71UPimU+84ZA6\nRo2N5RtFMJ27ToRI3xBuojSUCHQ/8SZbwWRXeIU4HyclKzmLpEcGATIcxDUDX6yKq4rFgsNf5wQ+\nCd/yh1v8mjztW7VPuLwaWxqh0+iGtPDgUCRYHY7UgHD477z/ezB4pT0ZjoGFDAdxTcB7QvA6jqrG\nKrg8LoweMbpDei0Akeb7zflvMHfyXFGxLIW3ceUuKy7e11V/kWsdfk2khjvwGkmFEHdU7hDuxM6u\nJ+08Bg7SqiKuKQoqCpD3cR6WTF8Ck8UEq92KnMwcAP6Lksligr5Yj7pLdWj3tiMxNhFWuxWTYyZj\nYvREVDZUIjYqFtkZ2R0ro4ex1tRAwg2D9FpJs9+EeypIADxw18JjJIHtfK9FSKuKIMLMlsNboIxU\n+uRHjL6+HYGL/LrSdSg5XoLGlka0e9sRKY9EvaMeANDc1ozL7ZcRrYiG3WWH0WyE0Wzsk+jhtUzg\nIi81GFLpdf67VA03lIwJ0f+Q4SCuCfgTaoOjAR7mgcFkgDJS2cH9BAArUlZgwqgJQg5EBhlcHhfi\nouOE0TFZTCioKBD6V/w7yDXVOTzJQKPSCJdeqNRb/ndFbUUH426ymGCuN0Or1grVXXIRDhxkOIir\nHh7fAHy9LxJjE0XRIO9lLkW66KjyVPAwDyaPnAyn2yk0r2xOG1weF4xmMhw9gbuhpFlo0oSEYIu/\n0+3scBz++eyMbD/3FLkIBwYyHMRVz4mLJ/x2FjanDRabBScunghZacyNjYd54GhzYPzI8VDHqBEh\njxCLFNdR4nUbRO8I5naSEqq9rHQ8MbCQ4SCuerYv3Y7p46YjOyMbypeUKF5Z3OVio1VrhYjijsod\nQpG1uKpYPDVz7SSeUSXtr0E7D3960l428PM8q8pis+Cb89/g3OVzQhBRna+GOkZN/TgGGDIcxDVB\ncVUxsjOyMUoxqlsLTGCGlTHLCMBXXR4ILwSk/hCh6WpXIf1cdz5vMBlgNBuDJjcQ/Q+p4xLXFD9N\n+2mPxyRPSBa/S9NvNSoNYqNixU6DP0ETwenOtaFdw/CAdhzEVYs0KC5VUe1p858NCzb4/c3dKBqV\nBnaXXfTOzkrOooUvBPya9caFF1jZbzQbYa43w2q3ijiTMlKJDQs20PUfKFg3cLvdTKvVsuXLlzPG\nGHvrrbdYSkoKk8vl7OjRo0HHnD59mul0OpaSksJSU1NZYWGheO+zzz5j8+bNY1qtls2dO5cdOXKk\nw/huTo0gukXmG5n9ctw1RWtYTllOvxz7aiKnLCfoddr26bawHJP+Da4wEGtnt3YchYWFSElJwaVL\nlwAAaWlpKCoqwmOPPRZyjEKhwLZt26DVauFwODBnzhzcfffdSE5OxtNPP41NmzZh8eLF2Lt3L55+\n+mmUlZWFww4SxIDC4xtEcAK1vzhS3bCe9tAgaZHBp0vDUVtbi9LSUjz//PPYunUrACA5ObmLUYBa\nrYZarQYAxMTEYObMmThz5gySk5Nx3XXXwW63AwBsNhuuv/76oMcwGAzid51OB51O1+X3EkQwpI2e\nwgllT4VGmhnF274Cfb9m3HColCqY683XfPGfyWSCyWQa0O/sUqvqH//xH/Hcc8+hubkZ+fn5eO+9\n98R7CxcuxKuvvoqbb7650y+xWCzIzMzE119/jZiYGFitVixYsAAymQxerxeffvopEhIS/CdGWlUE\nMey5d/e9KFpZ5CcZYq43+8WepDGi7uw+gmmBkT7YFQZdq6qkpASTJk1Cenp6ry2aw+HAihUrUFhY\niJiYGADAww8/jNdeew333nsv3n77baxduxYHDhzo1fEJghi6lFmuuKClVePZGdlCKqQ7HfxC9eq4\n1nYXQ4VODcfhw4exZ88elJaWwul0orm5GatXr8bOnTu7dfD29nbcd999WLVqFbKyrrgKjhw5gg8+\n+AAAsGLFCjzyyCN9OAWCIIY6gVLqQPfScwPHcwJ3F2Q8BpZODUdubi5yc3MBAOXl5cjPz+9gNEJt\niRhjePjhh5GSkoLsbP/t5/Tp01FeXo7MzEwcPHgQSUlJfTkHgiCGEPfuvlfsNOwuO2JyY+D2uqFR\nafDdhe/E5yw2S9hiT2Q4BpYe1XHIZDIAQFFREZ544gk0NjZi2bJlSE9Px969e1FXV4dHH30U77//\nPj755BPs2rULN910E9LT0wEAmzdvxpIlS/D73/8ev/zlL+FyuTBy5Ej8/ve/D/+ZEQQxKBStLBK/\nq/JUsG2wib+l3fy4lAiXQ+/u4k9GYvChRk4EQfQbgYaDB7F5zIIC2uFnINZOkhwhCKLfWKhZ6Pc3\n3y0YzcaBnwwRNshwEATRb/xLxr/4/c0Nh8VmIZfTMIa0qgiC6DekVd7cPcUl6EPpV1Fl+NCHdhwE\nQQwowdr1SiGF4aEP7TgIgggrnRXrGXQGUUVOgfHhC2VVEQTRb+iMOpj0pg7GZHbcbCFDzwsDpe/n\nZOb4xlNleI8ZdMkRgiCIcBBoALr6m3YjQxsyHARBhJVAVVypwKG0NSwxfCHDQRBEWOmrrhQZlaEP\nZVURBDGg9NRwUJbV0IMMB0EQYaGgoqDDa93dPawrXdfhNW4wpK4vYmhAhoMgiLDAmzNJ6Wz3IDUM\nJcdLOoztrQQ70f+Q4SAIYsAIZTiCfY43etpYvvFK0ycyIEMCCo4TBNFrCioK/NrA6ow6AN1rA1v6\nfSmMZiNsThvsLjs0BRo43U5MiZ2CpTcuRbm1HACQGJvoO7akERQF0AcXKgAkCCIs8GK/QIIV/wFA\nZUMlMhMzUW4tR7QiGr/+4a/9MrK4BDtVmvcMKgAkCGLYE6z6mxsE/l+j2RjSMAS2nCUGH4pxEAQR\nFvrSBnZ50nIA/vEOabEguaaGFuSqIgii35FKqu+o3IGczBxYbBbotXoA6OCeInoPuaoIgrgqkO4a\nNCpNp8aB+nEMfchwEAQxqAQGz6XNnciADE3IcBAEMaAEGoNgBoKMxtCGYhwEQQw6fNdhsphQbi3v\n0K+D6D4DsXaS4SAIYsggNSDBakKIrqHgOEEQ1wyBu45gfTyIoQHVcRAEMaTQqDSDPQWiC8hVRRDE\nkCOUfAnRNQOxdtKOgyCIIQftOoY2ZDgIghhy8Ipy6skxNOmW4fB4PEhPT8c999wDAHj77beRmpqK\niIgIHDt2LOiYmpoaLFy4EKmpqZg1axZee+01v/f/7d/+DTNnzsSsWbPwzDPP9PE0CIK4muDBcDIc\nQ5NuZVUVFhYiJSUFly5dAgCkpaWhqKgIjz32WMgxCoUC27Ztg1arhcPhwJw5c3DXXXdh5syZKCsr\nw549e/Dll19CoVDg/Pnz4TkbgiAIot/p0nDU1taitLQUzz//PLZu3QoASE5O7vLAarUaarUaABAT\nE4OZM2eirq4OM2fOxL//+7/j2WefhUKhAABMnDixL+dAEMRVRDAJEgCisRNAKbqDTZeG48knn8SW\nLVvQ3Nzc6y+xWCz44osvMH/+fADA999/j48++gjPPfcclEol8vPzMXfu3A7jDAaD+F2n00Gn0/V6\nDgRBDA8CjQIXRCTl3OCYTCaYTKYB/c5ODUdJSQkmTZqE9PT0Xk/M4XBgxYoVKCwsRExMDADA7Xaj\nqakJFRUV+Pzzz3H//ffj5MmTHcZKDQdBENcuBRUFgz2FIUvgQ/XGjRv7/Ts7DY4fPnwYe/bswdSp\nU/Hggw/i4MGDWL16dbcP3t7ejvvuuw+rVq1CVtaVJi/x8fH4yU9+AgCYN28e5HI5Lly40MtTIAji\naoXvPIqrisk1NYTo1HDk5uaipqYGp06dwu7du3HHHXdg586dfp8JVWjCGMPDDz+MlJQUZGf7N63P\nysrCwYMHAQDHjx9HW1sbxo8f35fzIAjiKocMx9ChR3UcMpkMAFBUVISEhARUVFRg2bJl+NGPfgQA\nqKurw7JlywAAn3zyCXbt2oWysjKkp6cjPT0de/fuBQCsXbsWJ0+eRFpaGh588MEOxoggiGuPwHTb\ngooC6Iv10Bl1KLeWQ2fUQWfUkdtqCECSIwRBDAqBnf6CBb/5ayRB0n1IHZcgiKuWUC1iA9NxAcBi\ns1BL2SEEGQ6CIAYEvvBLDUAwIwH4p+QadAaolCoyGkMIMhwEQYQVbgikC31BRQFsThsAn/tJp9F1\nMBQAQtZpZGdkB32dGBzIcBAEEVaCGQ6eTst/RFGfxFAE06KiXcbQhAwHQRBhxWKzCFl07oqqaqxC\nubVcdPez2Cwd4hbBjAQZjqEJZVURBNFnTBYTjGYjLDaL0JSaNnYaHG0OqGPUqGyoRGJsIjQqDZSR\nSuxbtY+C3f3EQKydZDgIgggrOqO/OyrwNdKc6l8oHZcgiGGBNDuK7zh4EDzQFUW7jOEP7TgIgggr\n+mI99Fp9h6wqyowaGKjnOEEQww6NStNhV0FG4+qCDAdBEGGFXFFXP+SqIgiCuIogVxVBEAQx5CDD\nQRAEQfQIMhwEQRBEjyDDQRAEQfQIMhwEQRBEjyDDQRAEQfQIMhwEQRBEjyDDQRAEQfQIMhwEQRBE\njyDDQRAEQfQIMhwEQRBEjyDDQRAEQfQIMhwEQRBEjyDDQRAEQfQIMhwEQRBEjyDD0U+YTKbBnkKf\nGM7zH85zB2j+g81wn/9A0C3D4fF4kJ6ejnvuuQcA8PbbbyM1NRURERE4duxY0DE1NTVYuHAhUlNT\nMWvWLLz22msdPvPqq69CLpfj4sWLfTiFoclw/59vOM9/OM8doPkPNsN9/gNBtwxHYWEhUlJSIJPJ\nAABpaWkoKirC7bffHnKMQqHAtm3b8PXXX6OiogK/+93v8O2334r3a2pqcODAASQmJvbxFAiCIIiB\npEvDUVtbi9LSUjzyyCOiHWFycjKSkpI6HadWq6HVagEAMTExmDlzJurq6sT7v/rVr/Db3/62L3Mn\nCIIgBgPWBStWrGDHjh1jJpOJLV++3O89nU7Hjh492tUh2KlTp9iUKVPYpUuXGGOMFRcXs+zsbMYY\nYxqNhl24cKHDGAD0Qz/0Qz/004uf/iYSnVBSUoJJkyYhPT29134/h8OBFStWoLCwEDExMWhpaUFu\nbi4OHDggPsOCNFYP9hpBEAQx+HTqqjp8+DD27NmDqVOn4sEHH8TBgwexevXqbh+8vb0d9913H1at\nWoWsrCwAQHV1NSwWC2bPno2pU6eitrYWc+bMwblz5/p2JgRBEMSAIGPdfLQvLy9Hfn4+3nvvPfHa\nwoULkZ+fjzlz5nT4PGMMa9aswfjx47Ft27aQx506dSqOHj2KcePG9WL6BEEQxEDTozoOnlVVVFSE\nhIQEVFRUYNmyZfjRj34EAKirq8OyZcsAAJ988gl27dqFsrIypKenIz09Hfv27Qt5TIIgCGKY0O9R\nlAAeeOABptVqmVarZRqNhmm1WsYYY5999pl4PS0tje3evTvo+AsXLrBFixaxG2+8kd11112sqalJ\nvJebm8umT5/OZsyYwfbv3z+g8/+///s/NmfOHJaWlsbmzJnDDh48GHS82WxmGRkZLC0tjd1zzz2s\nubmZMcZYa2srW7lyJUtLS2MzZ85kmzdvHjZzZ4yxyspKlpGRwVJTU1laWhpzOp3Dav6MMWa1Wll0\ndDTLz88P+9z7c/7dHT9U58/Y8Lh3P/vsMzZv3jym1WrZ3Llz2ZEjRxhjw+PeDTV3xnp37w644ZDy\n1FNPsU2bNjHGGGtpaWEej4cxxtjZs2fZ+PHjmdvt7jBm/fr17JVXXmGMMZaXl8eeeeYZxhhjX3/9\nNZs9ezZra2tjp06dYtOmTRPHG4j5f/HFF+zs2bOMMca++uordv311wcdM3fuXPbRRx8xxhj705/+\nxH7zm98wxhh744032MqVKxljvmuh0WiY1WodFnNvb29nN910E/vyyy8ZY4xdvHhxWF17zn333cfu\nv//+fjMcUsI5/+6ODyfhnP9wuXczMzPZvn37GGOMlZaWMp1OxxgbHvduqLn39t4dNMPh9XpZQkIC\nO3HiRIf3Tp48yW644Yag42bMmMHq6+sZYz4DM2PGDMaY74klLy9PfG7x4sXs008/7YeZ++hs/l6v\nl40bN461tbV1eC82Nlb8fvr0aZaSksIYY2zfvn3snnvuYW63m50/f54lJSX57aaG8tzff/99tmrV\nqn6ZazDCPX/GGCsqKmLr169nBoOh3w1Hf8y/O+PDRbjnP1zu3ZUrV7L/+Z//YYwx9uc//5n97Gc/\nY4wNj3s31Nx7e+8OmlbVoUOHEBcXh2nTponXjhw5gtTUVKSmpmLr1q1BxzU0NCAuLg4AEBcXh4aG\nBgC++Ep8fLz4XHx8PM6cOTOg8+e8++67mDNnDhQKRYf3UlNT8b//+78AfNItNTU1AIDFixdjzJgx\nuO6666DRaLB+/XqoVKphMffjx49DJpNhyZIlmDNnDrZs2dIv8+6v+TscDvz2t7+FwWDo13lzwj3/\n7o4PF+Ge/3C5d/Py8vDUU09hypQpWL9+PXJzcwEMj3s3cO6bN28GAHz//fe9u3d7bGq6waJFi9is\nWbM6/OzZs0d85vHHH2dbt24NOv7bb79liYmJzGazdXhPpVL5/T127FjGGGPr1q1ju3btEq8//PDD\n7N133x3w+X/11Vds2rRp7OTJk0GPXVVVxe6++242Z84ctnHjRjZ+/HjGGGNvvvkm+8lPfsLcbjc7\nd+4cmzFjRshjDLW5b9myhU2dOpVduHCBtbS0sB/84Afsww8/7PHcB2v+Tz31FHvrrbcYY4zl5OT0\naccxGPPv7vihOv/hcu/eeeed7C9/+QtjjLG33nqLLVq0iDE2PO7dUHPv7b07KK6q9vZ2FhcXx86c\nORPyM3fccQf761//2uH1GTNmCJ9eXV2dcFVt3rzZLyi1ePFiVlFREeaZ+wg1/5qaGpaUlMQOHz7c\nreN89913bP78+Ywxxn7xi1+wN998U7y3du1asZiFk3DO/ZZbbmGMMbZ79262Zs0a8d6mTZvYli1b\nwjZnKf1x7W+77Tam0WiYRqNhKpWKjRs3jv3ud78L+9wZ65/r35vxvaU/5j9c7t3Ro0eL371eLxsz\nZgxjbHjcu6Hm3tt7d1AMx969e0VwhnPq1CnW3t7OGGPMYrGwhIQEZrfbO4xdv3698Idu3ry5Q3Dc\n5XKJGInX6x2w+Tc1NbGbbrqJFRUVdTr23LlzjDHGPB4Pe+ihh9gbb7zBGGOssLCQ/fznP2eMMeZw\nOFhKSgr729/+Nizm3tTUxG6++WbW0tLC2tvb2aJFi1hpaWnY595f85diMBjYq6++Grb5BtJf1787\n48NBf8x/uNy76enpzGQyMcYY++CDD9jcuXMZY8Pj3g0194sXL/bq3h0Uw6HX69l//Md/+L325ptv\nstTUVKbVatm8efPY3r17xXuPPPKI2H1cuHCB3XnnnUHTcV9++WU2bdo0NmPGDJFBMFDz37RpE4uO\njhYpc1qtlp0/f17Mn2t6FRYWsqSkJJaUlMSeffZZMd7pdLKf/exnbNasWSwlJaXfArT9MXfGGNu1\naxdLTU1ls2bNEsZ8OM2f09+Goz/m39n44TB/xob2vcvXns8//5zdcsstbPbs2SwjI4MdO3aMMTa0\n792u5s5Y7+7dbleOEwRBEARAHQAJgiCIHkKGgyAIgugRZDgIgiCIHkGGgyAIgugRZDgIgiCIHkGG\ngyAIgugR/w8E3hxzzFk3OwAAAABJRU5ErkJggg==\n" + } + ], + "prompt_number": 5 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "\n", + "plt.figure(1)\n", + "\n", + "compute_g_10 = ptp.compute_g(10)\n", + "compute_g_100 = ptp.compute_g(100)\n", + "compute_g_1000 = ptp.compute_g(1000)\n", + "plt.plot(0,compute_g,'ro')\n", + "plt.plot(0,compute_g_100,'ro')\n", + "plt.plot(0,compute_g_1000,'ro')\n", + "\n", + "plt.show()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD9CAYAAAC1DKAUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF39JREFUeJzt3X9sVXf9x/HXZbcOpRGWBdrk3i6dvXf0ttgWaFNn4jd3\nCnZW1+CopjKEbJ02RChMnIuaOLo4oE6ihRopRjFkpusfJrQZl47hvELc2jrp1NhuXpA7b8taCaFD\nwqD0+vn+MXfT0nJpufe2tZ/nI1lyzzmfc877ncPu6/b8uNdhjDECAFhn3kwXAACYGQQAAFiKAAAA\nSxEAAGApAgAALEUAAIClbhkA7e3tys3NldfrVX19/bjlb775pu6//37Nnz9fe/bsic2PRCJ64IEH\nlJ+fr2XLlmnv3r3JrRwAkBBHvOcAotGoli5dquPHj8vlcqmkpETNzc3y+XyxMefPn9fbb7+tw4cP\n66677tL27dslSQMDAxoYGFBRUZEuX76slStX6vDhw2PWBQDMnLh/AXR1dcnj8Sg7O1tpaWmqqqpS\na2vrmDGLFy9WcXGx0tLSxszPzMxUUVGRJCk9PV0+n0/nzp1LcvkAgNvljLewv79fWVlZsWm3263O\nzs4p7yQcDqu7u1ulpaVj5jscjilvCwAgJeNLHOL+BZCMN+jLly+rsrJSDQ0NSk9PH7fcGDNn/3v6\n6adnvAb6oz8b+5vLvRmTvG/viRsALpdLkUgkNh2JROR2uye98evXr2vt2rVav3691qxZc/tVAgCS\nLm4AFBcXKxQKKRwOa3h4WC0tLaqoqJhw7I2pZIxRdXW18vLytG3btuRVDABIirjXAJxOpxobG1VW\nVqZoNKrq6mr5fD41NTVJkmpqajQwMKCSkhJdunRJ8+bNU0NDg3p6evTGG2/o+eefV0FBgZYvXy5J\n2rVrlx588MHUdzVL+P3+mS4hpejvf9tc7m8u95ZMcW8DTfnOHY6kns8CABsk672TJ4EBwFIEAABY\nigAAAEsRAABgKQIAACxFAACApQgAALBU3AfBABudOHJEx/bulfPaNY3ceac+W1ur//v852e6LCDp\nCABglBNHjuilrVv17JkzsXnf++9rQgBzDaeAgFGO7d075s1fkp49c0Yv79s3QxUBqUMAAKM4r12b\ncP4dV69OcyVA6hEAwCgjd9454fzo/PnTXAmQegQAMMpna2v1vZycMfO+m5Oj1Vu2zFBFQOrwbaDA\nDU4cOaKX9+3THVevKjp/vlZv2cIFYMwqyXrvJAAA4H8MXwcNAEgIAQAAliIAAMBSBAAAWIoAAABL\nEQAAYCkCAAAsRQAAgKUIAACwFAEAAJYiAADAUgQAAFiKAAAAS90yANrb25Wbmyuv16v6+vpxy998\n803df//9mj9/vvbs2TOldQEAMyfu10FHo1EtXbpUx48fl8vlUklJiZqbm+Xz+WJjzp8/r7fffluH\nDx/WXXfdpe3bt096Xb4OGgCmblq+Drqrq0sej0fZ2dlKS0tTVVWVWltbx4xZvHixiouLlZaWNuV1\nAQAzxxlvYX9/v7KysmLTbrdbnZ2dk9rwZNfdsWNH7LXf75ff75/U9gHAFsFgUMFgMOnbjRsADofj\ntjc82XVHBwAAYLwbPxzX1dUlZbtxTwG5XC5FIpHYdCQSkdvtntSGE1kXAJB6cQOguLhYoVBI4XBY\nw8PDamlpUUVFxYRjb7wgMZV1AQDTL+4pIKfTqcbGRpWVlSkajaq6ulo+n09NTU2SpJqaGg0MDKik\npESXLl3SvHnz1NDQoJ6eHqWnp0+4LgBgdoh7G2jKd85toAAwZdNyGygAYO4iAADAUgQAAFiKAAAA\nSxEAAGApAgAALEUAAIClCAAAsBQBAACWIgAAwFIEAABYigAAAEsRAABgKQIAACxFAACApQgAALAU\nAQAAliIAAMBSBAAAWIoAAABLEQAAYCkCAAAsRQAAgKUIAACwFAEAAJYiAADAUgQAAFiKAAAAS90y\nANrb25Wbmyuv16v6+voJx9TW1srr9aqwsFDd3d2x+bt27VJ+fr4+/vGPa926dbp27VryKgcAJCRu\nAESjUW3evFnt7e3q6elRc3Ozent7x4wJBAI6ffq0QqGQDhw4oE2bNkmSwuGwfv7zn+vUqVP661//\nqmg0qhdeeCF1nQAApiRuAHR1dcnj8Sg7O1tpaWmqqqpSa2vrmDFtbW3auHGjJKm0tFRDQ0MaHBzU\nRz/6UaWlpenKlSsaGRnRlStX5HK5UtcJAGBKnPEW9vf3KysrKzbtdrvV2dl5yzH9/f1asWKFtm/f\nrnvuuUcf/vCHVVZWplWrVo3bx44dO2Kv/X6//H7/bbYCAHNTMBhUMBhM+nbjBoDD4ZjURowx4+ad\nOXNGP/nJTxQOh7Vw4UJ96Utf0q9//Ws98sgjY8aNDgAAwHg3fjiuq6tLynbjngJyuVyKRCKx6Ugk\nIrfbHXdMX1+fXC6XXn/9dX3yk5/U3XffLafTqYcfflivvvpqUooGACQubgAUFxcrFAopHA5reHhY\nLS0tqqioGDOmoqJChw4dkiR1dHRo0aJFysjI0NKlS9XR0aH33ntPxhgdP35ceXl5qesEADAlcU8B\nOZ1ONTY2qqysTNFoVNXV1fL5fGpqapIk1dTUqLy8XIFAQB6PRwsWLNDBgwclSUVFRdqwYYOKi4s1\nb948rVixQl//+tdT3xEAYFIcZqIT+NO1c4djwusHAICbS9Z7J08CA4ClCAAAsBQBAACWIgAAwFIE\nAABYigAAAEsRAABgKQIAACxFAACApQgAALAUAQAAliIAAMBSBAAAWIoAAABLEQAAYCkCAAAsRQAA\ngKUIAACwFAEAAJYiAADAUgQAAFiKAAAASxEAAGApAgAALEUAAIClCAAAsBQBAACWIgAAwFK3DID2\n9nbl5ubK6/Wqvr5+wjG1tbXyer0qLCxUd3d3bP7Q0JAqKyvl8/mUl5enjo6O5FUOAEhI3ACIRqPa\nvHmz2tvb1dPTo+bmZvX29o4ZEwgEdPr0aYVCIR04cECbNm2KLdu6davKy8vV29urv/zlL/L5fKnp\nAgAwZXEDoKurSx6PR9nZ2UpLS1NVVZVaW1vHjGlra9PGjRslSaWlpRoaGtLg4KDeffddnTx5Uo89\n9pgkyel0auHChSlqAwAwVc54C/v7+5WVlRWbdrvd6uzsvOWYvr4+3XHHHVq8eLEeffRR/fnPf9bK\nlSvV0NCgj3zkI2PW37FjR+y13++X3+9PoB0AmHuCwaCCwWDStxs3ABwOx6Q2YowZt97IyIhOnTql\nxsZGlZSUaNu2bdq9e7eeeeaZMWNHBwAAYLwbPxzX1dUlZbtxTwG5XC5FIpHYdCQSkdvtjjumr69P\nLpdLbrdbbrdbJSUlkqTKykqdOnUqKUUDABIXNwCKi4sVCoUUDoc1PDyslpYWVVRUjBlTUVGhQ4cO\nSZI6Ojq0aNEiZWRkKDMzU1lZWfr73/8uSTp+/Ljy8/NT1AYAYKringJyOp1qbGxUWVmZotGoqqur\n5fP51NTUJEmqqalReXm5AoGAPB6PFixYoIMHD8bW37dvnx555BENDw8rJydnzDIAwMxymBtP4E/n\nzh2OcdcPAADxJeu9kyeBAcBSBAAAWIoAAABLEQAAYCkCAAAsRQAAgKUIAACwFAEAAJYiAADAUgQA\nAFiKAAAASxEAAGApAgAALEUAAIClCAAAsBQBAACWIgAAwFIEAABYigAAAEsRAABgKQIAACxFAACA\npQgAALAUAQAAliIAAMBSBAAAWIoAAABLEQAAYKlbBkB7e7tyc3Pl9XpVX18/4Zja2lp5vV4VFhaq\nu7t7zLJoNKrly5froYceSk7FAICkiBsA0WhUmzdvVnt7u3p6etTc3Kze3t4xYwKBgE6fPq1QKKQD\nBw5o06ZNY5Y3NDQoLy9PDocj+dUDAG5b3ADo6uqSx+NRdna20tLSVFVVpdbW1jFj2tratHHjRklS\naWmphoaGNDg4KEnq6+tTIBDQ448/LmNMiloAANwOZ7yF/f39ysrKik273W51dnbeckx/f78yMjL0\nxBNP6LnnntOlS5duuo8dO3bEXvv9fvn9/im2AABzWzAYVDAYTPp24wbAZE/b3Pjp3hijF198UUuW\nLNHy5cvjFj46AAAA49344biuri4p2417CsjlcikSicSmI5GI3G533DF9fX1yuVx69dVX1dbWpnvv\nvVdf+cpX9Morr2jDhg1JKRoAkLi4AVBcXKxQKKRwOKzh4WG1tLSooqJizJiKigodOnRIktTR0aFF\nixYpMzNTO3fuVCQS0dmzZ/XCCy/o05/+dGwcAGDmxT0F5HQ61djYqLKyMkWjUVVXV8vn86mpqUmS\nVFNTo/LycgUCAXk8Hi1YsEAHDx6ccFvcBQQAs4vDzODtOQ6Hg7uDAGCKkvXeyZPAAGApAgAALEUA\nAIClCAAAsBQBAACWIgAAwFIEAABYigAAAEsRAABgKQIAACxFAACApQgAALAUAQAAliIAAMBSBAAA\nWIoAAABLEQAAYCkCAAAsRQAAgKUIAACwFAEAAJYiAADAUgQAAFiKAAAASxEAAGApAgAALEUAAICl\nCAAAsBQBAACWumUAtLe3Kzc3V16vV/X19ROOqa2tldfrVWFhobq7uyVJkUhEDzzwgPLz87Vs2TLt\n3bs3uZUDABISNwCi0ag2b96s9vZ29fT0qLm5Wb29vWPGBAIBnT59WqFQSAcOHNCmTZskSWlpafrx\nj3+sv/3tb+ro6NBPf/rTcesCAGZO3ADo6uqSx+NRdna20tLSVFVVpdbW1jFj2tratHHjRklSaWmp\nhoaGNDg4qMzMTBUVFUmS0tPT5fP5dO7cuRS1AQCYKme8hf39/crKyopNu91udXZ23nJMX1+fMjIy\nYvPC4bC6u7tVWlo6bh87duyIvfb7/fL7/VPtAQDmtGAwqGAwmPTtxg0Ah8MxqY0YY2663uXLl1VZ\nWamGhgalp6ePW3d0AAAAxrvxw3FdXV1Sthv3FJDL5VIkEolNRyIRud3uuGP6+vrkcrkkSdevX9fa\ntWu1fv16rVmzJikFAwCSI24AFBcXKxQKKRwOa3h4WC0tLaqoqBgzpqKiQocOHZIkdXR0aNGiRcrI\nyJAxRtXV1crLy9O2bdtS1wEA4LbEPQXkdDrV2NiosrIyRaNRVVdXy+fzqampSZJUU1Oj8vJyBQIB\neTweLViwQAcPHpQk/eEPf9Dzzz+vgoICLV++XJK0a9cuPfjggyluCQAwGQ5z4wn86dy5wzHu+gEA\nIL5kvXfyJDAAWIoAAABLEQAAYCkCAAAsRQAAgKUIAACwFAEAAJYiAADAUgQAAFgq7ldBADY6ceSI\nju3dK+e1axq58059trZW//f5z890WUDSEQDAKCeOHNFLW7fq2TNnYvO+99/XhADmGk4BAaMc27t3\nzJu/JD175oxe3rdvhioCUocAAEZxXrs24fw7rl6d5kqA1CMAgFFG7rxzwvnR+fOnuRIg9QgAYJTP\n1tbqezk5Y+Z9NydHq7dsmaGKgNTh9wCAG5w4ckQv79unO65eVXT+fK3esoULwJhVkvXeSQAAwP8Y\nfhAGAJAQAgAALEUAAIClCAAAsBQBAACWIgAAwFIEAABYigAAAEsRAABgKX4PALgBPwgDW/AXQAoF\ng8GZLiGl5mJ/H/wgzA+OHZP/97/XD44d00tbt+rEkSMzXVrSzcXj94G53Fsy3TIA2tvblZubK6/X\nq/r6+gnH1NbWyuv1qrCwUN3d3VNady6b6/8I52J/o38QJvjfeXP1B2Hm4vH7wFzuLZniBkA0GtXm\nzZvV3t6unp4eNTc3q7e3d8yYQCCg06dPKxQK6cCBA9q0adOk1wVmG34QBjaJGwBdXV3yeDzKzs5W\nWlqaqqqq1NraOmZMW1ubNm7cKEkqLS3V0NCQBgYGJrUuMNvwgzCwSdyLwP39/crKyopNu91udXZ2\n3nJMf3+/zp07d8t1pfe/1nQuq6urm+kSUmou9vfsqNex7l56ST+Yg/9W5+Lx+8Bc7i1Z4gbAZN+c\nb/d7qfktAACYOXEDwOVyKRKJxKYjkYjcbnfcMX19fXK73bp+/fot1wUAzJy41wCKi4sVCoUUDoc1\nPDyslpYWVVRUjBlTUVGhQ4cOSZI6Ojq0aNEiZWRkTGpdAMDMifsXgNPpVGNjo8rKyhSNRlVdXS2f\nz6empiZJUk1NjcrLyxUIBOTxeLRgwQIdPHgw7roAgFnCpNiFCxfMqlWrjNfrNatXrzYXL16ccNzR\no0fN0qVLjcfjMbt37x6zbO/evSY3N9fk5+ebb3/726kueUqS0Z8xxvzoRz8yDofDXLhwIdUlT0mi\n/X3rW98yubm5pqCgwHzxi180Q0ND01V6XLc6HsYYs2XLFuPxeExBQYE5derUlNadabfb3z//+U/j\n9/tNXl6eyc/PNw0NDdNZ9qQkcuyMMWZkZMQUFRWZL3zhC9NR7pQl0t/FixfN2rVrTW5urvH5fOa1\n116Lu6+UB8CTTz5p6uvrjTHG7N692zz11FPjxoyMjJicnBxz9uxZMzw8bAoLC01PT48xxphXXnnF\nrFq1ygwPDxtjjPnXv/6V6pKnJNH+jHn/f7qysjKTnZ096wIg0f6OHTtmotGoMcaYp556asL1p9ut\njocxxhw5csR87nOfM8YY09HRYUpLSye97kxLpL933nnHdHd3G2OM+fe//23uu+++WdVfIr19YM+e\nPWbdunXmoYcemra6JyvR/jZs2GB+8YtfGGOMuX79+i0/cKX8qyBGPyewceNGHT58eNyYeM8M/Oxn\nP9N3vvMdpaWlSZIWL16c6pKnJNH+JOmb3/ymfvjDH05bzVORaH+rV6/WvHnv/zMrLS1VX1/f9BV/\nE3P9+Zbb7W9wcFCZmZkqKiqSJKWnp8vn8+ncuXPT3sPNJNKb9P5NKoFAQI8//visvAsxkf7effdd\nnTx5Uo899pik90/DL1y4MO7+Uh4Ag4ODysjIkCRlZGTEDsRoN3uWQJJCoZBOnDihT3ziE/L7/Xr9\n9ddTXfKUJNpfa2ur3G63CgoKpqfgKUq0v9F++ctfqry8PHXFTtJk6p3K8y0T9TqTbre/G8M5HA6r\nu7tbpaWlqS14ChI5dpL0xBNP6Lnnnot9KJltEjl2Z8+e1eLFi/Xoo49qxYoV+trXvqYrV67E3V9S\nvg109erVGhgYGDf/2WefHTPtcDgmfLYg3vMGIyMjunjxojo6OvTHP/5RX/7yl/WPf/wj8aKnIFX9\nvffee9q5c6defvnl2LyZ+FSSyuM3elsf+tCHtG7dutsvNElS/XzLTLvd/kavd/nyZVVWVqqhoUHp\n6elJrS8Rt9ubMUYvvviilixZouXLl8/a7wpK5NiNjIzo1KlTamxsVElJibZt26bdu3frmWeeuel2\nkhIAo9/AbpSRkaGBgQFlZmbqnXfe0ZIlS8aNife8gdvt1sMPPyxJKikp0bx583ThwgXdfffdySh9\nUlLV35kzZxQOh1VYWCjp/T9PV65cqa6urgm3kyqpPH6S9Ktf/UqBQEC//e1vk1v4bZrrz7fcbn8u\nl0uSdP36da1du1br16/XmjVrpqfoSUqkt9/85jdqa2tTIBDQ1atXdenSJW3YsCF2G/tskEh/xhi5\n3W6VlJRIkiorK7V79+74O0zStYubevLJJ2NXsnft2jXhRcDr16+bj33sY+bs2bPm2rVrYy587N+/\n33z/+983xhjz1ltvmaysrFSXPCWJ9jfabL0InEh/R48eNXl5eeb8+fPTWnc8kzkeoy+0vfbaa7EL\nbZM9ljMpkf7+85//mK9+9atm27Zt0173ZCTS22jBYHBW3gWUaH+f+tSnzFtvvWWMMebpp5++5V2T\n03Ib6Gc+85lxtxH29/eb8vLy2LhAIGDuu+8+k5OTY3bu3BmbPzw8bNavX2+WLVtmVqxYYX73u9+l\nuuQpSbS/0e69995ZFwCJ9ufxeMw999xjioqKTFFRkdm0adO09zCRierdv3+/2b9/f2zMN77xDZOT\nk2MKCgrMn/70p7jrzja329/JkyeNw+EwhYWFsWN29OjRGenhZhI5dh8IBoOz8i4gYxLr74033jDF\nxcWTvu3aYcz/6IlOAEBCZuelcABAyhEAAGApAgAALEUAAIClCAAAsBQBAACW+n/PDSsShdV0FQAA\nAABJRU5ErkJggg==\n" + } + ], + "prompt_number": 28 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "from point_pattern import PointPattern\n", + "from point import Point\n", + "from utils import *\n", + "%pylab inline\n", + "import matplotlib.pyplot as plt\n", + "import pysal as ps\n", + "\n", + "shapefile = ps.open(ps.examples.get_path('new_haven_merged.shp'))\n", + "dbf = ps.open(ps.examples.get_path('new_haven_merged.dbf'))\n", + "\n", + "point_list=[]\n", + "\n", + "for geometry, attributes in zip(shapefile, dbf):\n", + " point_list.append(Point(geometry[0],geometry[1],attributes[0]))\n", + "\n", + "ptp=PointPattern(point_list)\n", + "\n", + "#Run a few tests to explore the dataset.\n", + "nn = ptp.average_nearest_neighbor_distance_KDTree()\n", + "\n", + "smallest, largest = ptp.critical_points(ptp.permutations(n=1000))\n", + "plt.figure(1)\n", + "\n", + "\n", + "plt.plot(0,nn,'ro')\n", + "plt.plot([-1,1],[smallest,smallest],'b+-')\n", + "plt.plot([-1,1],[largest,largest],'b+-')\n", + "\n", + "plt.axis([-1,1,0,largest+0.01])\n", + "plt.show()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].\n", + "For more information, type 'help(pylab)'.\n" + ] + }, + { + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD5CAYAAADFqlkBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFF1JREFUeJzt3H9I2/nhx/HX577KHbMbucGMJZHZmThta6O3HGGDHjl6\n1tbR4Ljj5sGB9CyIW+fKldGOL+3i9u2d/iHHtW6cHSLKoNg/1irUZtexy8pui4HV/jO7mZbKYk5l\n1Anrxs5W3t8/bgtN1Y+pWn+0zwcEkvj+fPL+hI95xo/5xDLGGAEAnmrPrPcEAADrjxgAAIgBAIAY\nAABEDAAAknLWewKLsSxrvacAAJvOcj8guqH/MjDGcFmFy49+9KN1n8OTdOH55PncqJeV2NAxAACs\nDWIAACAGT4NgMLjeU3ii8HyuLp7PjcEyKz3Q9JhYlrXiY2AA8DRZyesmfxkAAIgBAIAYAABEDAAA\nIgYAAGURg0gkotLSUnm9XrW1tS04prm5WV6vVz6fT8PDw5KkZDKpl19+WTt27NDOnTt1+vTp9Phw\nOCy3263KykpVVlYqEoms0uYAAJbD9ruJ5ubmdPjwYf3617+Wy+XSiy++qFAopLKysvSYwcFB3bx5\nU4lEQkNDQ2pqalIsFlNubq7ee+89VVRU6O7du/ra176mvXv3qrS0VJZl6e2339bbb7/92DcQALA0\n278M4vG4PB6PioqKlJubq7q6OvX392eMGRgYUH19vSQpEAhoZmZGU1NTKigoUEVFhSRpy5YtKisr\nUyqVSi/HOQQAsHHY/mWQSqVUWFiYvu12uzU0NLTkmPHxcTmdzvR9Y2NjGh4eViAQSN935swZ9fb2\nyu/3q729XQ6HY97jh8Ph9PVgMMiZigDwgGg0qmg0uirrso1Btl8j/fC7/AeXu3v3rl577TW9//77\n2rJliySpqalJJ0+elCSdOHFCR48eVVdX17z1PhgDAECmh98kt7S0LHtdtoeJXC6Xkslk+nYymZTb\n7bYdMz4+LpfLJUm6d++eXn31Vb355puqra1Nj8nPz5dlWbIsS4cOHVI8Hl/2BgAAVs42Bn6/X4lE\nQmNjY5qdnVVfX59CoVDGmFAopN7eXklSLBaTw+GQ0+mUMUYNDQ3avn27jhw5krHMxMRE+vqFCxdU\nXl6+WtsDAFgG28NEOTk56ujoUHV1tebm5tTQ0KCysjJ1dnZKkhobG1VTU6PBwUF5PB7l5eWpu7tb\nkvTxxx/rF7/4hXbt2qXKykpJ0rvvvqt9+/bp2LFjun79uizL0rZt29LrAwCsD761FACeEHxrKQBg\nRYgBAIAYAACIAQBAxAAAIGIAABAxAACIGAAARAwAACIGAAARAwCAiAEAQMQAACBiAAAQMQAAiBgA\nAEQMAAAiBgAAEQMAgIgBAEDEAAAgYgAAEDEAAIgYAABEDAAAIgYAABEDAICIAQBAxAAAIGIAABAx\nAACIGAAARAwAAMoiBpFIRKWlpfJ6vWpra1twTHNzs7xer3w+n4aHhyVJyWRSL7/8snbs2KGdO3fq\n9OnT6fHT09OqqqpSSUmJ9u7dq5mZmVXaHADActjGYG5uTocPH1YkEtHIyIjOnTunGzduZIwZHBzU\nzZs3lUgkdPbsWTU1NUmScnNz9d577+lPf/qTYrGYfvrTn+rPf/6zJKm1tVVVVVUaHR3Vnj171Nra\n+pg2DwCQDdsYxONxeTweFRUVKTc3V3V1derv788YMzAwoPr6eklSIBDQzMyMpqamVFBQoIqKCknS\nli1bVFZWplQqNW+Z+vp6Xbx4cdU3DACQvRy7H6ZSKRUWFqZvu91uDQ0NLTlmfHxcTqczfd/Y2JiG\nh4cVCAQkSVNTU+mfO51OTU1NLfj44XA4fT0YDCoYDGa3VQDwFIhGo4pGo6uyLtsYWJaV1UqMMYsu\nd/fuXb322mt6//33tWXLlgUfY7HHeTAGAIBMD79JbmlpWfa6bA8TuVwuJZPJ9O1kMim32207Znx8\nXC6XS5J07949vfrqq3rzzTdVW1ubHuN0OjU5OSlJmpiYUH5+/rI3AACwcrYx8Pv9SiQSGhsb0+zs\nrPr6+hQKhTLGhEIh9fb2SpJisZgcDoecTqeMMWpoaND27dt15MiRecv09PRIknp6ejJCAQBYe5Z5\n+BjPQy5fvqwjR45obm5ODQ0N+uEPf6jOzk5JUmNjoySlP3GUl5en7u5uvfDCC/rd736nl156Sbt2\n7UofBnr33Xe1b98+TU9P6/XXX9df//pXFRUV6fz583I4HJkTs6x5h58AAItbyevmkjFYL8QAAB7N\nSl43OQMZAEAMAADEAAAgYgAAEDEAAIgYAABEDAAAIgYAABEDAICIAQBAxAAAIGIAABAxAACIGAAA\nRAwAACIGAAARAwCAiAEAQMQAACBiAAAQMQAASMpZ7wnYsaz1ngEAPB02dAyMWe8ZAMDmsZI30Bwm\nAgAQAwAAMQAAiBgAAEQMAAAiBgAAEQMAgIgBAEDEAAAgYgAAUBYxiEQiKi0tldfrVVtb24Jjmpub\n5fV65fP5NDw8nL7/rbfektPpVHl5ecb4cDgst9utyspKVVZWKhKJrHAzAAArYmzcv3/fFBcXm9u3\nb5vZ2Vnj8/nMyMhIxphLly6Z/fv3G2OMicViJhAIpH929epVc+3aNbNz586MZcLhsGlvb7d7aCPJ\nfPbtRFy4cOHCJbuLbF9X7dh+UV08HpfH41FRUZEkqa6uTv39/SorK0uPGRgYUH19vSQpEAhoZmZG\nk5OTKigo0O7duzU2NrZYhLIIVVY9AwBoZV9UZxuDVCqlwsLC9G23262hoaElx6RSKRUUFNg+8Jkz\nZ9Tb2yu/36/29nY5HI55Y8LhcPp6MBhUMBi0XScAPE2i0aii0eiqrMs2BlaWmXn4Xf5SyzU1Nenk\nyZOSpBMnTujo0aPq6uqaN+7BGAAAMj38JrmlpWXZ67L9B7LL5VIymUzfTiaTcrvdtmPGx8flcrls\nHzQ/P1+WZcmyLB06dEjxeHw5cwcArBLbGPj9fiUSCY2NjWl2dlZ9fX0KhUIZY0KhkHp7eyVJsVhM\nDodDTqfT9kEnJibS1y9cuDDv00YAgLVle5goJydHHR0dqq6u1tzcnBoaGlRWVqbOzk5JUmNjo2pq\najQ4OCiPx6O8vDx1d3enl3/jjTf029/+Vnfu3FFhYaF+/OMf6+DBgzp27JiuX78uy7K0bdu29PoA\nAOvDMtl8rGcdWJaV1SeOAACfWcnrJmcgAwCIAQCAGAAARAwAACIGAAARAwCAiAEAQMQAACBiAAAQ\nMQAAiBgAAEQMAAAiBgAAEQMAgIgBAEDEAAAgYgAAEDEAAIgYAABEDAAAIgYAABEDAICIAQBAxAAA\nIGIAABAxAACIGAAARAwAACIGAAARAwCAiAEAQMQAAKAsYhCJRFRaWiqv16u2trYFxzQ3N8vr9crn\n82l4eDh9/1tvvSWn06ny8vKM8dPT06qqqlJJSYn27t2rmZmZFW4GAGAlbGMwNzenw4cPKxKJaGRk\nROfOndONGzcyxgwODurmzZtKJBI6e/asmpqa0j87ePCgIpHIvPW2traqqqpKo6Oj2rNnj1pbW1dp\ncwAAy2Ebg3g8Lo/Ho6KiIuXm5qqurk79/f0ZYwYGBlRfXy9JCgQCmpmZ0eTkpCRp9+7dev755+et\n98Fl6uvrdfHixVXZGADA8uTY/TCVSqmwsDB92+12a2hoaMkxqVRKBQUFi653ampKTqdTkuR0OjU1\nNbXguHA4nL4eDAYVDAbtpgsAT5VoNKpoNLoq67KNgWVZWa3EGLOs5f47drHxD8YAAJDp4TfJLS0t\ny16X7WEil8ulZDKZvp1MJuV2u23HjI+Py+Vy2T6o0+lMH0qamJhQfn7+I08cALB6bGPg9/uVSCQ0\nNjam2dlZ9fX1KRQKZYwJhULq7e2VJMViMTkcjvQhoMWEQiH19PRIknp6elRbW7uSbQAArJBtDHJy\nctTR0aHq6mpt375d3/72t1VWVqbOzk51dnZKkmpqavSVr3xFHo9HjY2N+tnPfpZe/o033tA3vvEN\njY6OqrCwUN3d3ZKk48eP68qVKyopKdFvfvMbHT9+/DFuIgBgKZZ5+ID/BmFZ1rz/RQAAFreS103O\nQAYAEAMAADEAAIgYAABEDAAAIgYAABEDAICIAQBAxAAAIGIAABAxAACIGAAARAwAACIGAAARAwCA\niAEAQMQAACBiAAAQMQAAiBgAAEQMAAAiBgAAEQMAgIgBAEDEAAAgYgAAEDEAAIgYAABEDAAAIgYA\nABEDAICIAQBAxAAAoCxiEIlEVFpaKq/Xq7a2tgXHNDc3y+v1yufzaXh4eMllw+Gw3G63KisrVVlZ\nqUgksgqbAgBYNmPj/v37pri42Ny+fdvMzs4an89nRkZGMsZcunTJ7N+/3xhjTCwWM4FAYMllw+Gw\naW9vt3tos8TUAAAPWcnrpu1fBvF4XB6PR0VFRcrNzVVdXZ36+/szxgwMDKi+vl6SFAgENDMzo8nJ\nySWX/WzeAICNIMfuh6lUSoWFhenbbrdbQ0NDS45JpVL65JNPbJc9c+aMent75ff71d7eLofDMe/x\nw+Fw+nowGFQwGMx6wwDgSReNRhWNRldlXbYxsCwrq5U86rv8pqYmnTx5UpJ04sQJHT16VF1dXfPG\nPRgDAECmh98kt7S0LHtdtjFwuVxKJpPp28lkUm6323bM+Pi43G637t27t+iy+fn56fsPHTqkAwcO\nLHsDAAArZ/s/A7/fr0QiobGxMc3Ozqqvr0+hUChjTCgUUm9vryQpFovJ4XDI6XTaLjsxMZFe/sKF\nCyovL1/t7QIAPALbvwxycnLU0dGh6upqzc3NqaGhQWVlZers7JQkNTY2qqamRoODg/J4PMrLy1N3\nd7ftspJ07NgxXb9+XZZladu2ben1AQDWh2U26Md6LMviE0cA8AhW8rrJGcgAAGIAACAGAAARAwCA\niAEAQMQAACBiAAAQMQAAiBgAAEQMAAAiBgAAEQMAgIgBAEDEAAAgYgAAEDEAAIgYAABEDAAAIgYA\nABEDAICIAQBAxAAAIGIAABAxAACIGAAARAwAACIGAAARAwCAiAEAQMQAACBiAACQlLPeEwA2i6uX\nLunD06eV8+mnuv/ss9rb3KyXvvnN9Z4WsCqIwVMgGo0qGAyu9zQ2tauXLulX3/++Tt26paikoKT/\nvXVLkgjCCrF/bgxLHiaKRCIqLS2V1+tVW1vbgmOam5vl9Xrl8/k0PDy85LLT09OqqqpSSUmJ9u7d\nq5mZmVXYFCwmGo2u9xQ2vQ9Pn9ap/7z4R/9z36lbt3TlzJl1m9OTgv1zY7CNwdzcnA4fPqxIJKKR\nkRGdO3dON27cyBgzODiomzdvKpFI6OzZs2pqalpy2dbWVlVVVWl0dFR79uxRa2vrY9o8YHXkfPrp\ngvf/z7//vcYzAR4P2xjE43F5PB4VFRUpNzdXdXV16u/vzxgzMDCg+vp6SVIgENDMzIwmJydtl31w\nmfr6el28ePFxbBuwau4/++yC988999wazwR4PGz/Z5BKpVRYWJi+7Xa7NTQ0tOSYVCqlTz75ZNFl\np6am5HQ6JUlOp1NTU1MLPr5lWY+4OVhMS0vLek9h0zv1wPX0s/mrX+n/2E9XjP1z/dnGINsXY2NM\nVmMWWp9lWQven806AQCrw/YwkcvlUjKZTN9OJpNyu922Y8bHx+V2uxe83+VySfrsr4HJyUlJ0sTE\nhPLz81e+JQCAZbONgd/vVyKR0NjYmGZnZ9XX16dQKJQxJhQKqbe3V5IUi8XkcDjkdDptlw2FQurp\n6ZEk9fT0qLa29nFsGwAgS7aHiXJyctTR0aHq6mrNzc2poaFBZWVl6uzslCQ1NjaqpqZGg4OD8ng8\nysvLU3d3t+2yknT8+HG9/vrr6urqUlFRkc6fP/+YNxMAYMtsEOfPnzfbt283zzzzjPnjH/+46LjL\nly+br371q8bj8ZjW1tY1nOHmcufOHfPKK68Yr9drqqqqzN///vcFx335y1825eXlpqKiwrz44otr\nPMuNLZt97Xvf+57xeDxm165d5tq1a2s8w81lqefzo48+Ml/4whdMRUWFqaioMD/5yU/WYZabw8GD\nB01+fr7ZuXPnomMedd/cMDG4ceOG+ctf/mKCweCiMbh//74pLi42t2/fNrOzs8bn85mRkZE1nunm\n8IMf/MC0tbUZY4xpbW01x44dW3BcUVGRuXPnzlpObVPIZl+7dOmS2b9/vzHGmFgsZgKBwHpMdVPI\n5vn86KOPzIEDB9ZphpvL1atXzbVr1xaNwXL2zQ3zRXWlpaUqKSmxHZPNeQ/4zKOcy2H45NY8yz3H\nZrGPST/tsv3dZV/Mzu7du/X8888v+vPl7JsbJgbZWOycBsz3KOdyvPLKK/L7/fr5z3++llPc0LLZ\n1xYaMz4+vmZz3EyyeT4ty9Lvf/97+Xw+1dTUaGRkZK2n+cRYzr65pl9UV1VVlf5I6YPeeecdHThw\nYMnlOQkt02LP56lTpzJuL3YuhyR9/PHH2rp1q/72t7+pqqpKpaWl2r1792OZ72ay3HNs2EcXls3z\n8sILLyiZTOpzn/ucLl++rNraWo2Ojq7B7J5Mj7pvrmkMrly5sqLlsznv4Wli93z+91yOgoIC23M5\ntm7dKkn60pe+pG9961uKx+PEQMs/x+a/59IgUzbP5+c///n09f379+s73/mOpqen9cUvfnHN5vmk\nWM6+uSEPEy123DCb8x7wmWzO5fjXv/6lf/zjH5Kkf/7zn/rwww9VXl6+pvPcqFZyjg3my+b5nJqa\nSv/ux+NxGWMIwTIta99cnf9tr9wvf/lL43a7zXPPPWecTqfZt2+fMcaYVCplampq0uMGBwdNSUmJ\nKS4uNu+88856TXfDu3PnjtmzZ8+8j5Y++HzeunXL+Hw+4/P5zI4dO3g+H7LQvvbBBx+YDz74ID3m\nu9/9rikuLja7du2y/Ug0ln4+Ozo6zI4dO4zP5zNf//rXzR/+8If1nO6GVldXZ7Zu3Wpyc3ON2+02\nXV1dK943LWP49z0APO025GEiAMDaIgYAAGIAACAGAAARAwCAiAEAQNL/AxnsWomlDkNnAAAAAElF\nTkSuQmCC\n" + } + ], + "prompt_number": 10 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/point_pattern.py b/point_pattern.py new file mode 100644 index 0000000..a503cf0 --- /dev/null +++ b/point_pattern.py @@ -0,0 +1,281 @@ +from point import Point +import math +import random +import numpy as np +import scipy.spatial as ss + +class PointPattern(): + + def __init__(self,points): + self.points = points[:] + + + def average_nearest_neighbor_distance(self,mark=None,points_list=None): + + """ + Given a set of points, compute the average nearest neighbor. + + Parameters + ---------- + points_list : list + A list of points + mark : str + + Returns + ------- + mean_d : float + Average nearest neighbor distance + + References + ---------- + Clark and Evan (1954 Distance to Nearest Neighbor as a + Measure of Spatial Relationships in Populations. Ecology. 35(4) + p. 445-453. + """ + mean_d = 0 + if points_list==None: + points_tmp=self.points + else: + points_tmp=points_list + + if mark!=None: + points_tmp=[ point for point in self.points if point.mark==mark] + + length=len(points_tmp) + + #print(self.points) + nearest_distances=[] + for i in range(length): + distance=[] + for j in range(length): + if i==j: + continue + else: + distance.append(self.euclidean_distance((points_tmp[i].x,points_tmp[i].y),(points_tmp[j].x,points_tmp[j].y))) + nearest_distances.append(min(distance)) + + mean_d=float(sum(nearest_distances)/len(nearest_distances)) + return mean_d + + + def average_nearest_neighbor_distance_KDTree(self,mark=None): + + """ + Given a set of points, compute the average nearest neighbor. + + Parameters + ---------- + points : list + A list of points + mark : str + + Returns + ------- + mean_d : float + Average nearest neighbor distance + + References + ---------- + Clark and Evan (1954 Distance to Nearest Neighbor as a + Measure of Spatial Relationships in Populations. Ecology. 35(4) + p. 445-453. + """ + + if mark==None: + points_tmp=self.points + else: + points_tmp=[ point for point in self.points if point.mark==mark] + + nearest_distances=[] + point_list=[(point.x,point.y) for point in points_tmp] + point_stack = np.vstack(point_list) + kdtree = ss.KDTree(point_stack) + + for p in point_stack: + nearest_neighbor_distance,nearest_neighbor = kdtree.query(p, k=2) + nearest_distances.append(nearest_neighbor_distance[1]) + + nn_distances = np.array(nearest_distances) + + return np.mean(nn_distances) + + def average_nearest_neighbor_distance_np(self,mark=None): + + """ + Given a set of points, compute the average nearest neighbor. + + Parameters + ---------- + points : list + A list of points + mark : str + + Returns + ------- + mean_d : float + Average nearest neighbor distance + + References + ---------- + Clark and Evan (1954 Distance to Nearest Neighbor as a + Measure of Spatial Relationships in Populations. Ecology. 35(4) + p. 445-453. + """ + + if mark==None: + points_tmp=self.points + else: + points_tmp=[ point for point in self.points if point.mark==mark] + + nearest_distances=[] + point_list=[(point.x,point.y) for point in points_tmp] + ndarray = np.array(point_list) + + for i in range(len(ndarray)): + distance=[] + for j in range(len(ndarray)): + if i==j: + continue + else: + distance.append(ss.distance.euclidean(ndarray[i], ndarray[j])) + nearest_distances.append(np.min(distance)) + + return np.mean(nearest_distances) + + def number_of_coincident_points(self): + + length=len(self.points) + number_of_coincident_points=0 + for i in range(length): + for j in range(length): + if i==j: + continue + else: + if self.check_coincident((self.points[i].x,self.points[i].y),(self.points[j].x,self.points[j].y)): + number_of_coincident_points+=1 + return number_of_coincident_points + + def list_marks(self): + + points_marks=set() + for point in self.points: + if point.mark!=None: + points_marks.add(point.mark) + + return list(points_marks) + + def list_subset_of_points(self,mark): + + return [point for point in self.points if point.mark==mark] + + def create_random_marked_points(self, n=None,domain_min=0,domain_max=1,marks=[]): + + if n == None: + n=len(self.points) + + randoms_list=np.random.uniform(domain_min, domain_max, (n,2)) + points_list=[] + for i in range(n): + if len(marks)!=0: + points_list.append(Point(randoms_list[i][0], randoms_list[i][1],random.choice(marks))) + else: + points_list.append(Point(randoms_list[i][0], randoms_list[i][1])) + return points_list + + def create_realizations(self, k): + + return self.permutations(k) + + def permutations(self,p=99, n=100): + """ + Return the mean nearest neighbor distance of p permutations. + + Parameters + ---------- + p : integer + n : integer + + Returns + ------- + permutations : list + the mean nearest neighbor distance list. + + """ + permutation_list=[] + for i in range(p): + permutation_list.append(self.average_nearest_neighbor_distance(points_list=self.create_random_marked_points(n))) + return permutation_list + + def critical_points(self,permutations): + """ + Return the mean nearest neighbor distance of p permutations. + + Parameters + ---------- + permutations : list + the mean nearest neighbor distance list. + Returns + ------- + smallest : float + largest : float + + """ + + return min(permutations),max(permutations) + + def compute_g(self,nsteps): + + ds = np.linspace(0, 1, nsteps) + min_list=[] + for i_index,di in enumerate(ds): + tmp_list=[] + for j_index,dj in enumerate(ds): + if i_index != j_index: + tmp_list.append(np.abs(di-dj)) + + min_list.append(np.min(tmp_list)) + + return np.mean(min_list) + + + + def check_coincident(self,a, b): + """ + Check whether two points are coincident + Parameters + ---------- + a : tuple + A point in the form (x,y) + + b : tuple + A point in the form (x,y) + + Returns + ------- + equal : bool + Whether the points are equal + """ + return a == b + + def euclidean_distance(self,a, b): + """ + Compute the Euclidean distance between two points + + Parameters + ---------- + a : tuple + A point in the form (x,y) + + b : tuple + A point in the form (x,y) + + Returns + ------- + + distance : float + The Euclidean distance between the two points + """ + distance = math.sqrt((a[0] - b[0])**2 + (a[1] - b[1])**2) + return distance + + diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..1a389b5 --- /dev/null +++ b/utils.py @@ -0,0 +1,246 @@ +import math +import random +from analytics import average_nearest_neighbor_distance +from point import Point + +def generate_random_points(n): + """ + Return n random points. + + Parameters + ---------- + n : integer + + Returns + ------- + points_list : list + The random points list + + """ + return [Point(random.uniform(0,1), random.uniform(0,1)) for i in range(n)] + + +def create_random_marked_points(n, marks=[]): + + points_list=[] + for i in range(n): + points_list.append(Point(random.uniform(0,1), random.uniform(0,1),random.choice(marks))) + + return points_list + + + +def permutation(p=99, n=100): + """ + Return the mean nearest neighbor distance of p permutations. + + Parameters + ---------- + p : integer + n : integer + + Returns + ------- + permutations : list + the mean nearest neighbor distance list. + + """ + permutation_list=[] + for i in range(p): + permutation_list.append(average_nearest_neighbor_distance(generate_random_points(n))) + return permutation_list + + +def permutation_mark(p=99, n=100 ,marks=None,mark=None): + """ + Return the mean nearest neighbor distance of p permutations. + + Parameters + ---------- + p : integer + n : integer + marks : list + + Returns + ------- + permutations : list + the mean nearest neighbor distance list. + + """ + permutation_list=[] + for i in range(p): + permutation_list.append(average_nearest_neighbor_distance(create_random_marked_points(n,marks),mark)) + return permutation_list + +def critical_points(permutations): + """ + Return the mean nearest neighbor distance of p permutations. + + Parameters + ---------- + permutations : list + the mean nearest neighbor distance list. + Returns + ------- + smallest : float + largest : float + + """ + + return min(permutations),max(permutations) + + +def is_observed_distance_significant(smallest,largest,observed_distance): + """ + Returns True is the observed distance is significant. + + Parameters + ---------- + smallest : float + largest : float + + Returns + ------- + bool + + """ + + if observed_distance>=smallest and observed_distance<=largest: + return False + else: + return True + + + + +def manhattan_distance(a, b): + """ + Compute the Manhattan distance between two points + + Parameters + ---------- + a : tuple + A point in the form (x,y) + + b : tuple + A point in the form (x,y) + + Returns + ------- + distance : float + The Manhattan distance between the two points + """ + distance = abs(a[0] - b[0]) + abs(a[1] - b[1]) + return distance + + +def shift_point(point, x_shift, y_shift): + """ + Shift a point by some amount in the x and y directions + + Parameters + ---------- + point : tuple + in the form (x,y) + + x_shift : int or float + distance to shift in the x direction + + y_shift : int or float + distance to shift in the y direction + + Returns + ------- + new_x : int or float + shited x coordinate + + new_y : int or float + shifted y coordinate + + Note that the new_x new_y elements are returned as a tuple + + Example + ------- + >>> point = (0,0) + >>> shift_point(point, 1, 2) + (1,2) + """ + x = getx(point) + y = gety(point) + + x += x_shift + y += y_shift + + return x, y + + +def check_coincident(a, b): + """ + Check whether two points are coincident + Parameters + ---------- + a : tuple + A point in the form (x,y) + + b : tuple + A point in the form (x,y) + + Returns + ------- + equal : bool + Whether the points are equal + """ + return a == b + + +def check_in(point, point_list): + """ + Check whether point is in the point list + + Parameters + ---------- + point : tuple + In the form (x,y) + + point_list : list + in the form [point, point_1, point_2, ..., point_n] + """ + return point in point_list + + +def getx(point): + """ + A simple method to return the x coordinate of + an tuple in the form(x,y). We will look at + sequences in a coming lesson. + + Parameters + ---------- + point : tuple + in the form (x,y) + + Returns + ------- + : int or float + x coordinate + """ + return point[0] + + +def gety(point): + """ + A simple method to return the x coordinate of + an tuple in the form(x,y). We will look at + sequences in a coming lesson. + + Parameters + ---------- + point : tuple + in the form (x,y) + + Returns + ------- + : int or float + y coordinate + """ + return point[1] diff --git a/view.py b/view.py new file mode 100644 index 0000000..083cb08 --- /dev/null +++ b/view.py @@ -0,0 +1,83 @@ +import sys +from PyQt4 import QtGui + + +class Assignment9(QtGui.QMainWindow): + + def __init__(self): + + super(Assignment9, self).__init__() + self.initUI() + + + + def initUI(self): + + #create the exit aciton + exitAction = QtGui.QAction(QtGui.QIcon('exit.png'), '&Exit', self) + exitAction.setShortcut('Ctrl+Q') + exitAction.setStatusTip('Exit application') + exitAction.triggered.connect((QtGui.qApp.quit)) + + + #set statusBar + self.statusBar() + + #create the textEdit + textEdit = QtGui.QTextEdit() + self.setCentralWidget(textEdit) + + + #self.setCentralWidget(a_botton) + + menubar = self.menuBar() + fileMenu = menubar.addMenu('&File') + fileMenu.addAction(exitAction) + + self.setGeometry(300, 300, 300, 200) + + #set the titile + self.setWindowTitle('Assignment9') + #make the window in the center + + self.center() + self.show() + + + def closeEvent(self, event): + """ + + Show a message box when the user close the window. + """ + + reply = QtGui.QMessageBox.question(self, 'Message', + "Are you sure to quit?", QtGui.QMessageBox.Yes | + QtGui.QMessageBox.No, QtGui.QMessageBox.No) + + if reply == QtGui.QMessageBox.Yes: + event.accept() + else: + event.ignore() + + def center(self): + """ + + Set the window on the center of deskp. + """ + + qr = self.frameGeometry() + cp = QtGui.QDesktopWidget().availableGeometry().center() + qr.moveCenter(cp) + self.move(qr.topLeft()) + + + +def main(): + + app = QtGui.QApplication(sys.argv) + ex = Assignment9() + sys.exit(app.exec_()) + + +if __name__ == '__main__': + main() \ No newline at end of file