{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Quadrat Based Statistical Method for Planar Point Patterns\n", "\n", "**Authors: Serge Rey , Wei Kang and Hu Shao **\n", "\n", "## Introduction\n", "\n", "In this notebook, we are going to introduce how to apply quadrat statistics to a point pattern to infer whether it comes from a CSR process.\n", "\n", "1. In [Quadrat Statistic](#Quadrat-Statistic) we introduce the concept of quadrat based method.\n", "2. We illustrate how to use the module **quadrat_statistics.py** through an example dataset **juvenile** in [Juvenile Example](#Juvenile-Example)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quadrat Statistic\n", "\n", "In the previous notebooks, we introduced the concept of Complete Spatial Randomness (CSR) process which serves as the benchmark process. Utilizing CSR properties, we can discriminate those that are not from a CSR process. Quadrat statistic is one such method. Since a CSR process has two major characteristics:\n", "1. Uniform: each location has equal probability of getting a point (where an event happens).\n", "2. Independent: location of event points are independent.\n", "\n", "We can imagine that for any point pattern, if the underlying process is a CSR process, the expected point counts inside any cell of area $|A|$ should be $\\lambda |A|$ ($\\lambda$ is the intensity which is uniform across the study area for a CSR). Thus, if we impose a $m \\times k$ rectangular tessellation over the study area (window), we can easily calculate the expected number of points inside each cell under the null of CSR. By comparing the observed point counts against the expected counts and calculate a $\\chi^2$ test statistic, we can decide whether to reject the null based on the position of the $\\chi^2$ test statistic in the sampling distribution. \n", "\n", "$$\\chi^2 = \\sum^m_{i=1} \\sum^k_{j=1} \\frac{[x_{i,j}-E(x_{i,j})]^2}{\\lambda |A_{i,j}|}$$\n", "\n", "There are two ways to construct the sampling distribution and acquire a p-value:\n", "1. Analytical sampling distribution: a $\\chi^2$ distribution of $m \\times k -1$ degree of freedom. We can refer to the $\\chi^2$ distribution table to acquire the p-value. If it is smaller than $0.05$, we will reject the null at the $95\\%$ confidence level.\n", "2. Empirical sampling distribution: a distribution constructed from a large number of $\\chi^2$ test statistics for simulations under the null of CSR. If the $\\chi^2$ test statistic for the observed point pattern is among the largest $5%$ test statistics, we would say that it is very unlikely that it is the outcome of a CSR process at the $95\\%$ confidence level. Then, the null is rejected. A pseudo p-value can be calculated based on which we can use the same rule as p-value to make the decision:\n", "$$p(\\chi^2) = \\frac{1+\\sum^{nsim}_{i=1}\\phi_i}{nsim+1}$$\n", "where \n", "$$ \n", "\\phi_i =\n", " \\begin{cases}\n", " 1 & \\quad \\text{if } \\psi_i^2 \\geq \\chi^2 \\\\\n", " 0 & \\quad \\text{otherwise } \\\\\n", " \\end{cases}\n", "$$\n", "\n", "$nsim$ is the number of simulations, $\\psi_i^2$ is the $\\chi^2$ test statistic calculated for each simulated point pattern, $\\chi^2$ is the $\\chi^2$ test statistic calculated for the observed point pattern, $\\phi_i$ is an indicator variable.\n", "\n", "We are going to introduce how to use the **quadrat_statistics.py** module to perform quadrat based method using either of the above two approaches to constructing the sampling distribution and acquire a p-value.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Juvenile Example" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import libpysal as ps\n", "import numpy as np\n", "from pointpats import PointPattern, as_window\n", "from pointpats import PoissonPointProcess as csr\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import the quadrat_statistics module to conduct quadrat-based method. \n", "\n", "Among the three major classes in the module, **RectangleM, HexagonM, QStatistic**, the first two are aimed at imposing a tessellation (rectangular or hexagonal shape) over the minimum bounding rectangle of the point pattern and calculate the number of points falling in each cell; **QStatistic** is the main class with which we can calculate a p-value, as well as a pseudo p-value to help us make the decision of rejecting the null or not." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import pointpats.quadrat_statistics as qs" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['HexagonM',\n", " 'QStatistic',\n", " 'RectangleM',\n", " '__all__',\n", " '__author__',\n", " '__builtins__',\n", " '__cached__',\n", " '__doc__',\n", " '__file__',\n", " '__loader__',\n", " '__name__',\n", " '__package__',\n", " '__spec__',\n", " 'math',\n", " 'np',\n", " 'plt',\n", " 'scipy']" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dir(qs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Open the point shapefile \"juvenile.shp\"." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "juv = ps.io.open(ps.examples.get_path(\"juvenile.shp\"))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "168" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(juv) # 168 point events in total" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[94., 93.],\n", " [80., 95.],\n", " [79., 90.],\n", " [78., 92.],\n", " [76., 92.],\n", " [66., 93.],\n", " [64., 90.],\n", " [27., 70.],\n", " [58., 88.],\n", " [57., 92.],\n", " [53., 92.],\n", " [50., 90.],\n", " [49., 90.],\n", " [32., 90.],\n", " [31., 87.],\n", " [22., 87.],\n", " [21., 87.],\n", " [21., 86.],\n", " [22., 81.],\n", " [23., 83.],\n", " [27., 85.],\n", " [27., 84.],\n", " [27., 83.],\n", " [27., 82.],\n", " [30., 84.],\n", " [31., 84.],\n", " [31., 84.],\n", " [32., 83.],\n", " [33., 81.],\n", " [32., 79.],\n", " [32., 76.],\n", " [33., 77.],\n", " [34., 86.],\n", " [34., 84.],\n", " [38., 82.],\n", " [39., 81.],\n", " [40., 80.],\n", " [41., 83.],\n", " [43., 75.],\n", " [44., 81.],\n", " [46., 81.],\n", " [47., 82.],\n", " [47., 81.],\n", " [48., 80.],\n", " [48., 81.],\n", " [50., 85.],\n", " [51., 84.],\n", " [52., 83.],\n", " [55., 85.],\n", " [57., 88.],\n", " [57., 81.],\n", " [60., 87.],\n", " [69., 80.],\n", " [71., 82.],\n", " [72., 81.],\n", " [74., 82.],\n", " [75., 81.],\n", " [77., 88.],\n", " [80., 88.],\n", " [82., 77.],\n", " [66., 62.],\n", " [64., 71.],\n", " [59., 63.],\n", " [55., 64.],\n", " [53., 68.],\n", " [52., 59.],\n", " [51., 61.],\n", " [50., 75.],\n", " [50., 74.],\n", " [45., 61.],\n", " [44., 60.],\n", " [43., 59.],\n", " [42., 61.],\n", " [39., 71.],\n", " [37., 67.],\n", " [35., 70.],\n", " [31., 68.],\n", " [30., 71.],\n", " [29., 61.],\n", " [26., 69.],\n", " [24., 68.],\n", " [ 7., 52.],\n", " [11., 53.],\n", " [34., 50.],\n", " [36., 47.],\n", " [37., 45.],\n", " [37., 56.],\n", " [38., 55.],\n", " [38., 50.],\n", " [39., 52.],\n", " [41., 52.],\n", " [47., 49.],\n", " [50., 57.],\n", " [52., 56.],\n", " [53., 55.],\n", " [56., 57.],\n", " [69., 52.],\n", " [69., 50.],\n", " [71., 51.],\n", " [71., 51.],\n", " [73., 48.],\n", " [74., 48.],\n", " [75., 46.],\n", " [75., 46.],\n", " [86., 51.],\n", " [87., 51.],\n", " [87., 52.],\n", " [90., 52.],\n", " [91., 51.],\n", " [87., 42.],\n", " [81., 39.],\n", " [80., 43.],\n", " [79., 37.],\n", " [78., 38.],\n", " [75., 44.],\n", " [73., 41.],\n", " [71., 44.],\n", " [68., 29.],\n", " [62., 33.],\n", " [61., 35.],\n", " [60., 34.],\n", " [58., 36.],\n", " [54., 30.],\n", " [52., 38.],\n", " [52., 36.],\n", " [47., 37.],\n", " [46., 36.],\n", " [45., 33.],\n", " [36., 32.],\n", " [22., 39.],\n", " [21., 38.],\n", " [22., 35.],\n", " [21., 36.],\n", " [22., 30.],\n", " [19., 29.],\n", " [17., 40.],\n", " [14., 41.],\n", " [13., 36.],\n", " [10., 34.],\n", " [ 7., 37.],\n", " [ 2., 39.],\n", " [21., 16.],\n", " [22., 14.],\n", " [29., 17.],\n", " [30., 25.],\n", " [32., 26.],\n", " [39., 28.],\n", " [40., 26.],\n", " [40., 26.],\n", " [42., 25.],\n", " [43., 24.],\n", " [43., 16.],\n", " [48., 16.],\n", " [51., 25.],\n", " [52., 26.],\n", " [57., 27.],\n", " [60., 22.],\n", " [63., 24.],\n", " [64., 23.],\n", " [64., 27.],\n", " [71., 25.],\n", " [50., 10.],\n", " [48., 12.],\n", " [45., 14.],\n", " [33., 8.],\n", " [31., 7.],\n", " [32., 6.],\n", " [31., 8.]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "juv_points = np.array([event for event in juv]) # get x,y coordinates for all the points\n", "juv_points" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Construct a point pattern from numpy array **juv_points**." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pp_juv = PointPattern(juv_points)\n", "pp_juv" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Point Pattern\n", "168 points\n", "Bounding rectangle [(2.0,6.0), (94.0,95.0)]\n", "Area of window: 8188.0\n", "Intensity estimate for window: 0.02051783097215437\n", " x y\n", "0 94.0 93.0\n", "1 80.0 95.0\n", "2 79.0 90.0\n", "3 78.0 92.0\n", "4 76.0 92.0\n" ] } ], "source": [ "pp_juv.summary()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHwVJREFUeJzt3X+sXOV95/H3995rb3Ag8oWQYGNs44U6bazSYBe8TTcbQpU14JQoIgkhi0wWyq6WpmlJVQhdERqFbCJVTaiwonpNE6+WX4lNyo8NCEqdLWxlJ74mKSRshOVyzcUOGHOdZGuX++u7f8wZMx7PjzMz5+dzPi8J4Zk7M+c5zzn3e7/n+zzPGXN3RESk/IbyboCIiCRDAV1EJBAK6CIigVBAFxEJhAK6iEggFNBFRAKhgC6lYGa3mNnmvNshUmQK6JIpM3vRzI6a2f8zs1fM7BtmdnK397n7l9z9upjbuM3M/ufgrW37+d8zs+uannMzOyetbYrEoYAuefiQu58MnA/8JvBfc25PrsxsJO82SBgU0CU37v4y8CiwCsDMFpvZQ2b2upntMbPfq7+2Mes2s+VRRrzBzPaZ2Wtm9qfRz9YBtwAfj64CftRq29GVwufM7CdmNhldKbwl+tmomT1iZgejnz1iZkuin90O/Fvgzujz7zSzv48+9kfRcx+PXrvezH5oZofN7B/M7Nebtn+Tmf0j8M9mNhI998dm9o9m9nMzu7/eJpE4FNAlN2Z2FnAp8Ez01L3ABLAYuAL4kpld3OEjfhtYCVwM3Gpmv+rujwFfAu5395Pd/bwO7/8k8O+Bfw38Cm9eKQwB3wCWAUuBo8CdAO7+p8BTwO9Hn//77v6+6H3nRc/db2bnA38N/CfgNOCvgIfM7F81bP8TwGXAQnefiZ77GLAOOBv4deCaDu0XOY4CuuThb8zsMPA08L+pBe6zqAXom9z9X9z9h8Bm4OoOn/Nn7n7U3X8E/AjoFLxbudPdX3L314HbqQVY3P2Qu29z9yPu/svoZ/+ux8/+PeCv3H2nu8+6+xbgDWBtw2v+Mtr+0abn9kdtehj4jR63KxWm2p3k4cPu/reNT5jZYuD1KIDWjQNrOnzOzxr+fQToOrja5KWmbS2O2rIA+Cq1THk0+vkpZjbs7rMxP3sZsMHMPt3w3Pz6Nlpsv655nxa3eI1IS8rQpSj2A6ea2SkNzy0FXu7js+LeQvSspm3tj/79WWqlnAvd/W1AvaRiPXz+S8Dt7r6w4b8F7n5vH+0UiUUBXQrB3V8C/gH4b2b2lmgA8Vrg7j4+7hVguZl1O79vMLMlZnYqtYHU+6PnT6FWNz8c/ezzLT5/RZfn/jvwn83sQqt5q5ld1vQHSyRRCuhSJJ8AllPLlL8DfN7dn+jjc74d/f+Qme3u8Lp7gMeBvdF/X4ye/xpwEvAasAN4rOl9dwBXRDNg/jJ67jZgSzSj5WPuvotaHf1OYBLYgwY4JWWmL7iQKjKzF4Hrmmv5ImWmDF1EJBAK6CIigVDJRUQkEMrQRUQCkenCore//e2+fPnyLDcpIlJ6Y2Njr7n76d1el2lAX758Obt27cpykyIipWdm43Fep5KLiEggFNBFRAKhgC4iEggFdBGRQCigi4gEQgFdRCQQCugiUgpj45Ns3L6HsfHJvJtSWPrGIhEpvLHxST65eQdTM3PMHxni7uvWsnrZaPc3VkxpAvpjjz3G9PR03s0QkRw8PgFvTBuOMTU9yze/+zT7l+TdqvjmzZvHunXrUt9OaQL69PQ0H/rQh/JuhojkYPH4JE9s3sH0zBzzRoa55tJyZegPP/xwJtspTUAXkepavWyUu69by469h1i74rRSBfMsKaCLVMzY+GQpA+PqZaOlam/d2Pgkj0/UrjLSbr8CukiFaHAxW/X+fmPaeGLzjtT7W9MWRRLQz5S6PKbh7dh7iKmZOeYcpmfm2LH3UGbbbtRu30Obmljvb8cy6W9l6CID6ifrzStTXrviNOaPDEWDi0OsXXFa6tts1m7fQ7x6qPf31PQs80aGU+9vZegiA+on680rU64PLt74wZW5Bcx2+16Uq4ck1fv7smWeSX8rQxcZUD9Zb56Zct6Di+32vQhXD2lYvWyU/UvIpM8V0KUnZZ0h0YvGfQS67m8/U+rSnoZX5OPUbt9bPV/k/SgiBXSJLcQaZ7PGfRwZHgJ3Zua86/72k/WmlSmX4Ti12/fG58uwH0WjGrrEFmKNs1nzPk7Peun2N5TjFMp+ZEkBXWKr1ziHjaBqnI2a93HesJVuf0M5TqHsR5ZUcpHYqrD8unkfoXsNvWhCOU6h7EeWFNClJ3nPkMhC8z6mtb+9Dr72okzHqdPAZ3NNXcG9MwV0kRz0O/gamrgDnxogjUc19IJrXAqd17Lo0JZjF0Evg68h93/cgU8NkMajDL3AipDFKTNKR+MimuHo2M7O+QmDf6H3f9zFRKEuOkqaAnqB7dh7iDem53BgamYOA5w3M5QsfrFbZUbttqsaZ3xxB1976f8yijvwqQHSeBTQC2x0wXy84XGUyGWaocTNjELPJNMQZ/C1Cplp3AHcMg305kUBvcAmj0wxZDDnMGTw8d9cyuKFJw2UofSaRcfNjELNJAe96hj0/cpMpRcK6AXWnJ195PwlA/1C95tFx8mMQswkB73qSOqqRZmpxKWAXmBJZ2dpZtEhZpKD9leoVy1SXAroBZdkdrZ2xWmMDBnTs87wkHXMohtLBT/92S959LkDXLJqEVdduDSRtpRB81XH6IL5bNy+J/YfrBCvWqTYFNCrxqK5MmZtX9JYKhgymJmrPf/UC68BtAzqIQ6KNl51jC6Yzxce+XFP+xfiVYsUmxYWlVivC0527D3E9ExtGuRMzEUc9WBed/8P9rXcZtoLP/JaXLN62Sg3XHQOk0emjtu/bbsnYrWn/n4Fc8mCMvSS6icjbpwGORc9bqWxVGB2fFD/yYFf8OzLPz9hm2mWF4qQ/R+3EGjI2Do2wcxsOFcjEgYF9JLqZ8Bt8sjUscVJFj1upblUUK+hnzRvmL99/pWW2+y3vBBnWl8ag4uDTN/cf/go935/X2UHO+P2nRaaZU8BvaT6yYgbM3SnfYYOxw/Grl42ylUXLmVsfJK/f+Fg2232OoAbN/NOOvsfdPrm2Pgk23ZPVHKwUzfTKjYF9BJpznh6zYibFyq1y9DbSfo7H5sz7227J1p+Vi/72qk99Z/tP3z0zVsqTLffbi/90I8yZrDtxkqa9yPu6yRZsQK6mf0RcB21xO5Z4FPAIuA+4FRgN3C1u/cWISS2dhlPL78YSWS6SX7nYy916Tj72qk9zTN3GscSvr3rJWZ7vOnZoNNJy5rBtprK2Wo/4r5OktV1louZnQn8AbDG3VcBw8CVwFeAr7r7ucAkcG2aDa26JGaR1DPLGz+4MpFfqEHb1Niej645i5nZwfavU3safzbbMMhrwEwO3xta1tvBNp9DzbN/6vsR93WSrLgllxHgJDObBhYAB4APAFdFP98C3AZ8PekGSk1SdeSkFyollfG3qkv3WpLo1J7mqwHMmJ2t/dupBfXh4eLd9CxP7fq/+Rxqtx9xXyfJ6RrQ3f1lM/tzYB9wFHgcGAMOu/tM9LIJ4MzUWimFXKSSZJta3U6210v0Tu1p9fn1BUO3PfRc7UXurT42FUU8no3iloTi7kfR9zcUXQO6mY0ClwNnA4eBbwOXtHhpy98GM7seuB5g6dLqLBtPQxFu0tSctSXZpsbP2rh9T19TFTu1p/lnq5eNsnH7HmbmHAdm5zzTKYhFOJ7t9DJVNO5+FHl/QxGn5PI7wD+5+0EAM3sA+C1goZmNRFn6EmB/qze7+yZgE8CaNWuyS4EkcVkO5GVVkihD6SMP6pdyihPQ9wFrzWwBtZLLxcAuYDtwBbWZLhuAB9NqpBRDlncPzOoSXaWA1tQv5RSnhr7TzLZSm5o4AzxDLeP+X8B9ZvbF6Lm70myo5C/rrC2rS3SVAlpTv5RPrFku7v554PNNT+8FLki8RVJYSS8sEkmKzsMarRSVniS5sEgkCToP36Tb50pPGm9jW9bFMRIWnYdvUoYusTVnQreuf7dmQkjuNCPnTQroEltzJjR5ZEozISR3mpHzJgV0ia1VJlSEmRBpDohpsK0cinAeFoECusRWxEwozQExDbZJ2SigS0+Klgmludhpx95Dx903vWrfTCTlo1kuUmr1MtCwkfiAWNzvYBUpCmXoUmpploEG/YYnkawpoEvppVUG0nQ4KRsFdJE2ijgILNKJArpIB0UbBBbpRIOiIiKBUECXlhrv2RLyNkOlvqwmlVzkBHksqNEinuSoL6tLGXrFxMnckrh7Xa8Zou6Ylxz1ZXUpQ6+QuJnboNP1+skQNUUwOerL6lJAD1zjzaXiLpMfdLpeP8vxNUUwOerL6lJAD9gg9y8fZLpevxmipggmR31ZTQrogWh1m9ek7l/e6y1kq5gh6ja7UgQK6AFoV7NO4v7l/c6YqFKGqFklUhSa5RKAdrMa6pnyjR9c2XeQSWPGRGhzpDWrRIpCGXoAOtWsB82Uk54xEWI2q1klUhQK6AFIs2ad9Gen+YUUeanimIEUkwJ6INKsWSf52aFms1UaM5DiUkCXTCmbFUmPArpkTtmsSDo0y0VEJBAK6CIpCm2KphSbSi4iKQlxiqYUmzJ0kZQMsuBImb30Qxm6SEr6naKpzF76pYAukpJ+p2iGuPhKsqGALpKifqZopr34SneGDJcCukjBpLn4SuWcsCmgC1CurC1uW5Pepyz7KK3FVyrnhC1WQDezhcBmYBXgwH8EfgrcDywHXgQ+5u4aki+hMmVtcdua9D6VqY86CfVeOlITd9riHcBj7v4u4DzgeeBm4El3Pxd4MnosJVSm+3l3a2t9ut+23ROJ7lOZ+qiTJO6RL8XVNUM3s7cB7wOuAXD3KWDKzC4H3h+9bAvwPeCmNBop6SpT1taprY1Z9MjwECNDxuycJ7JPZeqjbnQvnXDFKbmsAA4C3zCz84Ax4DPAO939AIC7HzCzd7R6s5ldD1wPsHTp0kQaLckq0x0QW7W1Xtt++fDRY1n07OwcV16wlMULT0pkn5q3C7Bx+57C95dUS5yAPgKcD3za3Xea2R30UF5x903AJoA1a9Z4X62U1JUpa2tsa6es/CPnL0n8yz7qf0BCqKdLeOIE9Algwt13Ro+3Ugvor5jZoig7XwS8mlYjRdpprG0nnZXH2aZmikiRdA3o7v4zM3vJzFa6+0+Bi4GfRP9tAL4c/f/BVFsq0kJzbTvprDzONstcT5ewxJ2H/mngbjObD+wFPkVthsy3zOxaYB/w0XSaKNJeHvX/Mo05SLXECuju/kNgTYsfXZxsc0R6l2b9v91ion63WaYFXFI+Wikq0oYWJ0nZ6H7okrmy3Ou702KifvYhlMVJUlzK0CVTZcpS2w1+9rsPGkyVtCmgS6aas9RtuycKW1NuN/jZLtPuth8aTJW0KaBLphqz1OEhY+vYBDOzxc3WWw1+Nmfaowvmx87Yy7SAS8pHAV0y1Zil7j98lHu/v690C3SaM20tNJKiUECXzDUuod+2e6KUNeXmTFu1cSkCBXTJTSg15VD2Q8pPAb0DLQJJXyg15VD2Q8pNAb2NMk2vExGBCi0s6nUhiBaBZK/TMSrLYqSQVKnP4+5r0fukEhl6P9m2FoFkq9Mx0tVS9qrU53l9T20aKpGh95Nt67sXs9XpGOlqKXud+vyenfu4+q6d3LNzX44tHFw9236g4ftnpzqcX60WxRUtW69Eht5vtq2Brux0Oka6Wspeuz6/Z+c+bvnOswA89cJrAFx1Yfm+WrIx2x4ymIu+S23OYXTB/JbvKcOiuEoEdE0ry0Y/s4Ia39PuGOn4ZW/1slFuXf9uHn3uAJesWnSszx997sBxr3v0uQOlDOiN2bY3fDHmEDB5ZKrle8qwKK4SAR2Ubaetn/piq/fccNE5LV+r45etsfFJvvDIj5mameMHL77OyjNOYfWyUS5ZtehYZg5wyapFObayf83ZNmbMzna/Aiz6orjKBHRJVz/L3/NcMq81Bp21Ozb1bLyeuZcxO4cTr/qg+83VOr2/KOeQArokop86d1618TLMVshbp2Nz1YVLSxvIGzVf9fV6DhTxqlEBXRLRT8aSV5ajm2l1V9QMVDpTQJfE9JOx5JHlaNZMPP0cm8ZSFvRWxpDBKaBL5Sj7TEdjKWskGmgs2rS+0CmgSyaKNghZxPpn2R1Xypp1wHFU1sqSArqkToOQ1dDvVEBJjgK6pE6DkNXQbSrgPTv3lX66Y9EpoEvqNAhZHe2mAoZyy4Ciq8TNuSRfzTc6Awp3UyNJV6tbBkjylKFLJhqXTKueXj2h3DKg6BTQJVOqp1dTKLcMKDoFdMmU6unl1+8U1JVnnMLkkSlWnnFKiq2rNgV0yZQW9ZRbvyUzldqyoUFRydzqZaPccNE5+oUuoX6/PSqJb50q+vd5FoEydBGJrd+S2aClNmX48Sigi0hs/ZbMBi21pT2YXrRbU/RLAV1EetLvfXAGuX9OmoPpIWX/CugJCOWvu8ig0vpdSHMwPaSptLEDupkNA7uAl919vZmdDdwHnArsBq5299bfrhqwkP66iwwi7d+FtO6QGdJU2l5muXwGeL7h8VeAr7r7ucAkcG2SDSuLXkbvNUovIUtiJksemm9NUeaELFaGbmZLgMuA24EbzcyADwBXRS/ZAtwGfD2FNhZa3L/uyuQldGXOdEO5P37cksvXgD8B6ku8TgMOu/tM9HgCOLPVG83seuB6gKVLw1vuG7e2F1KdTqQVLRrLX9eAbmbrgVfdfczM3l9/usVLvdX73X0TsAlgzZo1LV/Tzdj4JI9PwOLxyUKeJHH+uieRvWjwVYouq0xXvwutxcnQ3wv8rpldCrwFeBu1jH2hmY1EWfoSYH8aDayXKt6YNp7YvKO0pYpBsxeVbERq9LvQXtdBUXf/nLsvcfflwJXA37n7J4HtwBXRyzYAD6bRwHqpwrFSDbS00s+S9/pA6rbdE6UccAqVBrjTEadfyzr4moVB5qHfBNxnZl8EngHuSqZJx6uXKqamZ5k3MlyqgZZBHfct6sNDjAwZs3NeugGn0ChDTEfcfi3z4Gvaegro7v494HvRv/cCFyTfpOPVSxXf/O7TXHNp/784Zay5NWYis7NzXHnBUhYvPKlU+1BUg5wPGuBOR7vMu/k4afC1vVKsFF29bJT9SxgomJcxo2rORD5y/pJStLvoBj0flCGmo7lfRxfMb3ucQplmmLRSBPRBJZ1RZZXtr142yq3r333sW150Ag+mftxePnx0oPMhjeNSxivIpDVn3roS6l0lAnqSGVWW2f7Y+CRfeOTHTM3M8YMXX2flGafohO5TkuMRSR+Xsl5BpqE589aVUG8qEdCTrLllmTUoQ0lOkuMRSR8XHefWVCvvXSUCOiRXc8uyftqqprhx+x6d3H1Icjwi6XOgl8+rWmlGtfLeVCagJyXLrKFxW6ML5h+7zK/6ZXk/kjxuSZ8DcT9PpRnpRgG9SZwMKMusob6tjdv36LJ8QEket6TPgU6fl9RgroRPAb1BkTMgTZWrJi0uk14ooDco8uCUBoh6F0K9udtgbgj7KMlRQG9Q9CxYA0TxFflqqxedBnND2UdJjgJ6A2XB6ck6kyzy1VYvOp2ToeyjJEcBvYmy4OTlkUkW/WqrF+3OyZD2UZKhgC6pyyOTrMLVVhX2UXqjgC6pyyuTrMLVVhX2UeJTQJfUKZMUyYYCeomUeYpalTLJMh8nKTcF9JLQFLVy0HGSPHX9TlFJn75HMRxlPk76ntTyU4aeM32PYljKepx0ZREGBfScxZ3SV4WBxRBqz2U9TlqkFAYF9Jz1ktGFPLAYUoZYxuNU1isLOZ4Ces7KmtElrYoZYppXJL1+ts7DMCigF0AZM7qkVS1DTPOKpN/P1nlYfprlIj1JayZEPUO88YMrjwWge3bu4+q7dnLPzn2JbisJg/ZDmrNhyjzTRgajDF1iS7vO3Zgh3rNzH7d851kAnnrhNQCuunBpYtsaRBL9kOYVSdWuduRNCugSW5Z17kefO3DC46IE9CT6Ic2aterh1aWALi21GlTLMvO7ZNWiY5l5/XGntmWhvt3RBfMT6Yc0a9ZpfnYI00tDpYAuJ2hXUsgy86tn448+d4BLVi069jiv6Y3N2711/buZPDJVuaAW0vTSECmgywk6lRSynAlx1YVLTyiztBvwS/uPTPN2J49MccNF56SyrX5lkTlXcXppmSigywmKPKjW3LbRBfMzyRiL3CeQXeZc9H6oOgV0OUGRB9Wa25ZVxljkPoETM+dtuyc04FpBCujSUpEXmTS3LauMsch90pg5Dw8ZW8cmmJlNf3qpFIsCupRaFTLGOLXxxn7Yf/go935/n+rcFaSALqUXcsbYS2283g9j45Ns2z2hOncFKaCLFFg/YwRVuGqR1hTQRQqs31klIV+1SHtdA7qZnQX8D+AMYA7Y5O53mNmpwP3AcuBF4GPuru+uEkmQsm3pRZwMfQb4rLvvNrNTgDEzewK4BnjS3b9sZjcDNwM3pddUkWpSti1xdb19rrsfcPfd0b9/CTwPnAlcDmyJXrYF+HBajRQRke56uh+6mS0H3gPsBN7p7gegFvSBd7R5z/VmtsvMdh08eHCw1oqISFuxA7qZnQxsA/7Q3X8R933uvsnd17j7mtNPP72fNoqISAyxArqZzaMWzO929weip18xs0XRzxcBr6bTRBERiaNrQDczA+4Cnnf3v2j40UPAhujfG4AHk2+eiIjEFWeWy3uBq4FnzeyH0XO3AF8GvmVm1wL7gI+m00TJgr60QKT8ugZ0d38asDY/vjjZ5kge9KUFImHoaZaLhEnfEi8SBgX0ghsbn2Tj9j2Mjae3CLe+vHzYKOXNnAbtoyz6WCQLupdLgWVVCinz8vJB+0jlJgmJMvQCy7IUsnrZKDdcdE7pgtmgfaRyk4REAb3Ayl4KycKgfaQ+lpCo5FJgzaUQgI3b95SuLJKmTuWiXr/pR/0qZaeAXnCN30KjWm9rre5G2M83/YiUnUouJaFab2/UX1JFCugloVpvb9RfUkUquZSEar29UX9JFSmgl4hqvb2J21+6j42EQgFdKk2DzRIS1dCl0jR4KiFRQJdK0+CphEQlF6m0boOnqq9LmSigS+W1GzxVfV3KRiUXkTZUX5eyUUAXaUP1dSkblVykq0515JBrzFqcJGWjgC4ddaojV6HGrMVcUiYquUhHnerIqjGLFIsCunS0dsVpjAwPYcDw8PF1ZNWYRYpFJRfpzv34/0dUYxYpFgV06WjH3kPMzDkOzM4523ZPHBfAVWMWKQ4FdOmoXlaZnpljeMjYOjbBzGy4g6AiZVaagD5v3jwefvjhvJtRSf/lXfDCL2DyX+D/vGI4xtT0LN/87tPsX5J360SKb968eZlspzQBfd26dXk3ofLGxifZtXkH0zNzzBsZ5ppLlaGLFElpArrkT4OgIsWmgC490SCoSHFpHrqISCAU0EVEAqGALiISCAV0EZFAKKCLiARCAV1EJBDmTTdcSnVjZgeB8Q4veTvwWkbNKTr1RY36oUb9UFPVfljm7qd3e1GmAb0bM9vl7mvybkcRqC9q1A816oca9UNnKrmIiARCAV1EJBBFC+ib8m5AgagvatQPNeqHGvVDB4WqoYuISP+KlqGLiEifFNBFRAJRmIBuZuvM7KdmtsfMbs67PVkxs7PMbLuZPW9mPzazz0TPn2pmT5jZC9H/K3HPWjMbNrNnzOyR6PHZZrYz6of7zWx+3m1Mm5ktNLOtZvZ/o/Pi31TxfDCzP4p+J54zs3vN7C1VPB96UYiAbmbDwEbgEuDXgE+Y2a/l26rMzACfdfdfBdYCN0T7fjPwpLufCzwZPa6CzwDPNzz+CvDVqB8mgWtzaVW27gAec/d3AedR649KnQ9mdibwB8Aad18FDANXUs3zIbZCBHTgAmCPu+919yngPuDynNuUCXc/4O67o3//ktov75nU9n9L9LItwIfzaWF2zGwJcBmwOXpswAeArdFLgu8HM3sb8D7gLgB3n3L3w1TwfKD2BTwnmdkIsAA4QMXOh14VJaCfCbzU8Hgieq5SzGw58B5gJ/BOdz8AtaAPvCO/lmXma8CfAHPR49OAw+4+Ez2uwnmxAjgIfCMqPW02s7dSsfPB3V8G/hzYRy2Q/xwYo3rnQ0+KEtCtxXOVmk9pZicD24A/dPdf5N2erJnZeuBVdx9rfLrFS0M/L0aA84Gvu/t7gH8m8PJKK9EYweXA2cBi4K3USrLNQj8felKUgD4BnNXweAmwP6e2ZM7M5lEL5ne7+wPR06+Y2aLo54uAV/NqX0beC/yumb1IreT2AWoZ+8LokhuqcV5MABPuvjN6vJVagK/a+fA7wD+5+0F3nwYeAH6L6p0PPSlKQP8BcG40gj2f2uDHQzm3KRNRnfgu4Hl3/4uGHz0EbIj+vQF4MOu2ZcndP+fuS9x9ObXj/3fu/klgO3BF9LIq9MPPgJfMbGX01MXAT6jY+UCt1LLWzBZEvyP1fqjU+dCrwqwUNbNLqWVkw8Bfu/vtOTcpE2b228BTwLO8WTu+hVod/VvAUmon90fd/fVcGpkxM3s/8Mfuvt7MVlDL2E8FngH+g7u/kWf70mZmv0FtYHg+sBf4FLXkq1Lng5n9GfBxajPBngGuo1Yzr9T50IvCBHQRERlMUUouIiIyIAV0EZFAKKCLiARCAV1EJBAK6CIigVBAFxEJhAK6iEgg/j/vVQ/e7Wlv1wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pp_juv.plot(window= True, title= \"Point pattern\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rectangle quadrats & analytical sampling distribution\n", "\n", "We can impose rectangle tessellation over mbb of the point pattern by specifying **shape** as \"rectangle\". We can also specify the number of rectangles in each row and column. For the current analysis, we use the $3 \\times 3$ rectangle grids." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "q_r = qs.QStatistic(pp_juv,shape= \"rectangle\",nx = 3, ny = 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the plot method to plot the quadrats as well as the number of points falling in each quadrat." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuUVPWV6PHv7m4bQ4ThJdqKSHwEW3BC7BYwiQzqyBofd0bBQJT4ABxm7iJKhgmME3Nzx1mio0tHvYYh1/ExLDXRCJmlImJmAL3ETDPQYLQBHwwBaUAFaQaNCv3Y94+qaquLqupTVed99metXlDV9fid8zu1e5/9+/1OiapijDEm+qqCboAxxhh3WEA3xpiYsIBujDExYQHdGGNiwgK6McbEhAV0Y4yJCQvoxhgTExbQTWyIyAgRURGpCbotxgTBArrxlIjcKCJvisinIvK+iPyTiPxB0O3KR0Qmikirg8eNFZEVInJQRA6IyH+KyAwf2veKiNzk9fuY6LKAbjwjIn8N3A3MB/4AGA+MAH4lIsf43BYRkYqPdxE5H1gNvAqcAQwG/idwaaWvbUzFVNV+7Mf1H6A/8AkwNef+44APgRvSt/8FuCPr9xOB1qzbtwL/BXwMbAGuyvpdNXAvsB/YDswBFKhJ//4VYCHwGvAZqQA8A9iafr3twF+kH/vl9GO60u3+BDgpz3b9GljUy7b/ObANOAA8n3kdUn/MutuX1cab0v+/Mf369wJtwO+AS9O/Wwh0Ap+n2/aToPvYfsL3Yxm68co3gGOBX2bfqaqfAC8Bkxy+zn8BF5DK8G8HnhSRuvTv/hy4Avg60Ahcnef51wGzgX7ATlJ/TK4g9QdnBnC/iJyrqr8nlWXvUdXj0j97sl9IRPoC5wNLCzVWRC4C7gKmAnXp93za4bYCjAPeBoYA9wCPioio6m3AWuB76bZ9r4TXNAlhAd14ZQiwX1U78vxuL3C8kxdR1WdVdY+qdqnqM8C7wNj0r6cCD6jqLlU9QCqQ5voXVd2sqh2q2q6qL6rqf2nKq8CvSP3BcGIgqc/M3iKPmQ48pqobVfUw8LfA+SIywuF77FTVf1bVTmAJqT8KJzh8rkk4C+jGK/uBIQVmnNQB+5y8iIhcLyKvpwcgDwKjSf2xADgJ2JX18J15XiL794jIpSLSlB7MPAhclvV6vWkjVZKpK/KYk7LbkT4j+Qg42eF7vJ/13E/T/z3O4XNNwllAN175D+AwMDn7ThH5MqnSxqvpu34P9M16yIlZjz0V+Gfge8BgVR0AtACSfshe4JSs5w7P047u60OLSB9gGaka9Qnp11uR9XpFryWdDrD/AUwp8rA9wKlZ7/llUgOnu0ltKxTYXgfsWtemKAvoxhOq+t+kat4PicifiMgx6bLDs6Sy96fSD30duExEBonIicD3s17my6SC2D6A9NTA0Vm//wVwi4gME5GBpAZQi6kF+qRfr0NELqVnLf8DYHAv0yoXADeKyHwRGZxu19dEJFMn/xkwQ0TGpP+A3AmsU9UdqrqPVGD/rohUi8hM4PRe2pztA+C0Eh5vEsYCuvGMqt4D/JBURvwxqVkbfYE/Tg9CAjwB/BbYQaqe/UzW87cA95HKij8AziE1YyXjn4GX08/fSM4AbJ72fAzcQuoPQRtwLalZKJnfvwX8HNieLvGclOc1fgNclP7ZLiIHgIdJZfqo6irgf5E6E9hLKmB/J+sl/pzUNM6PgFHAb4q1OceDwNUi0iYi/6eE55mEEFU7izP+SGektwPfVNX3gm6PMXFjAd34SkSuA9pVtZSpfMYYByygG2NMTFgN3RhjYsLXq9INGTJER4wY4edbGmNM5DU3N+9X1V4X4/ka0EeMGMGGDRv8fEtjjIk8Ecm3aO4oVnIxxpiYsIBujDExYQHdGGNiwgK6McbEhAV0Y4yJCQvoxhgTExbQjTGR0LyzjUVrttG8sy3opoSWr/PQjTGmHM0725j+SBNHOrqoraniqZvG03DqwKCbFTrRCegivT/GGBNLDcBb2XfcEVBDKuHDdbOiE9DBlx1iSiBifRJGMeyXTIbe3tHFMVHM0H1KSKMV0I0xidRw6kCeumk8Tds/Yvxpg6MVzH1kAd2YhGne2RbJwNhw6sBItTejeWcbDZl/PW6/zXIxJkEypYv7fvU20x9pshkjHsvsb8CX/W0B3RgXlDOlLohpeE3bP+JIRxddCu0dXTRt/8i3985WaNvjNjUxs7/Bn/1tJRdjKlTOlLqgpuGNP20wtTVV3YOL408b7Pl75iq07XGcmpjZ34Av+9sydGMqVE7WG1SmnBlcnDdpZGABs9C2h+XswU2Z/Q34sr8tQzemQuVkvUFmykEPLhba9jCcPXghs6/92Oe+fkl0Y2Ojlv2NRTGcWxtFPWZIjBgUyz7J3kbA0YyQcmaOeDbbRITmHQdCPZOl0Lbn3h/VGTlHqTB+iUizqjb2+jgL6Map3BrnW3dcFrs+yd7GmuoqUKWjS6NV0xXhrB+tiHwtOlY1dZ8CutXQjWO5Nc44yt3G9k6NZE03DrXoONbUvWYB3TiWqXFWS2rEPo5yt/GYaun+f5RqutnbEKV2Z8vti6huh5+s5GJKYjX0CIhADd0pq6Fnnm41dOM165OKePaHI2L94jRoRzq4+xTQbdqiMQGIxeCrC5wOfMZqgNRD8SyExkj2UuiglkXHbTl2GJQy+Brn/e904NMGSJ2xDD3EwpDFWWbkjexFNNXpvu3s0qMG/+K+/50uJorroiO3WUAPsabtH3G4vQslNQ1NAOWLDMWPD3a+zKjQ+0a6xumz3Ot7Q/4aein7P4qcXufcrofujAX0EBvYt5bsYZR0IudrhuI0M4p7JumF3CX4+fZXEjJTp5ciCPqSBVFgAT3E2j49QpVAl0KVwLTzhnPSgC9VlKGUmkU7zYzimklWetZR6fMtMzWlsIAeYrnZ2eRzh1X0gS43i3aSGcUxk6z0rMOtsxbLTI1TFtBDzO3szMssOo6ZZKX7K65nLSa8LKCHnJvZ2fjTBlNTJbR3KtVVUjSLzi4VvP3+x7zUspdLR9dx7bjhrrQlCnLPOgb2rWXRmm2O/2DF8azFhJsF9KSR9FwZkYIPyS4VVAlkrsO19t39AHmDehwHRbPPOgb2reXvl28uafvieNZiws0CeoSVOuDWtP0j2jtS0yA7ipQAsksFXTmrlZ9Z/x5tnx5JvWeB53hRXghqSmTmDGnRmm09tm/ZxlZH7bH6t/GTBfSA7Nq1i+uvv57333+fqqoqZs+ezdy5c5k/fz4vvPACtbW1nH766Tz++OMMGDDgqOeXkxFnT4PsSt/OJ7tUIFkZOsCWvYd4c/d/p66HXuA5bpcX/Mr+P//8cyZMmMDhw4fp6Ojg6quv5vbbb2fWrFms/c06dh/4PccMPIkT/8c8lja30tEZn7ORMJs5cybLly9n6NChtLS0dN//0EMP8ZOf/ISamhouv/xy7rnnngBbGQ629D8gNTU13HfffWzdupWmpiYWLVrEli1buOSSS2hpaeGNN97gq1/9KnfddVfe55ezFLrt0yNkCi2Svp1P9vdOPvMX3+DOq87hgjOHMOnsE+js0rzXQy/3uyqdLGv3Ytl3vvft06cPq1ev5re//S2vv/46K1eupKmpifvvv593trbwm//cyNhzRnLK+2vp6EzuMnSnlyJw65IFN954IytXruxx35o1a3juued444032Lx5Mz/4wQ8qeo+4sAw9IHV1ddTV1QHQr18/6uvr2b17N5MmTep+zPjx41m6dGne55eTEWdn6ErhDB16lgoaTh3IteOG07yzjf/37r7u9yz2HCecZt5uZ/+F3ldEOO644wBob2+nvb0dEaF///4AnDt8ACOH9KHmDwawN6GDnUFcTGvChAns2LGjx32LFy/m1ltvpU+fPgAMHTq0rNeOGwvoIbBjxw42bdrEuHHjetz/2GOPMW3atO7buXXkUgfcchcqFcrQCzlqufodldW2czPvQnXpUra1WHsyv9tz8LMvLqnQ3vN9xwzrT0NDA9u2bWPOnDndfTJjxgxWrFjB2WefzYsv3seN+w5XXNOP4qUSCp0t5W6H08eV65133mHt2rXcdtttHHvssdx7772cd955lW1cDDgK6CLyV8BNpBK7N4EZQB3wNDAI2Ahcp6qlRQjDJ598wpQpU3jggQe6M0GAhQsXUlNTw/Tp04HCGU8pHww3Mt3c96wkC+txgaoqKVqXdrKtxbLC3Jk72WMJz27YRWfWRc9ef/11Dh48yFVXXUVLSwujR4/m8ccfp7Ozk5tvvplnnnmGGTNmBLLIK2j5pnLm2w6njytXR0cHbW1tNDU1sX79eqZOncr27duRIrO3kqDXGrqInAzcAjSq6migGvgOcDdwv6qeCbQBs7xsaBy1t7czZcoUpk+fzuTJk7vvX7JkCcuXL+epp57qPkDdqCOXW+cuppI2Zbfn242nVFyXLraPsn/XmVX+F6Ajz6VrBwwYwMSJE3vUbqurq5k2bRrLli0ruW2ltDXMco+htk+P5N0Op48r17Bhw5g8eTIiwtixY6mqqmL//v1ubGKkOS251ABfEpF2oC+wF7gIuDb9+yXA3wGL3W5gXKkqs2bNor6+nnnz5nXfv3LlSu6++25effVV+vbt232/W3Vkt6fRuZXxN+9sY9nG1h6vVWpJotg+yj0bQITOztT/lVRQ188/ZtTg1Efis88+49///d9ZsGAB27Zt44wzzkBVeeGFFzjrrLNK3s5S2hoWhfZ/7jFUaDucPq4cV155JatXr2bixIm88847HDlyhCFDhlT0mnHg6CvoRGQusBD4DPgVMBdoUtUz0r8/BXgpncEXZF9B94Vf//rXXHDBBZxzzjlUVaVOlO68805uueUWDh8+zODBqQN+/Pjx/PSnPwVCWHN1+bsrc7+SrZxTdCc19OzL1Q7sW8vfPd+S+oKJj3Zw7G/+L32qoauri6lTp/KjH/2ICy64gEOHDqGqfO1rX2Px4sU9ymNubK+r/enCZ6WUkpDXXyF3zTXX8Morr7B//35OOOEEbr/9dq677jpmzpzJ66+/Tm1tLffeey8XXXRRydvpm7B8p6iIDASWAdOAg8Cz6dv/Oyegr1DVc/I8fzYwG2D48OENO3fuLHFTul8oVgE9qvz6kuhFa7Zx36/epkuhWmDepJHMufCMyL6Pr1z4rMRyvwQpRN8p+sfA71R1X/qFfwl8AxggIjWq2gEMA/bke7KqPgw8DKkM3WH7TQjlZm1v9f6UsvlVkohC6SMItl+iyUlAfw8YLyJ9SZVcLgY2AGuAq0nNdLkBeM6rRppwyB3I85Jf10Gx663kZ/slmnoN6Kq6TkSWkpqa2AFsIpVxvwg8LSJ3pO971MuGmuDlZm1e8+s6KHa9lfxsv0SPo0FRt9igaPTl1tDdHBQ1LkngZyV0EwZyhaiGbkw3NxcWGeOGqC7S8oJdnMuUJPeCS1FcHGPiJaqLtLxgGbpxLN8sF5sJYYJmM3K+YAHdOJZvlovNhDBBsxk5X7CAbhzLN8slDDMhvBwQC/1gmwHCcRyGgQV041i+y+cGzcsBMRtsM1FjAd2UJGyZkJffZdq0/aMe1013+3tSjXGbzXIxkZYpA1ULrg+IOf0OVmPCwjJ0E2leDohV+g1PxvjNArqJPK/KQDYdzkSNBXRjCrDpcCZqLKAbU0TYBoGNKcYGRY0xJiYsoJu8cq/ZEtf3jCvbl8lkJRdzlCAW1NgiHnfZvkwmy9ATxknm5sbV60rNEO2Kee6yfZlMlqEniNMsuNLpeuVk2zZF0F22L5PJAnrMZV9cyuky+Uqn65WzHN+mCLrL9mUyWUCPsdxM+cdXjHKcuVUyXa/cbNumCLrH9mUyWUCPiXyXec3NlNs+PVJW5lbqJWSTmG3bZXZNGFhAj4FCNet8mXKpmVu5s0+SlCHaDB0TFjbLJQYKzRDJZMrzJo0sO8h4MfskbnOkbYaOCQvL0GOgWM260kzZ7dknccxmbYaOCQsL6DHgZc3a7df28gspgpLEMQMTThbQY8LLmrWbrx3XbDZJYwYmvCygG19ZNmuMdyygG99ZNmuMN6IV0EWCboHJZX0STtYviRStgK7a+2OMf0SsT3oRyIIj65fw8ekPbLQCujEREscpmibcbGGRMR6pZMFR3BZfGX9Yhm6MR8qdommZvSmXBXRjPFLuFM04Lr4y/rCAboyHypmi6fXiK7syZHxZQDcmZLxcfGXlnHizgG6AaGVtTtvq9jb5uY+8Wnxl5Zx4cxTQRWQA8AgwGlBgJvA28AwwAtgBTFVVG5KPoChlbU7b6vY2RWkfFRPXa+mYFKfTFh8EVqrqWcDXgK3ArcAqVT0TWJW+bSIoStfz7q2tmel+yza2urpNUdpHxbhxjXwTXr1m6CLSH5gA3AigqkeAIyLyZ8DE9MOWAK8Af+NFI423opS1FWtrdhZdU11FTZXQ2aWubFOU9lFv7Fo68SXayxJhERkDPAxsIZWdNwNzgd2qOiDrcW2qetRRIiKzgdkAw4cPb9i5c2eZLbXlzF4qqz4cUJ/ktjVze/fBz3j6P9+jS6Fa4Dtjh3PSgC95UkMHwjvmYJ+V8KmwT0SkWVUbe32cg4DeCDQB31TVdSLyIHAIuNlJQM/W2NioGzZscLQBeRpiB2nYhKBPcrNyVLuzcq9KCqGvp4egX0wOnwK6k0HRVqBVVdelby8lVS//QETqVHWviNQBH5bdWmPKlF3b7uzscj0r7+09baaICZNeA7qqvi8iu0RkpKq+DVxMqvyyBbgB+If0v8952lJj8sitbU8+d5jnwTVO9XQTL72WXKC7jv4IUAtsB2aQmiHzC2A48B7wbVU9UOx1rOQSMyHpkyDm0Id63n5I+sVkCUsN3U0W0GMmAX0SycVJCeiXyAlRDd2YRLLFSSZq7HroxndRudZ3scVE5WxDXBYnmfCyDN34KkpZaqHBz3K3wQZTjdcsoBtf5Wapyza2hnZwsdBVDwtl2r1th5dXUTQGLKAbn2VnqdVVwtLmVjo6w5ut51smn5tpD+xb6zhjt2X3xksW0I2vsrPUPQc/4+fppfpRWqCTm2nbQiMTFhbQje8yWWrzzjaWbWyNZE05N9O22rgJA5uHbsrnQp+EeoFOCUK1HfZZCR+f5qHbtMUivJxeN3PmTIYOHcro0aO773v22WcZNWoUVVVVlP2HL2IaTh3InAvPCD4IArt27eLCCy+kvr6eUaNG8eCDD3b/7qGHHmLkyJGMGjWKBQsWHPXcMG1HnBTqk/nz53PWWWfxh3/4h1x11VUcPHgw4JaGhKr69tPQ0KBlg/KfW4YNOw7oyB+t0K/culxH/miFbthxwNXXf/XVV7W5uVlHjRrVfd+WLVv0rbfe0j/6oz/S9evXu/p+nvC5T7y2Z88ebW5uVlXVQ4cO6ZlnnqmbN2/W1atX68UXX6yff/65qqp+8MEHQTazdzHql0J98vLLL2t7e7uqqi5YsEAXLFgQZDN7V2GfABvUQYxNTIZearbt9SKQCRMmMGjQoB731dfXM3LkSFffJ0qK9ZEfi5Hq6uo499xzAejXrx/19fXs3r2bxYsXc+utt9KnTx8Ahg4d6lkbwiQMC8AK9cmkSZOoqUkNAY4fP57W1taK3sfptoZhnxSTiEHRchaC2CIQfxXroyAWI+3YsYNNmzYxbtw45s+fz9q1a7nttts49thjuffeeznvvPM8ff+ghXEBWHafZHvssceYNm1a2a8b1PfUeiERAb2caWW2CMRfxfrI72mBn3zyCVOmTOGBBx6gf//+dHR00NbWRlNTE+vXr2fq1Kls374dEfGsDUErts9/tu49XmrZy6Wj67h23HBf2pPbJxkLFy6kpqaG6dOnl/yamYHsPQc/697WI0WOrygsiktEQC8327ZFIP4p1kd+ni21t7czZcoUpk+fzuTJkwEYNmwYkydPRkQYO3YsVVVV7N+/n+OPP96zdgSt0D7/2br3+OG/vgnA2nf3A3ge1PP1CcCSJUtYvnw5q1atKvmPa3a2XSXQlZ6A0qUwsG9t3udEYVFcIgK6Zdv+KGfqXvZzCvWRX/2nqsyaNYv6+nrmzZvXff+VV17J6tWrmThxIu+88w5HjhxhyJAhnrQhLBpOHciPrxjVnYln9vlLLXt7PO6llr2eBvRCfbJy5UruvvtuXn31Vfr27Vvy62Zn29mzCauAtk+P5H1OFBbFJSKgQ/iy7WuuuYZXXnmF/fv3M2zYMG6//XYGDRrEzTffzL59+7j88ssZM2YML7/8ctBNdaSc+mK+58y58Iy8j/Wj/1577TWeeOIJzjnnHMaMGQPAnXfeycyZM5k5cyajR4+mtraWJUuWxLrcAqm++fvlmznS0cX6HQcYeWI/Gk4dyKWj67ozc4BLR9d52o5CfXLLLbdw+PBhLrnkEiA1MPrTn/7U8evmZtuI0NnZ+xlg2BfFJSagh83Pf/7zvPdfddVVPrfEHeXUuYNcMp/vbOJb3/oWWmDxx5NPPulLu8KiUN9ksnG/auiF+uSyyy6r6HVzz/qg94urFXt+WJJFC+jGFeXUuYOaSRSF2QpBK9Y3144b7ttgqJdyz/pKPQbCdtYPFtCNS8rJWILKcuxiWr0LawZqirOAblxTTsYSRJZjawycKadvsktZUFoZw1TOArpJHMs+vZFdyqpJDzSGbVpf3FlAN74I1dUICWf9M+p6lLI6FVAUK2v5yQK68ZwNQiZDuVMBjXssoBvP2SBkMvQ2FTCISwYkjQV04zkbhEyOQlMBg7hkQBIl5vK5JjiZzG3epJE8ddN4gFBfgtS4L98lA4z7LEM3vsheMm319OTx+5IBSWUB3fjK6unJ5PclA5LKArrxldXTo6/cKagjT+xH26dHGHliPw9bl2wW0I2vbFFPtJVbMrNSmz9sUNT4ruHUgcy58Az7QEdQud+168Z39Ib9+zzDwDJ0Y4xj5ZbMKi21WYbvjAV0Y4xj5ZbMKi21eT2YHrZLU5TLAroxpiTlXgenkuvneDmYHqfs3wK6C+Ly192YSnn1WfByMD1OU2kdB3QRqQY2ALtV9QoR+QrwNDAI2Ahcp6r5v101xuL0192YSnj9WfDqCplxmkpbyiyXucDWrNt3A/er6plAGzDLzYZFRSmj9zZKb+LMjZksQci9NEWUEzJHGbqIDAMuBxYC8yT1lecXAdemH7IE+DtgsQdtDDWnf90tkzdxF+VMNy7Xx3dacnkAWABklngNBg6qakf6ditwcr4nishsYDbA8OHxW+7rtLYXpzqdMfnYorHg9RrQReQK4ENVbRaRiZm78zxU8z1fVR8GHgZobGzM+5jeNO9soyHzbwgPEid/3d3IXmzw1YSdX5mufRbyc5KhfxP4UxG5DDgW6E8qYx8gIjXpLH0YsMeLBmZKFW8B0x9pimypotLsxUo2xqTYZ6GwXgdFVfVvVXWYqo4AvgOsVtXpwBrg6vTDbgCe86KBmVIFRGugJZ9ylrxnBlKXbWyN5IBTXNkAtzec7NeoDr76oZJ56H8DPC0idwCbgEfdaVJPmVIFELmBlkr1+Bb16ipqqoTOLk3cfggbyxC94XS/Rnnw1WslBXRVfQV4Jf3/7cBY95vUU6ZUwR1U9MGJYs0tOxPp7OziO2OHc9KAL0VqG8KqkuPBBri9USjzzu0nG3wtLBIrRbM7shxRzahyM5HJ5w6LRLvDrtLjwTJEb+Tu14F9awv2U1ymGbotEgG9Um5nVH5l+w2nDuTHV4zq/pYXO4Ark+m33Qc/q+h48KJfongG6bbczNvOhEqXiIDuZkblZ7bfvLONv1++mSMdXazfcYCRJ/azA7pMbo5HuN0vUT2D9EJu5m1nQqVJREB3s+bmZ9ZgGYp73ByPcLtfrJ/zs1p56US1rLU+ZWlsbNQNGzaU92TJt5bJGGMiooJYKyLNqtrY2+OilaH7+MenGD/rnZn3Gti3tvs0PzSn5SKh6RMn3Ow3t48BJ6/nuDQTsX5JBJ8S0mgFdB84+WD5OcKeea9Fa7bZaXmF3Ow3t4+BYq/n1mCuiT8L6FnCPDhlU+WSyRaXmVJYQM8S5sEpGyAqXRymAvY2mBuHbTTusYCeJexZsC2mcC7MZ1ulKLa4LC7baNxjAT2LZcHe8TuTDPPZVimKHZNx2UbjHgvoOSwLdl8QmWTYz7ZKUeiYjNM2GndYQDeeCyKTTMLZVhK20ZTGArrxXFCZZBLOtpKwjcY5C+jGc5ZJGuMPC+gREuUpaknKJKPcTybaLKBHhE1RiwbrJxOkXr9T1HjPvkcxPqLcT/Y9qdFnGXrA7HsU4yWq/WRnFvFgAT1gTqf0JWFgMQ6156j2ky1SigcL6AErJaOL88BinDLEKPZTVM8sTE8W0AMW1YzObUnMEL08Iyn1te04jAcL6CEQxYzObUnLEL0+Iynnte04jD4L6KYkPTI/F183X4b4s3Xv8VLLXi4dXce144a7+G6VqzS79vqMJGlnOybFArpxLDerfMvl18/OEH+27j1++K9vArD23f0AoQnqbmTXXp+RJOlsx3zBArpxLDer9NJLLXuPuh2WgO5Gdu11zdrq4clkAd3kla+kkJtVeunS0XXdmXnmdrG2+SH7C7vdyIC9rFl7+dpxmF4aVxbQzVEKlRRys0ru8K4NmWw8t4Ye1PTG3Pf98RWjaPv0SOKCWpyml8aRBXRzlGIlBT9nQlw7bvhRZZZCS+u9zhhz37ft0yPMufAMT96rXJnMeY6H75HE6aVRYgHdHCXMUwhz2zawb60vGWOY9wn0zJznpG8ncT8knQX0gMycOZPly5czdOhQWlpaAJg/fz4vvPACtbW1nH766Tz++OMMGDDA97aFeZFJbtvczBh37drF9ddfz/vvv09VVRWzZ89m7ty53f3SKdUcN/hkHvynh0O1T6Bn5gywbGOrJ/0X1LHR2dlJY2MjJ598MsuXL/flPaPIrrYYkBtvvJGVK1f2uO+SSy6hpaWFN954g69+9avcddddAbUu9cGdc+EZoQtc0LNtmYyxWqg4Y6ypqeG+++5j69atNDU1sWjRIrZs2dLdL+9u3cyffPPrvPTUYhe3xh3Z+wFgaXMr9/3qbaY/0uT61RPxU8tcAAAIm0lEQVSDODYefPBB6uvrfXu/qLKAHpAJEyYwaNCgHvdNmjSJmprUSdP48eNpbW0NommRkskY500aWXG5pa6ujnPPPReAfv36UV9fz+7duwPvFyeXtc3eDwAdndG8hG8+ra2tvPjii9x0001BNyX0rOQSUo899hjTpk0LuhmR4MVA7Y4dO9i0aRPjxo3rcb/f/VLKrJLs/RCnOvf3v/997rnnHj7++OOgmxJ6lqGH0MKFC6mpqWH69OlBNyWRPvnkE6ZMmcIDDzxA//79u+8Pol/K/cIMt85agpYZZ2pocPNCE/FlGXrILFmyhOXLl7Nq1SpEJOjmJE57eztTpkxh+vTpTJ48ufv+oPql3FklcbnQ1muvvcbzzz/PihUr+Pzzzzl06BDf/e53efLJJ4NuWjipatEf4BRgDbAV2AzMTd8/CPg34N30vwN7e62GhgYtG5T/3JD63e9+p6NGjeq+/dJLL2l9fb1++OGHAbaqBDHrk66uLr3uuut07ty5Pe4Pul827DigP1n9rm7YccDZE2LWLxlr1qzRyy+/POhmlKfCPgE2aC/xVVUdlVw6gL9W1XpgPDBHRM4GbgVWqeqZwKr0bePQNddcw/nnn8/bb7/NsGHDePTRR/ne977Hxx9/zCWXXMKYMWP4y7/8y6CbmSivvfYaTzzxBKtXr2bMmDGMGTOGFStWBN4vYZ5xZMJFUsG/hCeIPAf8JP0zUVX3ikgd8Iqqjiz23MbGRt2wYUOZLRUosa3GY9Yn4WT9Ej4V9omINKtqY2+PK2lQVERGAF8H1gEnqOpegPS/Qws8Z7aIbBCRDfv27Svl7YwxxpTAcUAXkeOAZcD3VfWQ0+ep6sOq2qiqjccff3w5bTTGGOOAo4AuIseQCuZPqeov03d/kC61kP73Q2+aaIwxxoleA7qk5mg9CmxV1X/M+tXzwA3p/98APOd+84wxxjjlZB76N4HrgDdF5PX0fT8E/gH4hYjMAt4Dvu1NE40f7EsLjIm+XgO6qv4aKLSS4mJ3m2OCYF9aYEw82NJ/U/bycmNMuFhADzknV9qrlJuXoA1CpfvIj31sjB/sWi4h5lcpJMxfaNGbSveRlZtMnFiGHmJ+lkKiury80n1k5SYTJxbQQyzqpRA/VLqPbB+bOCn5Wi6VsGu5lC57OiF4/+32JQlJnxSacul0KmbspmyGpF9MFp+u5WIBPSJCWesNcZ+Ecn/5JcT9klhhvDiXCY7Vektj+8skkQX0iLBab2lsf5kkspJLhISu1hvyPgnd/vJLyPslkXwqudg89AiJy/dE+sXp/kps4DexYwHdJFqiB09N7FgN3SSaDZ6aOLGAbhLNBk9NnFjJxSRab9exsfq6iRIL6CbxCg2eWn3dRI2VXIwpwOrrJmosoBtTgNXXTdRYycX0qlgdOc415ihfJ94kkwV0U1RvdeS415htMZeJEiu5mKJ6qyNbjdmY8LCAbooaf9pgaqqrEKC6+ug6stWYjQkPK7mY3mUuKpTn4kJWYzYmPCygm6Katn9ER5eiQGeXsmxj6xcBHKsxGxMmFtBNUZmpe+0dXVRXCUubW+noTA2CvhV044wxPUQroIsE3YLEaYDigdv6JJysXxIpOgHdLtgfuMwUxvaOLo6J6TRFY6IsOgHdBM4W2hgTbhbQTUlsENSY8LJ56MYYExMW0I0xJiYsoBtjTExYQDfGmJiwgG6MMTFhAd0YY2JC1McFOyKyD9hZ5CFDgP0+NSfsbF+k2H5Isf2QktT9cKqqHt/bg3wN6L0RkQ2q2hh0O8LA9kWK7YcU2w8pth+Ks5KLMcbEhAV0Y4yJibAF9IeDbkCI2L5Isf2QYvshxfZDEaGqoRtjjClf2DJ0Y4wxZbKAbowxMRGagC4ifyIib4vINhG5Nej2+EVEThGRNSKyVUQ2i8jc9P2DROTfROTd9L+JuGatiFSLyCYRWZ6+/RURWZfeD8+ISG3QbfSaiAwQkaUi8lb6uDg/iceDiPxV+jPRIiI/F5Fjk3g8lCIUAV1EqoFFwKXA2cA1InJ2sK3yTQfw16paD4wH5qS3/VZglaqeCaxK306CucDWrNt3A/en90MbMCuQVvnrQWClqp4FfI3U/kjU8SAiJwO3AI2qOhqoBr5DMo8Hx0IR0IGxwDZV3a6qR4CngT8LuE2+UNW9qrox/f+PSX14Tya1/UvSD1sCXBlMC/0jIsOAy4FH0rcFuAhYmn5I7PeDiPQHJgCPAqjqEVU9SAKPB1JfwPMlEakB+gJ7SdjxUKqwBPSTgV1Zt1vT9yWKiIwAvg6sA05Q1b2QCvrA0OBa5psHgAVAV/r2YOCgqnakbyfhuDgN2Ac8ni49PSIiXyZhx4Oq7gbuBd4jFcj/G2gmecdDScIS0PN9RXmi5lOKyHHAMuD7qnoo6Pb4TUSuAD5U1ebsu/M8NO7HRQ1wLrBYVb8O/J6Yl1fySY8R/BnwFeAk4MukSrK54n48lCQsAb0VOCXr9jBgT0Bt8Z2IHEMqmD+lqr9M3/2BiNSlf18HfBhU+3zyTeBPRWQHqZLbRaQy9gHpU25IxnHRCrSq6rr07aWkAnzSjoc/Bn6nqvtUtR34JfANknc8lCQsAX09cGZ6BLuW1ODH8wG3yRfpOvGjwFZV/cesXz0P3JD+/w3Ac363zU+q+reqOkxVR5Dq/9WqOh1YA1ydflgS9sP7wC4RGZm+62JgCwk7HkiVWsaLSN/0ZySzHxJ1PJQqNCtFReQyUhlZNfCYqi4MuEm+EJFvAWuBN/midvxDUnX0XwDDSR3c31bVA4E00mciMhH4gapeISKnkcrYBwGbgO+q6uEg2+c1ERlDamC4FtgOzCCVfCXqeBCR24FppGaCbQJuIlUzT9TxUIrQBHRjjDGVCUvJxRhjTIUsoBtjTExYQDfGmJiwgG6MMTFhAd0YY2LCAroxxsSEBXRjjImJ/w9At+2KU/dihgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "q_r.plot()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "33.107142857142854" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_r.chi2 #chi-squared test statistic for the observed point pattern" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_r.df #degree of freedom" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.890978545159614e-05" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_r.chi2_pvalue # analytical pvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the p-value based on the analytical $\\chi^2$ distribution (degree of freedom = 8) is 0.0000589, much smaller than 0.05. We might determine that the underlying process is not CSR. We can also turn to empirical sampling distribution to ascertain our decision." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Rectangle quadrats & empirical sampling distribution\n", "\n", "To construct a empirical sampling distribution, we need to simulate CSR within the window of the observed point pattern a lot of times. Here, we generate 999 point patterns under the null of CSR." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "csr_process = csr(pp_juv.window, pp_juv.n, 999, asPP=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We specify parameter **realizations** as the point process instance which contains 999 CSR realizations." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "q_r_e = qs.QStatistic(pp_juv,shape= \"rectangle\",nx = 3, ny = 3, realizations = csr_process)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.001" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_r_e.chi2_r_pvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The pseudo p-value is 0.002, which is smaller than 0.05. Thus, we reject the null at the $95\\%$ confidence level." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hexagon quadrats & analytical sampling distribution\n", "\n", "We can also impose hexagon tessellation over mbb of the point pattern by specifying **shape** as \"hexagon\". We can also specify the length of the hexagon edge. For the current analysis, we specify it as 15." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "q_h = qs.QStatistic(pp_juv,shape= \"hexagon\",lh = 15)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VFUTh38nCQFpUoOBQCIQQgw9UUKRYghKCQgI0kGa0hRBARXFBqhY4BNQmorSpEkXqQGkSUIRECGUJDSpoZeUne+P2Q2bsMne3b1tN/d9nn0gu3fPmXvv3jlz5syZEUQEAwMDAwPPxUtrAQwMDAwMlMVQ9AYGBgYejqHoDQwMDDwcQ9EbGBgYeDiGojcwMDDwcAxFb2BgYODhGIrewMDAwMMxFL2BxyOECBJCkBDCR2tZDAy0wFD0BpoghOgthDgkhLgrhPhPCDFNCPG41nLZQgjRRAhxVsJxzwgh1gohrgshrgkh/hJCvKKCfLFCiH5K92PgvhiK3kB1hBAjAHwO4G0AjwOIBBAEYL0QIp/KsgghhMvPgRCiHoDNALYCqAygJICBAFq42raBgcsQkfEyXqq9ABQFcBtAp2zvFwZwCUAv898/AfjU6vMmAM5a/T0awEkAtwD8A6Cd1WfeAL4EcAXAKQCDARAAH/PnsQDGAdgB4B5YMb8C4Ki5vVMAXjUfW8h8jMks920AZW2c158Apto59/4ATgC4BmClpR3wIJcpn5WM/cz/721u/0sAKQBOA2hh/mwcgAwA982yTdH6Hhsv/b0Mi95AbeoDKABgmfWbRHQbwO8Amkts5ySAZ8Ezgo8AzBVC+Js/6w+gNYDaACIAvGTj+z0ADABQBEASeJBpDR6IXgHwjRCiDhHdAVvl54mosPl13rohIURBAPUALMlJWCHEcwAmAOgEwN/c50KJ5woAdQEcA1AKwBcAZgshBBG9B2A7gCFm2YY40KZBHsFQ9AZqUwrAFSJKt/HZBQClpTRCRIuJ6DwRmYjoVwAJAJ4xf9wJwCQiOkNE18AKNjs/EdERIkonojQiWkNEJ4nZCmA9eCCRQnHws3Qhl2O6AfiBiPYR0QMA7wCoJ4QIkthHEhHNJKIMAHPAg0UZid81yOMYit5Aba4AKJVDBIw/gMtSGhFC9BRCHDAvfF4HUA08iABAWQBnrA5PstGE9ecQQrQQQuw2L6JeB9DSqj17pIBdO/65HFPWWg7zDOYqgHIS+/jP6rt3zf8tLPG7BnkcQ9EbqM0uAA8AtLd+UwhRCOwi2Wp+6w6AglaHPGF1bCCAmQCGAChJRMUAHAYgzIdcAFDe6rsVbMiRmZ9bCJEfwFKwD7yMub21Vu3lmsvbrHh3AeiQy2HnAQRa9VkIvGB7DnyuQA7nKwEj17hBrhiK3kBViOgG2Kf+rRDiBSFEPrP7YjHY2p9nPvQAgJZCiBJCiCcADLNqphBYuV0GAHMIYzWrzxcBeF0IESCEKA5euM0NXwD5ze2lCyFaIOtawUUAJe2Ef44E0FsI8bYQoqRZrppCCIsffj6AV4QQtcwDy3gAe4gokYgugxV+dyGEtxCiD4BKdmS25iKAig4cb5DHMBS9geoQ0RcA3gVb0LfAUSQFATQzL34CwC8ADgJIBPvLf7X6/j8AvgJb0RcBVAdH0FiYCeAP8/f3IdvCrw15bgF4HTxApADoCo6KsXz+L4AFAE6ZXUVlbbSxE8Bz5tcpIcQ1ADPAMwMQ0SYA74NnDhfAiryzVRP9weGmVwGEAdiZm8zZmAzgJSFEihDifw58zyCPIIiMWZ+Btpgt2I8ANCCiZK3lMTDwNAxFb6ALhBA9AKQRkSMhhwYGBhIwFL2BgYGBh2P46A0MDAw8HF1k8ytVqhQFBQVpLYaBgYGBWxEfH3+FiOxuMtSFog8KCkJcXJzWYhgYGBi4FUIIW5sBH8Fw3RgYGBh4OIaiNzAwMPBwDEVvYGBg4OEYit7AwMDAwzEUvYGBgYGHYyh6AwMDAw/HUPQGBgYGHo4u4ugNVIII2LULuH9fa0kepVQpoEYNraVwH/bvB1JStJbiUQoVAp55BhDC/rEGqmEo+rzERx8Bc+YATz6ptSSPcvgw8N13QIfcancYAOB7OGoU8NRTWkvyKAkJwLBhwIgRWktiYIWh6PMKc+YAP/8M7N4NlNFhqdF9+4DnnwcCAoC6dbWWRr9s3gyMHAnExgKhoVpL8yhnzgD16gFBQcagrSMMH31ewKIc1qzRp5IHgDp1gB9/BF58ETh1Smtp9Mk//wBdugALF+pTyQNA+fLAypXAa6+xUWGgCwxF7+lYlMOvv+pXOVho3Rp4/32gZUvg2jWtpdEX//0HtGoFfPkl0LSp1tLkTp06wE8/Ae3aASdPai2NAQxF79lYlMNXXwFNmmgtjTQGDWKF364d8OCB1tLogzt3gJgY4JVXgB49tJZGGq1a8aDdqpUxaOsAQ9F7KtbKoXt3raVxjC++4Cicvn05Uigvk5EBdOsGhIWx4nQnjEFbNxiK3hNxZ+UAAF5ewC+/ACdOAGPHai2NtowYAdy8CcyY4Z4hi198AZQuDfTpYwzaGmIoek9kxAjg1i33VQ4AULAgL+rNm8eLtHmR//0P2LABWLYM8PXVWhrnsAzap04BH3ygtTR5FiO80tOYPJmVw44d7qscLPj5caRQ48ZAhQpAVJTWEqnHihXA55/zfSxWTGtpXOOxx/h86tUDKlZkd6KBqhgWvSexYgVPldescX/lYKFqVWDRIo4cOnxYa2nUYe9eoF8/YPlyjkf3BCyD9ujRwMaNWkuT5zAUvaewdy/Qvz8re09RDhYaNwa++YYX9i5c0FoaZUlMBNq2BWbNAp5+Wmtp5KVqVWDxYqBr17wzaOsEQ9F7AhblMHMmEBGhtTTK0K0bW7kxMRxR5Ilcv857CEaN4vvpiTRqBEyaxGGXnj5o6whD0bs7FuUwerTnKgcL770HVK/ObpyMDK2lkZfUVE4ZEB0NvPGG1tIoS9euPPts3Rq4fVtrafIEhqJ3Z1JTgfbtgebNgddf11oa5RECmD4duHuXE2d5SrgeETBgAFC4MPD111pLow7vvQfUrOmZg7YOsavohRA/CCEuCSEOW71XQgixQQiRYP63uPl9IYT4nxDihBDibyFEHSWFz9NYlEPRorzzNa/g6wssWQJs2cIRRp7AJ58AR44A8+cD3t5aS6MOlkH7/n3PGrR1ihSL/icAL2R7bzSATUQUDGCT+W8AaAEg2PwaAOA7ecQ0eASLcpg3T3Pl0KdPH/j5+aFatWqZ7127dg3R0dEIDg5GdHQ0UuTMnV6sGEdwTJzIkSnuzNy5wA8/AKtWcS53jTh27Bhq1aqV+SpatCgmTZqkbKf58nneoK1XiMjuC0AQgMNWfx8D4G/+vz+AY+b/TwfQxdZxub3Cw8PJwAHmzCEKCiK6cEFrSYiIaOvWrRQfH09hYWGZ77399ts0YcIEIiKaMGECjRw5Uv6O9+4lKlWKaM8e+dtWgy1biEqXJjp8WGtJspCenk5lypShxMREdTpMTCQqW5Zo2TJ1+vMgAMSRFB0u6aBHFf31bJ+nmP9dDaCh1fubAETk0OYAAHEA4ipUqKDCJfEQtmwh8vMjOnJEa0mycPr06SyKvkqVKnT+/HkiIjp//jxVqVJFmY5XriTy9yc6dUqZ9pXin3/4Pm7apLUkj/DHH39Q/fr11e00Ls69B22NkKro5V6MtbXf3qbzjYhmEFEEEUWULl1aZjE8lKNHgZdfBhYs0Gd1ISsuXrwIf39/AIC/vz8uXbqkTEcxMcC773K4nh5L69ni4kWW94svgOee01qaR1i4cCG6dOmibqfh4ezCevFF4PRpdfvOAzir6C8KIfwBwPyv5Sk+C6C81XEBAM47L55BJjpXDpoyZAhXp2rfniOR9Mzdu0CbNpxuuFcvraV5hNTUVKxcuRIdO3ZUv3PLoN2ypfsM2m6Cs4p+JQDLr7QXgBVW7/c0R99EArhBRMauCFe5e5cfgp49dakcbFGmTBlcMG+IuXDhAvz8/JTt8MsveZG2Xz/9RnBkZHDK6JAQ4MMPtZbGJr///jvq1KmDMlpVIhsyBGjRggdtI7WxbNhNaiaEWACgCYBSQoizAMYC+AzAIiFEXwDJACzD/1oALQGcAHAXgPtkL0pP55J7aWlaS5IVk4m3/1et6lYpe9u0aYM5c+Zg9OjRmDNnDtoqvZnL25sjkBo3BmrXBsaNU7Y/Z/jlF7ZUFy7UbVbRBQsWqO+2yc7EiUDHjkBgIGdg1VvIaf78PKv2cp9tSHYVPRHldNcfSSVoXhwY7KpQqkMEDBwI7Nypvzwx6elcCPrdd3WrHLp06YLY2FhcuXIFAQEB+OijjzB69Gh06tQJs2fPRoUKFbB48WLlBSlYEChQAIiL4xS/PjpKzkrEWUVjYnSbVfTu3bvYsGEDpk+frq0g3t6Avz+7K7/8EihSRFt5snPiBN/HL7/UWhLpSFmxVfqleXjlhAlEtWoR3byprRw5MWMGkRBE8+drLYm+eeMNIm9vDrvUI1u2EHl5EY0dq7UkmhCXeI2mbE6guMRruR/43Xf8e//1V3UEc5SrV4lCQoimTNFaEnnDK5V+aaroFywgKl+e6OxZ7WSQwltvsZLYvVtrSfTJ1KmsHBYv1lqS3Jk5k+WcO1drSVQlLvEahYxZS0+OXk0hY9bmrOw3bODf+SefqCugo5w8yWG9q1ZpKoZURa+jua0G/Pkn54jZuBEoV07WptetW4c0Of39jRrh6a1bUfrZZ7F5+nTcL1VKvrbdnFL79yPyww/xb7duOJE/P+8y1StlyiC0bVtU6tkT2y9cwI2QEK0lUoX1Z4EHaQIEgdS0DPy09k+cD8h6TKFz59BkyBCce/ZZHKhZU5X7mC9fPrzwQvaN/xKoWBH47TdOzPbHH0AdnWd7kTIaKP3SxKI/fpyoTBmiP/5QpPmVK1fK32hGBlFYGFHJkkR37sjfvjvy779E+fIR9eihtSSO0bo1Uf78RElJWkuiChaLvmJOFn1KClGRIkR166oql8vP6dKlROXKaXYfodGGKffg8mWO1f3kE8786C54efFCo5cXWxAmk9YSacv161ycIzwc+PlnraVxjBUrgOBgjhC6e1exbuKTUjB1ywnEJ2kblx4eWBzz+kViePMQzOsXifDA4g8/TE/nTJbFivEsWyXik1Kw/ixcuzbt2wPDh/Melxs35BNOZvKeor9/n3ffdezIObHdjQIFgL//BpKTOd44r5KeDtSowcph+3atpXEcLy+uCubjw8pegUE7PikF3Wbtxlfrj6HbrN26UPaDm1bOquQBLkaSksK/a5UipSzXZnWScP3avPkmh/W+9JL+wrPN5C1FbzLxhqPy5YFPP9VaGud54gm2fDZtAoYO1VoabXj2WbboVVQOjmLXmi5QADh4EKbkZCTXbSy7It596ipS000wEZCWbsLuU1dlbd+CrfOUPJPo2ZNnqXv3qlrn2HJtCML1ayMEV83Kn5/DtEl/G/b0+YQoxbvvAufO8eKrG212sEmdOsCvv/LMJDQUGDRIa4nUo0cPID4eOHRIt0XQLRZjaroJvj5ej7orLMc9yI9Pu32BxT8Mw7aXegJLfrZ5nDNEViwJXx8vpKWbkM/HC5EVS8rSrjW2zhOApHPHuHG8yW3dOt4trCKWa5OaloF8Pt6uXxsfH94I16gRMGEC6xod4ebazgGmTweWLeP85QUKaC2NPHTowA/L0KG88p8XGDeOC3SsWaO6cnAEqdb07lNXcbB0RQxuOwpd49bg9lfy5YDP1S8uE7bOU9K5L14MvP8+MGUKl09UGcu1aRVI8l2bwoWB1atZ18yf73p7MpI3LPp16zh9wJ9/Ap4WlvjOO8C//3KY199/s3XvZsQnpWD3qauIrFgy9wdu0SJWDlOnaqIcLFjkLV7QFyl3U23KLdWathy3sWoDfNO0J4ZP+QRoVZ+TtDkpl7U84YHFFVHwFnI6T+v3ihf0xdQtJx7KFRfHJQSHDmVXh0aEBxbH+QDIe33KlmVlHxXFLuJnn5WvbRfwfEV/8CD7AZcvBypX1loaZZgzB0hIAJ55BkhKAkqU0FoiyUh1cSAujotKv/66psrBWl4TcV7u/PkeldtiMdobwLIc99o3EB+QU4O25OsoMzmdp+W94gV98fHqI5ly/RoTiJpNGgHNmnluVanq1dkl1bEjsHWrLmaenq3oz57lh2bKFKB+fa2lUZZt24BKlThM7fRp3S5QZsfWNP8RBXX+PFtG0dG86KUh1vICXGwhJ7mlWtNZjpszh3OpODhoS7qOCmHrPC3vTd1y4uGgeP8+qjR/lpOVrV2rimyaER3NbsZWrYBduwCNa254ro/+1i1W8kOHAp06aS2N8vj48Ozl5k2gQQOtpZGMZervLWDbxXH/Pg9eQUHsl9cYi7xe5vxyXshBbgewRKjM35PMkSo/L2cFX6MGh5E6IFeO11ElskfbZMpFJqz8aRhMJhPil3lAMIQU+vYFOnfm+gP37mkqiiAdhAJFRERQXFycfA2mp3N2uQoVgO+/1yTr46pVqxATE6N6v0hIAMLCeNo4b576/TtBjj56k4mnwRcvAmfOAI89pp2QVkjx0TvSli1X0IJOoajdsCanp96zxyG5XJHHFXJyH8UnpaBghxfx5IFdaPrqTKSU8FPNtWQPxZ9TIqBbN46v//VX2Qc4IUQ8EUXYO87zhlUiYLA5U/LUqZqn9l23bh1CQkJQuXJlfPbZZ8p3GBwM/P47h3p9/LHTzTxiZSqw2cbSBwDbG2liYoBTp4ADBzRV8mfOnEHTpk0RGhqKsLAw/Ln8ZwxuWhkrvn4bM998CX3bNkVQUBBq1arlcNs5uYJ2XkkH4uJg2r8fx6NiJF3/HDckqURO0Tbhkz9ByP4/0aXrBFwoUkrRmP7cuH//Pp555hnUrFkTYWFhGKtGfQchgB9/ZGNl9Gjl+8sB93DkOsLEiWwBbd+uuZ86IyMDgwcPxoYNGxAQEICnn34abdq0wVNK13uNigKmTeNFy5AQrjPrAFIXHF3B7uLh8OEcMrprFxAQkHNDKuDj44OvvvoKderUwa1btxAeHo7o6Gj8+uuvmceMGDECjz/+uMNtZ8Zzm6+1tSsoHsC3nT7E7PljMLXnUODnb3VhBeeEzQicGTOASZNwevIMHL1UDt4KxvTbI3/+/Ni8eTMKFy6MtLQ0NGzYEIULF1Z+5p0/PydAq18fePJJTYIJPEvRL1rEC687d+qiWMFff/2FypUro2LFigCAzp07Y8WKFcoregB49VXg2DGeNlasyDlhJLL71FU8SDNlVnXPbcHRWXJdPJw+HZg0CScnz8S628URmZSiqYLz9/fPLHRepEgRhIaG4ty5c5n3kYiwaNEibN682eG2raNWsruCpm45gW3la2JM80EY98c0bJgdAXw8TNZzk5NHInBO7mel9uGHqDS0H+Zp7FoSQqBw4cIAgLS0NKSlpUGoNeMvWZIXoBs25MXoli3V6deM5yj6nTu53uSGDZpbgBbOnTuH8uUf1koPCAjAHon+Vln4+muOsW/UiH33Eq9L8YK+sF65EXBsgU+KrzjHOPNNm4BBg3Bu2Ci0ulQWqeuPqRouaI/ExETs378fdevWzTxP38vHUKZMGQQHB9v9viOx7pZrtKh2C1RJOYfe40YAMQ0cGrTVJvNcEhKAF17g2eQHH2T9TEMyMjIQHh6OEydOYPDgwQhRM/SxUiXetNmmDbB+Pec4UgnPUPQnTnAWuTlzOEJDJ9ha6FbNgrCwejVQrRpQqxYnQitY0O5XUu6mwksg023TMLgUhjWrIukhlRrPbTP+OiGBE7V17ozlMX2Ruv6YJuGCOXH79m106NABkyZNQkJKRuZ5Xls/Da+0sF8T19FYd+trVOO1HyAGdXd40NaE69eBiAj+zelsh6i3tzcOHDiA69evo127dlkMMVWoVw/47jtef9q1izdVqYD7L8ZevcrToI8/1l02x4CAAJw5cybz77Nnz6Js2bLqCuHlBezbx4tC4eGSsiRah+rlz+clWckDD90+JgJS0x5ddLMOvzv23y3sPnUVx/679VA51K4NzJuHyIol4eMlIAB4ewnNwgUtpKWloUOHDujWrRvat2+f6XrKyMjA7X93Ii0w0u6itTNJxrIssK5axW64WrUUTW3sEpaUw0WK8CxbpxQrVgxNmjTBvn371O/8pZeAYcM4xv7mTVW6dG+LPjWVUw63awcMGKC1NI/w9NNPIyEhAadPn0a5cuWwcOFCzNfCwjFnSUSlSkCZMrxTOBfCAcTevI8LN+7B//HH8MT/Nkru6rkLN5HvxJXMv59NLgWsLgoA+O/mfRw8dAH5TYT95hnDswBuE+H+ie0oULQosGPHw8aEAECaR04REfr27YvQ0FAMHz4cwMPB8MapePiWDMD6pHSsPZW7m8nlJGNeXpzMrXx5jrMfNEjza/MIK1cC167xZi+dbdq7fPky8uXLh2LFiuHevXvYuHEjmjZtqo0wI0YAJ08CxYvznh8JM21X0NedcJTr14EjRziPjQ7x8fHBlClT8PzzzyMjIwN9+vRBWFiYNsKULctKPimJQy/tKIgnzC9rUtNNSM0wwdfbC74+tieD5R+ko839NBDY7VP4dD4gP//MCj9IR0vzZ1kggu+d67ja+DmUNCuH3aeuIj2DF4QzMpx33cgRW75jxw788ssvqF69emYI5fjx4zGvXz0MfnUmird8EccyHs5iJm08bnMWJDUtQq4UKMCZS9evBxYsALy9nTonxfjvP55x6DANx4ULF9CrVy9kZGTAZDKhU6dOqK2inzwLQvDu2e+/B/bvV3yTo3srej8/Xtzo1AnYsoU3CumMli1boqXKK+w2ad+eH8LERF71dxCp/uWVe5Lx7m+HMv8e3646utatAAA4Zm4jLd0EIYB0Ky9Sk5N78cOSj3HmnY9QfsJYh6zfnJS5q/lfHrYbZnO9BQB2r1v2sJ80E0wAdpy4gr2J12z25/KC5Lvvcprt7ds5gkNm7A2MdgfOI0dY0ffpA/zwg+zyuUKNGjWwf//+LO+t0qq+8N69HBm3Zw+nu1AY91b0ANCkCUeXtGoF7N7NRTkMsvLOO1y6butWp5Q8ID2XyuHzN3L8O7tFe+y/W5ix7SSSrt5FbKWn8WlUP7z/+UfAMzUQ3q6dJOs3N2XuSv4XRwYJy3lN2ngcO05cUW4Bec4c4LPPeAOOQko+t3OWdE3Cwvi31ro17+odOVJ2Od2exER2Oc+cqYqSBzxhMRYAunfnvBKtWwN37mgtjb6YMwf4/HOnlYNl8bR4QV9JuVSyO4Sy/229uNi1bgV81akW8ufjdufXa4eEdt1geqkjjvy+TdJOT2tlbnGbPJJnJReZc6qONGnj8SyLykv3nc11sTU8sDiGNaviVL4ZSdWY/vyTreTRo7lKmgJkHxizn3Nu1zoLLVtyZsrRo3mjkMFDrl9no3TkSFb2KuGSRS+EeBNAP/CemkMAXgHgD2AhgBIA9gHoQUSpLsppnzFjeLt8167sztGb71ILtm9/qBzsLMDaIrsF90HrMLu5XdrXCcDi+LOZLpf2dXIPA8y+YagNumJ2/N+o06Y5Du46iJoRucc5P6wUZNttktusILfqSPfTHvqVTAAWx51BhokcDxm1gyQrOSmJdzu3bQuMH2+3TWexdpd5ewksiT+L9IyHctm71lkYOhQ4epRzLu3dq2rMuG5JTeViQVFRwBtvqNq10xa9EKIcgNcBRBBRNQDeADoD+BzAN0QUDCAFQF85BJUgEO+ovH2bt8/ndRITOef3iy86rRyyW3gpd1PtWtjhgcWxoH8k3no+BAv6S/OJWyz3lLupSE03oVunT3GxSAlUbt6QHw47353XLxINgktlxv5bW6NADnl0bJyfdXUkawSAtAxiS9ZOWKSj+Wbshlzevs0+79BQNmBkJPtMwroiVceI8kjPeNTtZeta53g9pk1j12qDBrw+lJch4sjAQoWAb75RvXtXXTc+AB4TQvgAKAjgAoDnACwxfz4HgHrzE19fYOlSXqzy1KIGUrh9my2op57i6+Ekzqa+dTa5VmZ/3l5o1/9b5Delc2y9ndj/7G4TizX61fpj6DZrd44uEVvnl5mG2HyMlwC8rZ4SE/HOYbnI9RqbTKzkCxQA/vpLtj6BhzOJ7NfIcu/a1wmwKZfDLqr163lzV82adgdtj+bTT4HDhzWLlHLadUNE54QQXwJIBnAPwHoA8QCuE5ElifZZAOVsfV8IMQDAAACoUKGCs2I8SrFinFOifn1eeFTRD6YLrJWDi+kWZAkHlIB1JId1fz6DDnBStg4d7Pp6rWU9f/0eFvyVbHdRVEp1pJS7qTh//R7m70kGgS2jlLvyKaxcr3GzZsCFC1xIxle+wQWwv1Cdm1wO/S68vDj7aPnyPGgfOJA3ctFbM3cuMHs2B4sUKqSJCE4reiFEcQBtATwJ4DqAxQBsbU21GZdGRDMAzAA4H72zctgkMJBX/lu0AMqV03VuENmJipJVOSidn8SWj3pwU0vJx+Kc+6ZRI15nsJPm2SJrfFIKlu47Kyk0M7fqSNYySm3PGWxe41df5aph+/ZxGLHMSAlfze3eO/S7KFiQY8WrVOEw3+XLXRHdvYiNZVfyli2aRgS6shjbDMBpIroMAEKIZQDqAygmhPAxW/UBAM67LqYTRETwKNq2Le+2fPJJTcRQlVdf5eiM+HhFlIMS2A2BbNCAI4Z692br/pVX7LYp90xErZlNJpMmcejdb79xlSkFUP2cKlR4OGiPGsWRYJ7O0aOc1G3BAs33+Lii6JMBRAohCoJdN1EA4gBsAfASOPKmF4AVrgrpNG3acMRCq1as7Itrn/1QMb75hpXD8uWKKQclyG5ZFi/oi6lbTmRVPj17ciKvfv04jUOjRnbblXsmolrmxdWr2QKcOJGNFAVRPZtkgwYc7tuzJ1v3fftqXhVLMS5eZL3z+ec8y9YYV3z0e4QQS8AhlOkA9oNdMWsALBRCfGp+b7YcgjrN0KEcdtmhA7Buney+Tl2wahXnzpg4kQc3NyJ7eOXHq4/YDjX85BNOuRwdzf964gzt8GHO29Rs8NzXAAAgAElEQVS3L99PT6R7d66TMGAA/i1cBt0Oezu9c1m33L3Lz2GPHjwT1QEurYoQ0VgiqkpE1YioBxE9IKJTRPQMEVUmoo5E9EAuYZ3myy+BokWB/v05zMmTOHyY/Z4eoByOnL+Re6jh4sU8Ba5ThyOLPIkrV4DISA4imDlTa2mU5ZNPgPbtEdy9PUpfOe9QNk/dk5HBg1mVKsCHH2otTSZ5Y/nb25sLZR896lIdVd1hUQ4NGritcrAO81scdwY+3nbC9nbv5vqxNWtKSrmcU59y1MGVqx2kprK7rUwZXrTLCyxejPshoVjz4xso+uCuZuUFZWfkSM7eOWuWrjKLun+uG6kUKsQujshIIChIsW3kqmGtHJwoYacXrBdjM0yEl58pj3LFHsvZZ+vrC/z9N9/D557jqAYHcDXRmdztAOBiFPfuAceP56nQw0L79iK1fCC2LxyOEzsPuL/bZsoUDu3euZPrxOqIvPOrAlgprl3Lo667W06RkawcDh7UXDn06dMHfn5+qFatWuZ7Bw4cQGRkJGrVqoWIiAj8lcOGn+wbhjrUCbC/2apUKbbsd+xgd5wDOFP8I6d2ciuwIpmOHdn9tm8fYK5nqgXHjh1DrVq1Ml9FixbFpEmTlO3U1xe+Rw7h8VspCO/dXtm+lGbVKt6BvnatLoM+8paiB3gr+cKFQOfOwD//aC2Nc3TsyOlgNVYOFnr37o1169ZleW/kyJEYO3YsDhw4gI8//hgjc8hiaL3t3iGruFo1Dj+cPRv46ivJsjq72zc71nV1TXByt+z773Nagw0bNF9cDgkJwYEDB3DgwAHEx8ejYMGCaNeunfIdWwbtnTs5qsodiY/nnFLLl2t+H3Mi77hurGnalCNUWrXiuo3ulNp4zBhWDrGxuvlRNWrUCImJiVneE0LgprlM2o0bN3Itoeh0mF/r1rzQ/tZbQHCwpIgjueLHrevqegkndsvOnQuMG8cDlYRwUTXZtGkTKlWqhEAnU1o7jGXQbtOG90q8/bY6/cpBcjKHwc6YoVrKYWfIm4oe4Fje06e5SG9srGZbkx3il194ejh7NvDss1pLkyuTJk3C888/j7feegsmkwk7laofOnw4+7Y7dGDLSsIeAjnix10qC7hzJ68Rvf22pA1garNw4UJ06dJF3U5bt+a6EsOHc8SKwnsIZOHGDU7JPGIEh8XqGSLS/BUeHk6aYDIR9exJ1LYtUXq6rE2vXLlS1vbozz+JvLyIRo+Wt12ZOH36NIWFhWX+PXToUFqyZAkREf36668UFRWlrABNmxIVLEh08aKy/VgRl3iNpmxOoLjEa9K/lJRElD8/0YsvKieYCzx48IBKlixJ//33nzYCvPYakbc30cGDqnTn9HOamkrUrBnRkCGsRzQCQBxJ0LGaK3nSUtETET14wEpi2DBZm5VV0Scm6lo5ED2q6IsWLUom8wNgMpmoSJEiROSkcpRCRgbdC3ySbpcoTXHHNVJS9rh9m6hECaIaNYgyMmRvXo5ru3z5coqOjpZRKidQcdB26jk1mYj69CGKiZHdQHQUqYo+77puLPj6ss+7Xj3gxAnA31+WZmskJfFKvKsQ8eJxSIjs+ciVpGzZsti6dSuaNGmCzZs3Izg4WN6QxGzEn7mBfi9PROzkXqhWNQCXO3VF6SL6CnHDypVAvnxciEPmSCm5ru2CBQvUd9tkZ+NGdt+UKcOLnAqm9XXqOT1zhlMcbNvmNgWODEUP8K7ZSpU4rW+RIrI0WerOHR44XMVk4vKI9evragOGNV26dEFsbCyuXLmCgIAAvPL6SLQYOBYDhw6DrxehQIECmDFjhtM1XKXkQ9l96ipueBXAuSKl8NSVJOTfsgkopCNFT8Qb3J59VpE0HK7Ux7Vw9+5dbNiwAdOnT5ddPofw8uKkhCdPcj57BdOWOPWc3rzJqTjcYV3PjKHoAeDdd3lhJTmZ87jLwOZVqxATEyNLW1i6lEMqq1UDBg+Wp00ZWbBgQaYyzsxXk2KC74sT8IO1ZZmU4vACplRLNbJiSXy9+msEXzuL5wdOx/hRHfW3AWffPo7MGDoU+PZbWZt2aXHYTMGCBXH1qg7SEIwbx+ku1q9nhaogTj2nt28DjRsDEyaw7nADDEU/Ywa7RHbtkk3Jy06HDvzjf/11DiNs3lxribJgrYy9hECGiUB41LJ0JrRRqqUaPnca6vwTixVf/IjxL8XoT8kDnKPn11950K5aVdZBW/W0w0qxZAnvL5g6VXEl7zSFC3OW0chIDnHW2tUlgbyt6NetA8aO5SLaJXWeZ+OddzhrY+vWvBs2NFQzUbK7UqyVMUDw9hIgIpuWpXVooxSXjCRLddEi4P33IaZOxYsDHS+CrhQ2z69DBw6Rff11dhe+8IJs/amedlhu4uJ4I+PQocDAgVpLkzv+/qzso6K4VKLOw53zrqI/eJBj6ZcvBypXtn+8Hpgzh/2JzzzDefZLlFBdBFuulOzK+IPWYUi5m5qrApfqkrFrqcbFAV27Am+8oSvlkOv5jR7Ng3ZMDOft0XDQ1g3nz/PGseho96n3XL06MH8+z9C2buWACZ2SNxX92bNsGU+dyouc7sTWrWwJ1qzJG7581L2Ftlwpg5tWVswlA9i2VOOTUvD3X/+gV49m8IqO5sIrOsLu+f30ExdT0XDQ1g3373Od48BAYM0araVxjGbNeIZm2WVfurTWEtkk7+W6uXWLlfzQoTwSuxs+PmwF3rwJNGyoevc55YoJDyxuPxmZhHakEJ+Ugle+34Y2vVridOHSiJ8216lzURJJ57d1K7sMa9YE0tPVF1IPmExAeDj/Gx+veYI+p+jTh11ObdtyokEdkrcs+vR0oFMnXkRxp3wa2Xn8cXZZhIUB3bpxrn2VkGvRz5V2dp+4jCWzhgIAWr8yGUMSUxD+pL7WWCSdn2XQLl+eawrs2aO+oFrTpg1XgEtI4CLi7sonn/AMu2dPXnDX2YCVdxQ90cMohylTdBuTLpngYOD33zkCJyQE+OAD1bqWa9HP2Xa6jRuCAtcvoumrM0EFCui2YIWk8yta9OGg3bUr+3zzCsOHc0DErl28oOnOCAH88AOvMYweDXzxhdYSZSHvKPqJE9li2r5ddb+2YkRFAdOm8SJkSAhXnPd0hg9HsW2bcXTZOnQvEuTeoYQWsg/aY8dqLZHyTJ8OTJoELFgAPP201tLIQ/78nIWzfn0Ou9RRcICHaDw7LFrEVvzOnbLtfNUNr77KxZa7dQMqVvSch8YWVsohtE0zuGOsSo4hpVFRwHffAa+9xsq+c2fthFSazZuBQYO4pqqnGSclS3LxkYYNeXG5ZUutJQKQFxT9zp3AkCFc3MHdp4c58fXXrOwbNWJfpyee56ZNbq8c7IaUDhjAYZfdu/OgreP85k6TkMB7Bzp3VtXdqCqVKvEmzDZteHdv7dpaS+ThUTcnTgDt2wM//8yRDZ7MqlWsHGrVAu7e1VoaeUlIAFq0cHvlIKmM4ddfA88/D1OjRvhx0XbXC4/rievXOYdNrVqqBhBoQr16PEOLieEkaBrjuYr+yhWeNn38say7D3WLlxeHpwnxMFzNE7Aoh9q13V45SA0pjZ/yM04WKYO2vVuj77RYl5R9fFIKpm45of2AkZ7OxlaRIjzLzgu89BIwbBjH2JurrWmFZyr6+/eBF19ka37AAK2lUY8CBXjHb2Ii/7jcHYtyKFqUC4G7OVLr4+5OTEFM70kgAEtnDcXuE5ed6s/iKvpq/TF0m7VbW2XfuDFw7RqHk3pKMIQURozg0NmOHYG0NM3EcEnRCyGKCSGWCCH+FUIcFULUE0KUEEJsEEIkmP9VNyTCZOLybOXK8Y41DTlz5gyaNm2K0NBQhIWFYbIaW7vLluU82Rs2sDWhE5yyLBs3BlJSePCSQTnkJoM9+fr06QM/Pz9Uq1Yt872DBw+iXr16qF69Oho1ewFfrt5v9/ykbCyLrFgSKFAArfp8i3I3L6P7J4MknmFWJLmK1KB3b+Cvv/il8Q7goKAgVK9eHbVq1UJERITyHQrBmUq9vXmNicj+d5RASnWSnF4A5gDoZ/6/L4BiAL4AMNr83mgAn9trR9YKU++8Q1S/PtG9e/K16QQrV66k8+fPU3x8PBER3bx5k4KDg+nIkSPqCLBoEZEQRNOmqdNfLsQlXqOQMWvpydGrKWTMWmkVkHr2JMqXj+joUcVlkCLf1q1bKT4+PksVrYiICIqNjaW4xGv0ROthVKz+y9LPT4K8UzYn0JGVm7i03uuvO9VGyJi1VNGR6y4348dzCcx169Tv2waBgYF0+fLlzL9lL/mZEzdvEtWqRTRhgqzNQmKFKacteiFEUQCNAMw2DxipRHQdQFvzAGAZCF50tg+HmTmT81ivWKGLlMP+/v6oU6cOAKBIkSIIDQ3FuXPn1Om8Y0ferWeJONIQhy3LCROAuXN5gblqVUVliE9KwaSNx/EgLXf5GjVqhBLZrNFjx46hUaNG2H3qKvJVqIk7x3bKZjlbLP+nYp7jTVTffst7JhxsQ4qrSDGWLgXee4+TlD3/vLp9640iRTjb5bRpXDFOZVyZD1cEcBnAj0KImgDiAbwBoAwRXQAAIroghPBzXUwJ/PEH57Hevh0oVUqVLh0hMTER+/fvR926ddXr9L33gKNH2V9/6JBm2fUcKoqxZAnL/e23sioHaxm8vb1w7vo9zN+TzEVS0k0gsB/Tkbw71apVw8qVKxFZqxHuH9+B9FtXMtuOT0qRT7F26gQcP875mYKDHcrTrlnq4n37OAx20CA2NnSCEALNmzeHEAKvvvoq/GUqHSqJcuVY2TdrxiHQauaqkmL223oBiACQDqCu+e/JAD4BcD3bcSk5fH8AgDgAcRUqVHBt/nLwIFHp0kTbt7vWjoxYTwlv3bpFderUoaVLl2ojTGQkUZEiRCkp2vRPEgtXx8ezm2LIEMVkeHfZ3xT87hp6cvRqqvTOGgoatZoCR62mJ0evpu6zducqX/YC6EePHqXo6GiqU6cODRg2ih4r/Hhm24q4Srp1Y3fWv//K267cXLhA9NhjRFoXGbfBuXPniIjo4sWLVKNGDRo/frz6QqxbR1SmDNHx4y43BaVdNwDOAjhLRJZMTEsA1AFwUQjhDwDmfy/lMMDMIKIIIooo7Upqz3PnOBvl//6nSTZHe6SlpaFDhw7o1q0b2rdvr40Q27cDxYoBNWpoliXR7iLk+fN8/6KiZC+zZy1D2WKPId1EMBFgIi6S4i0AXx8vDGtWxSHrt2rVqli/fj3i4+MxYmAf+JULzGxbkcXPuXO5StXTT3PYqR65f58jpcqX5zw2OqNs2bIAAD8/P7Rr1w4JCQnqC/H88+xWbdECuOxcRJWjOO26IaL/hBBnhBAhRHQMQBSAf8yvXgA+M/+7QhZJbXHvHiv5QYN0uWWciNC3b1+EhoZi+PDh2gni48ORK+XLcxm0Dh20k8UWRDylLV+ec74oiDNFUnLi0qVL8PPzg8lkwqeffoqeffth0U3X6rba5c8/eWNc8eL8m9dZlkRs3szGxP79upPtzp07MJlMKFKkCO7cuYP169ejuVZlOfv356ydfn6850fhCneuxqwNBTBPCOEL4BSAV8CuzkVCiL4AkgEol/T97l3gv/90W9llx44d+OWXXzLDuQBg/PjxaKlF/ovHH+cFof/+Y/+pniDie1mliuLKwdn0yF26dEFsbCyuXLmCgIAAfPTRRziafAk//zADj+XzRpdOL+Gjt4agbfJ1SW1LKaNoEx8fHhDPnOENct7e0r+rBteucQoAHQRDZOfixYto164dACA9PR1du3ZFjRo1tBPoqaf436QkfSt6IjoA9tVnJ8qVdiVTsiRHZrRowQsdOssN0rBhQ8t6hPY0a8a78y5eZCtCb+zYwbl6Ro0CPv/c7uFOK0o4t0C5YMGCR/qfOGs3CnefAl8fL3TsFwkhMfW11DKKNnn1Vc7CevAgu+L0RnIyD9jt23OZTh1RsWJFHDx4MMt7q1at0kaYrVuBt94CDh/mFNUK4/5b1CIigNmzeSfsjh2cHtQgKwMGsJ9+3z59KnmAdw/+9BPQqxfP0Pr0yfFQlxSlTOQUrilFLkfKKGbhm284hHj5cn0qeQCoUIHdN88+K3nQznMcPcqRVPPnq6LkAU9Q9ABniUtK4jDCHTvYf2nAfPMNMGsWK4fq1bWWJnd69OAwwv79uWB7o0Y2D8uuKJfuO+tyxStHsRUyKlUuh8JNLaxezdvpJ07k37ueqV+fC9n37MnWfd++WkukHy5eZD31+ecceKASnqHoAY4xPnWKFxrXrQN8fbWWSHvcSTlY+OQTTtUbHc3/2pihZYmJ9xJYEn8W6RnqWvc5+fqlyOXwOsHhw0C7dqwwR4xQ+tTkoXt3Tp09YAAP2o0bay2R9ty9y89hjx6cFkJFhB58yBERERQXF+d6QxkZrOgff5zdABqWC1y1ahViYmI06x+HD3PGx969ebrvbtSpwzU4k5NtFoux+OjPX7+HBX8lw0SAtwCGNw/B4KaVNRBYIbmuXAGCgjgj6datssurOJ068U71o0c5WkhnqPacZmTwbvVChThtuky6SQgRT0R2k/boK/7JVby9OZXt0aOcntgDkZQc7MoVLoDeoIF7KnkA2L0beOwxzl1uI+WyJS6/fZ0ASal/1UJWuVJT2RdfpgywZYv8wqrBokVAtWo8cGucqldTRo7kiKRZszQxQD3HdWOhUCGOxImMZEuoVy+tJZINSYuQ1sph82ZtBJUDX19OaRsUBDz3HBAba/MwZ8MllUYWuerV470ix4/rLibdIXbt4vtYqxYXA3Lnc3GGqVO5vODOnVxXVgM884qXKcMXduRIt7GEpFjqkpKDWZTDwYPu/0CVKsWW/Y4dvECbAxYrGsAj11DLwhtSUhLnSMeO7H7bt483uSmAatfGMmhfugQ0bepyc/bk1k2xFYDXycaNY32kYZCI51n0FkJDOUtc586s7C2bE3SI1HBBu9EaL73EyuHffxVTDqpTrRrw22+8iFW1ao6LkbauISAt3FF3vP8+1xzdskWxcGHVQ1Qtg3bt2kC/fuzCcAJ7cush9PahMPFcG2P1as3Dvt3c5LND06YccdKqFe8I1SlS0/jmmnZ2zBhWiBs2aP6jkp3WrYGvvgLefhtYudLmIbauoW4KbzjC3LlsAc6alWN4qRzkdG3m70lGj9l7MH9PsvydWgbtH37g59IBLFb6sn1nM+VOtXFPrc8rNc2ESRuPa2PZJycDbdsC06cDamaszQHPtegt9OzJ0Rtt2rCft2BBrSV6BEfiqm3u6pw7l6tpzZ6tqHLQlDff5HC9Dh3YUsq2YSina+hwvLqW7NjBa0pvv82WoILYul7z9yTj3d8OAQC2J1wBAHStW0Hejlu35gLow4fzxjgJYb/WVrqXAEzmQEETAcULZg2jtpxXapoJJgA7TlzB3sRr6lr2N25wveoRI3iHsA7wfEUPAB98wDH23bpxvnOd5QdxaeHOohxGjVJcOajNI2kOvv+elX29ejx4m3f5Wo6zlaBMjwu1NklO5g00MTGq7CYNDyyOD1qH4ffDF9Cimj/CA4tj0sbjWY75/fAF+RU9wCUujx1jJbhvn91dvtZWunU0uBeAlLupWY61PEuTNh7HjhNXHN997CppaexCbdJEV6U884aiF4LDDF94gfNLfPON1hI9glMFIpKSWDm0acNVmTyIHH2tmzbxbsuaNYGkJMRfuJOrT9aVwhuu5NNxiNu3OSKlShX2zatAfFJKZtGVvYnXEPJEEbSo5p9pyQNAi2oKFuX47jubg7Ytsm+QgxDIyMh5lhYeWBzDmlXB3sRr6s7miIDXXuOEbpMmabqPJzt5Q9EDvPK/dCnHlv/vf8Drr2stkWvcvs0LWyEhfF4eRo75YLy8gAMHOINj3brY/fVi5/LG2EG1RT2TiWPMfX2BuDjVIqVsXV9L5JLFylfEmrdm48Ysg3ZOu9mzz3gt8uc2AGsSdjt+PP82t26VpZi9nOhLGqUpXpzDnOrX57hed0kLkB1r5bB3r+ZhlJMnT8bMmTNBROjfvz+GyTBlzXXdonBhfqBCQtDlizfxbcQg2S03pxOPOUrz5lw859QpVdN22Lq+ffr0werVq+Hn54dfDh8GALz88ss4cOgf3ElNh+n+HRQp+jhen/qbPMrTMmhXqMCZZw8cyPHQ7DMzKX2rWkZx/nz2Guzapc+INyllqJR+hYeHu1xSyyH++ouoVCmivXsV60LR6vJRUUQFCxJdvKhcHxI5dOgQhYWF0Z07dygtLY2ioqLouAwl0ogklB/cvp3Iy4vOD3rTfplCJ/oOGbOWKipVFpCIaOBALp24f7/8bUsg+/XdunUrxcfHZymXaLkOT45eTcWeaUclGnWXv1RiYiJR/vxE7drJ054DyPKcbt3KpUwPHXK9LQeBCqUE3Zenn+bRt21bnjK6EwMHcvTQjh26SDl89OhRREZGomDBgvDx8UHjxo3x22+/ydK2rQ1HWTbDNGwI/PAD/L+bhMHJO2S13nINZZWDyZN5cXnxYvbPa0D269uoUSOUKFEiyzGWmU2GiXDr6HY8FvKs/OGqgYHsxlmxAnj3XXnaVItjx3hz2/z5HD6qU/KW68aaF19kJd+yJSvNYsW0lsg+kydzXO5vv2mmHLJTrVo1vPfee7h69Soee+wxrF27FhERdnMsOYVNv3mvXvyw9enDlY1krBus2NT/9985XPSzzzgrpY6xuHhunvobPoWKoaBfeWRkmCAEsP7Ifyhe0FceX37DhsCPP3ISvpAQ90hdcvky648JE7iwj47Jmxa9hTfe4KiVDh04R4yeWbuWlcPnn/NMRCeEhoZi1KhRiI6OxgsvvICaNWvCR6GFqBw3QI0fz9ekWTP9z9COHOG1od69OUWHzrHMbCrdOoiBfXtiQf9IPBdaBukm4ODZG3j3t0Pyba7q2RN45x0etP/8U542leLePb6PXbrkWiRHL+RtRQ9wqGWhQlyiTQcpm21y5Agrst69eTONzujbty/27duHbdu2oUSJEggODlakH4t16S0Aby+B89fvPdz1uGwZp72oXZsjkvTItWu8SzIykneHugk1yxXB8T2bMOK1VxAeWBz30zKyfP774QvydTZuHM9yoqL0O2ibTJxTvmJFrp/gBuRd140Fb29gwQIufdagASdEk4GICxfkeZhNJuCPPzjeWKfK4dKlS/Dz80NycjKWLVuGXbt2KdKPxbpcuu8slsSfxYK/krF039mHPvQ9ezjsskgRHhh1FMcMgHPXlCrldnnlN27ciKpVqyIgIAAAlI+3X7KEB+ygIE5fki+fvO1b4dRzev48y7Rpk/5+YzlgKHqAU4cWK8YbN+7dk6XJgjduAA8euN6QycS77cwPmR7p0KEDrl69inz58mHq1KkormCWvvDA4th96irSM2yEPvr4AEWLcpbE06c1Dzt9hLt3udqS3uQy06VLF2zYtAUp166ijH9ZjP/0E/Tt2xcLFy5Ely5dMo/rWrcCkq/ewboj/+GFsCeUibf39+dwy9OnFQ07deo5vXaN18gUHIDkxlD0RMDgwazsExJk2+iwTc7KNZs2cbx11aqczkFnbN++XdX+coyxj4kBzp4FzpzR58CYkMDFoLt25SgNnfHWZ9Owf9ZuFDUvdtdqxhlAf/rppyzHxSel4KddiUhNN+GnXYmIDntC0qK15J3Gw4cD69fzDO2ZZ5w/IQk49Zw+eMClLkePBr74QhnBZMZQ9BMn8g9q+3bd7WbLJCoKmDaNQytDQoCXX9ZaIk2xuetx+HB2ce3apU8lDwDBwRxx07w538exY7WWKAtSN4k5s5lM8k7j6dM5fcCCBYoreafJn58j3+rX50yxAwdqLZFddKrZVGLRImDKFK78YqMuqa549VUOI+zWjReBnn5aa4k0JUvoo7Vy0Pt1iYriPC+vvcYzNB0N2lKzqDqSbdWCpMFh82Zg0CDgww9zvS6q5SDKjZIlORKuYUPeB9CypTZySCTvKvqdO4EhQzh/u14twOx8/TUr+0aN2A3gLnIryaZNkpSDrhgwgIvD6GzQlpofxpk8MnYHh4QETjr48su5uid1VVikUiWO9mrThl1NtWtrI4cEXFb0QghvAHEAzhFRayHEkwAWAigBYB+AHkSkryD1Eyc4ReqcOZxQyZ1YtYp34NWqxalts+XX14W1oxYJCUCLFlxFTIdrF7mi00Hb1iYxW78pRzeT5To4XL8ORETwb9rO2oVqOYikUq8ez9BiYthtWL68drLkghwW/RsAjgIoav77cwDfENFCIcT3APoC+E6GfuTh6lWeZn38MSsJd8PLi3N4ly8PhIdzjL05ikNKmTWPGQQsyqF2bWDePK2lcY5Vq4Dq1XMctPWAnBa0zcEhPZ3Pv0gRnmXbwRm3keK89BKQmMihoH/+yZFfOsOlOC8hRACAVgBmmf8WAJ4DsMR8yBwAL7rSh6zcv8+pD9q14+mzu1KgABf/tvy4zORWOs/ywH61/hi6zdqtj8LJzpKezjOxokU5fYW74uXF1bK8vHjQNpm0lugRFC/H2KQJG19//y0pGELxHETOMmIE++s7duRwaJ3hakDvJAAjAVh+oSUBXCeidPPfZwGUs/VFIcQAIUScECLu8uXLLoohAZOJKzD5+ytepCM+KQXrz0JZZVq2LEcKbdjAqRGQdedodmvH2Qc2SxIxvdC4MccyHzyo30gpqRQowPHiSUlZBm29kNtvyhXik1Jw9IUOMO3ZA/z1F5AtmVpu2Ep2pzlCcJ0LHx9eM9LZLnunFb0QojWAS0QUb/22jUNtnjERzSCiCCKKKF26tLNiSGfMGJ4ez5mj6IYVi+W8OkkobzlHRLBPc/Jk4LvvcrV2nHlgdTkL6NWLFYODykHXlC0LbNvGg/Ybb2gtTRaUsKDjk1KwtfebCFn/G/p3HIv4gk/IIKkO8PEBFptnkhgAACAASURBVC7kAjIqlIN0BFfMoQYA2gghWgIoAPbRTwJQTAjhY7bqAwCcd11MF5k5k9PB7toFPPaYol1ZLGeCUGexqFMn4PhxjiAKDkZ4s2Y2+3M0UiI+KQWTNh7HgzQTCDpZ+JowgQuhr13LeW10iNPrIBERrCQ6deIY+0GDlBPSQZzJ4pnbdbj403wMi/0ZY5sNQGxgbdTR+nclJ0WKAKtX8yJtUBAHCugApxU9Eb0D4B0AEEI0AfAWEXUTQiwG8BI48qYXgBUyyOk869cD77/Pbo5SpRTv7mEV+gzk8/FWZ7FozBgO12vZEjh0iBWFDaQ+sNYLcASe9mm+8LV0KfDee7zv4fnntZMjF1xeuHzpJU6SNXQob66KjlZOWAXJ9Trs24cWn7yBeRGtMD8iRvvflRKUK8fKvlkzjqaSMXW2syjhwxgFYLgQ4gTYZz9bgT6k8fffQPfunCRJoYyK2bFYzq0CyaEH3WVf+Ny5vKD39NMckeIC1v58LwE0CC6l7cLXvn0cXz14sKqWrqP3RJaFy/fe49S3rVpx+KUbkv06LNt3FlO3nMDBvf8CDRtCPPccQpf8rL8FVTmpUQP45RcevBMStJZGng1TRBQLINb8/1MAtN+7fP480Lo1L5CoPKKGBxbH+QBpdS0B1y3BzGny3JUIbxrOP7JTp5xeqIysWBI+XgJpGQQfL4Fhzapo9zBeuMD3LyoK+PZbVbqMT0rBsn1nsTjuDNJNJPmeWF83by/hsKWaeR/HfYvwkyd50E5Odo+iOFZYh0B6e3thcdwZeKU+QOdpr+D+E2VRYN06hHt5eaaCt+b553mG1rIlh46qsRaZA24espADt26xRTRwoG58ZLnhyiaQ7IPEgtVbUbthTU677Eq6YCEAkLZpWO/f5zDK8uU5R4wKWK6nZW0CcPCeOHndHhnsZRq0tcB6Pejc9XtYuDsRa396E14mE+ZMXYZXdZq9UxH69+cMnG3b8i5uhdcIc8Lzrnh6Oiv3iAjOLqcDMjIyULt2bbRu3drm566EsGUfJHZezQD27uX47B49nJLXkgaYAKRnmDBp43H1I25MJr6HJhOwf79qqX0fLqYzAnxP/vhuLPz8/FDNqi7oyy+/jFq1aqFWrVoICgpC26j6mdctI8Mx180jg33yDXY9Xr/OgzZ0GuqaA5YQyA51AjDzt3F4MuU8Xuz/LSKe0nbnaJ8+fR65j4rz6adAhQocMabRXgn3MROkQMQLWenpnO1RJ0UBJk+ejNDQUNy8edPm587kDrFgc6dgYHFgzRrOHVK1Kvt9HSBzQdmseP5MuIK9idfU9ae2bQucPMl+ahV3jGZ3O7wUHoAOdQJwJ6kQCo8egZ49e2Ye++uvv2b+f8SIEbhDvoh1ctemzftYrBgP2tWr42q7TuhW7RV95HhxgPAp40En92Lp5IX4OsZ2RJia9O7dG0OGDMlyHxXHywv46SdenH3nHW1CL4lI81d4eDjJwsSJRNWrE924IU97LrBy5UoiIjpz5gw999xztGnTJmrVqpXT7cUlXqMpmxMoLvGa9M+mTSMSgmjRIqf66z5rNwWNWk2Bo1ZTxdGracrmBGfFd4y33iLy8iLavVud/rKR0/U8ffo0hYWFPXK8yWSigIAAOn78eK73yZl+4xKv0fIvfqQMIejLZ7urfy9cYeZM/v3Nnau1JFmwvo+W51QVLl8mCg4m+v572ZoEEEcSdKznWPRLlvDGoZ07dZVrYtiwYfjiiy9w69Ytp9uwt1ibY9jkwIFsEXfpwnmzIyIk9xkeWBzDmlXB3sRr6uYVmTUL+Oorzl9Tt67y/ZnJHvftiOW5fft2lClTJrNWrly5YB7e99KIf34gPlz3HU6XKIeN1RujeEFfTN1yQr95i2JjObX2Bx9wlk4DDu+2pDauUEHVXFueoeh37WKltn69rrLHrV69Gn5+fggPD0dsbKzT7biUsW/SJFb2zz7LrpCyZSX364pLySk2b36oHKxK1ymNq1FPCxYsyFJqTy6s7/u8Wi3xvPd1/G/VRPzxYgO8ufqIft04J09yxEnHjpw+2uAhlSvznpAXX+Sd0LVqqdKt+y/Gnjz5MOWwzvJB79ixAytXrkRQUBA6d+6MzZs3o3v37g6343K+kTVreJdezZocyeIAquUVOXmS1xTsKAclFiRdiX9PT0/HsmXL8LICufCz3/cCU6fAq3lzRL/aEY+nXFYu0Zgr3LzJ+zmqV+edvgaP0qABMHXqw9KXKuDeFv2tWxyj+sEHuqzwMmHCBEwwJ1CLjY3Fl19+iblz5zrcjsuWtSW1cdmyHN4VHa2bhWoAvIi+fTuHEuaiHJQqOuFK6tuNGzeiatWqCFAgn7zN+75mDdKCQ7Bnai/sCKwBk7cPaux4HPhMJ4Wq9+wBChUCdu/WWhJ906kTZ58tX54T2lVQoMC6Fe6t6E0mIDWVf1gejjP5RrLg7c0vgK0uvSn6jAy791GpohNSBtIuXbogNjYWV65cQUBAAD766CP07dsXCxcuVMRtYy1bFnm8vPBYEb5OlQsQChX0QuG0e0DaPcVkcIi0NCBfPtXCYR0l+31s166d48XB5cISTebgLNsppKzYKv1yKermn3+I/PyINm92vg0FUHU1Xwq1axMVK0Z065ZDX3MlisQhDh0i8vEh6tcvV1lCxqyliqNXU8iYtcrLJBOyXsMOHYh8fYlOnXK9LSW4fJmoUCGiRo20lkQSmj2nK1cSPfEE0cmTLjWDPBN1ExrK0/3OnYEtW4CnntJaIv3RsSNXovr3X6BwYclfU7U+Z7VqwG+/cf3NqlW5kEM2VF8clgFZr+H77/M12rKFo6j0SKlS7LapXRvo14+jqAyyEh8P9OnDic8qVlSlS33OrxylaVNg4kROe3DxotbS6IsxY7iA8YYNNpVDboubilcXyk7r1hxa+fbbwMqVNg/RZdGJXJDtGs6dC4wbx4qzUSOXZFJ8h61l0P7hB+DLL5Xpw11JTubNgNOnqxo+7P4WvYWePTmnREwMx/DqsP6m6vzyCzB+PDB7tk3lYM/a1KQ+55tv8syjQwe2fGrUUL5PBZHlGu7Ywdvn336bq6S5gGqztNatuQD68OFAlSo8U8vr3LjBQSPDh3OkoIp4jqIHOPrm1CneoLFkycPFx7zIjh1A7965Kgd7i5uauUqmT+diKvXq8eDt5yfpa3osfu7yNUxO5sydMTGybJ1XakHbJsOG8aDdvj1Hfbn5oO0SaWmcsrhJk8zSn2riWYpeCK4m9cILwFtvAd98o7VE2pCUxMqhTZtclYMUa9PlaB9n2bSJLcGaNfl8fH1zPVzV9QQHsXcNcxygbt/mDTVVqrD7zQ5SBjrVZ2nff+/UoO1REAGvvcb1gSdN0iTizbMUPcAKYdkyoH59zl8+dKjWEqnL7du8EBYSwjvwckHXi5teXlw0u3x59mXu35/r4apaqjKS4wBlMgF16vDvOS7Obrii1IFOk3u+caNDg7bHMX48/5a3btUs3bRnLMZmp1gxzikxYUKOi3oeicnESt7Xl7MeSohl1sPiZo6Lg4UL8wPyzz92fZrZd5HeupeGHrP3YP6eZAUlfxTZqlI1bw6cOwccPChJMTqy6Kv6PbcM2g8eqLoAqQvmzQNmzOAIGwci3uTG8yx6C0FBwIoVHImzdq1DCb3cluhorqx1+rTbWE12LdHAQHbjNG4MvPsuW0c2sLZUb91Lw/fbTgEAtidcAQB0ravszkNJ52IDm66UQYM4oCAuDihTRlLfmiycO0LhwjwrCwnhQVuCK8rt2baN/fGbNwP+/pqK4pkWvYWnn2affdu2PGX0ZAYO5Knhjh269oNmt3glWaING3Ko3mefcU6jHLBYqkcuZM37//vhCzb7lpP4pBRM2njc4VBKywCVWT91+c/s11682KGEV4+046S1rmjoZWAgu3FWrOBB25M5doz3r8yfz+GmGuO5Fr2Ftm05p0TLlqwE3az+piQmT+ZIld9+Uy0bnjPYsnglW6K9evHD06cPUKlSrnWAW1Tzz7TkLX8ruVibvfygl4PJ5zIXa9euZQvws8+Adu0clsPVhXNVFrQbNgR+/JEjwkJC+L56Gpcusb6ZMIGLjegAz1f0APDGGxx22aED1x51E7eGJNasYeXw+ec8qOkYW9b74KaVpS8Ojh/P4XpRURzJERho8zCLm+b3wxfQopo/utatgKlbTmT2nZrG5RHlKnpuXX7QC0CDyqUcb/vIEb5/vXsDI0dK/pqcIaWqLWj37Cl50HY77t3jaLeuXfn8dELeUPQAb95o357znf/wg76SejnLkSOc1/qVVzheXufkZL07ZIkuW8YLzrVrc4x5DgtcXetWyOKXzyyPmGaCCcCOE/KVR8x+Xg4r+WvXeJEyMpJ/mxLJboF/0DoMKXdTnVb6qvr5x41jZW9n0HYrTCau01ypEvDxx1pLkwXBeXG0JSIiguLi4pTv6M4d3rDQpg3nDVGQVatWKZsV78oVXnCuU4cXfdwEWSzQ1FRWDIUKsZKQmCnR4kffceIKTAR4C2B48xAMblrZOTmyte3UeaWn87nky8ezTgeyPk7dcgJfrT8GE/FMwstLwETkkttF9U1ntWuza/XMGU2iUmR9Tt9+m9M0b9gA5M8vT5t2EELEE5HdSBPPXozNTqFCwKpVnBLAibzwuiEtjWOS/fw4OkNj+vTpAz8/P1SzWnS6du0aoqOjERwcjOjoaKSk8OKeLKF9vr4cdnjhgkM+UEt5RJeKuOTStlPnFRnJBsjff9tU8rktjlqHlFqUvCs5da5fv44JI/rj24Gt0f2F+ti1a5fDbTjMnj28kahmTbaI3ZVp01i3LF+umpJ3hLyl6AHgiSfYrz1iBEepuCP16rFyOHBAF3m/e/fujXXr1mV577PPPkNUVBQSEhIQFRWFzz77TN5O/fy4hOS2bbzrUCJyRafIQufOwKFDnNPHRp1ji2vmq/XH0G3W7keUvfW5fNy2mssD2BtvvIEXXngB//77Lw4ePIjQ0FCnT00ylkH7v//YjeOOrFkDfPIJL6aXKKG1NLaRksvY1gtAeQBbABwFcATAG+b3SwDYACDB/G9xe225lI/eWTZu5Dz2//yjSPOK5bnu2JHzkbuYx1puTp8+TWFhYZl/V6lShc6fP09EROfPn6cqVaoo0/GKFURCEH3zjTLtK8XYsUReXkRbtuR4yJTNCfTk6NUUOGo1VRy9mqZsTsi1SVfy3t+4cYOCgoLIZDI5/F1Z+PtvrkfQv7+q3br8nO7bR1SqFNGuXfII5CCQmI/eFXMwHcAIIgoFEAlgsBDiKQCjAWwiomAAm8x/64+oKI5UcafUxh98wGkN/vhDtTzWznLx4kX4mzeJ+Pv749KlS8p01KYNp6gePpx3H7oD8+bxYt306bxmlAOO1gp2xS126tQplC5dGq+88gpq166Nfv364c6dOw634zTVq/Nve9YsDpxwB86c4WRz33/PLjg9I2U0kPICsAJANIBjAPzN7/kDOGbvu5pY9BY++IDomWeI7tyRtVnZLfq5c9lynTVL3nZlIrtF//jjj2f5vFixYsoK0K8fW4SHDinbj6vs2sWW/FtvSTpcrQpfe/fuJW9vb9q9ezcREb3++us0ZswYRfu0yZdf8u981SpVunP6Ob1+nahaNZZXQyDRopdLyQcBSAZQFMD1bJ+l5PCdAQDiAMRVqFBB6euRMyYTUY8eRO3aEaWny9asrIp+505WDm+/LV+bMqOZ68aaxo25jN3ly8r35QxJSUT58xPFxGgtySNcuHCBAgMDM//etm0btWzZUhth+vdXbdB26jlNTSWKjiYaNIj1h4ZIVfQux9ELIQoDWApgGBHdFBLj04loBoAZAIdXuiqH01hSGz/3HE8fS8oThVH/6lXgiy9cb8hkAv76i3faydGeSrRp0wZz5szB6NGjMWfOHLRVYzPX5s1cRat0ac5eqoOF6kyIOCd75cocmaEznnjiCZQvXx7Hjh1DSEgINm3ahKe0Kss5YwbH2FevDjzzjKIbHJ16Ti9f5pDYyZPdZj+OS4r+/+2deXRV1fXHP4dEFJEKVsoUgjI0P3EWsAqKLkQGoYA2MoiC1kqLIFSoSKtLpWKx+mt/cQCtiAMFgkUGGVVEGayAQBSohBQIAmGQoEBs0ZJh//7Y72kKAV7y7n13eOez1l3Je4R393nn3u85d5999jbGnIaK/FQRiWYp+sIY00BE9hpjGgAuOWcdRARKSzWW2aGLqiw11ZnPioacFRfH/1ku0a9fP5YuXcqBAwdIS0tjzJgxjB49mt69ezNp0iTS09OZMWOG+4aUlWk/gn73fhJ6UFHwcT8+99xz9O/fn6NHj9K0aVNeffVV74yJfk/Vq7sq9FW6T1NTdf+DD/YgxUws0/6KDsAAk4GsY95/Ghgd+X008NSpPstTH31pqUifPiK9e+vvDuGo62bNGpGUFJFhw5z7zASTEF9z27YitWqJfPmle+eIh927RWrUEOnc2WtL/M2AAeq6cSkirjxVuk+Li0VuuknXhQLiuolnytMOuAPoYIz5NHLcBDwJ3GiM2YIuzjocQO0wDz+sq+evv+6/GWCU1q0hO1sLqbzwgtfWVJpTxYM7wsCBmoN/zRr/xjI3bKhx/++9p/mXLMczbpxuZpw/HxIRx18VUlPhjTc0jbQD5R0TQZVdNyLyITqrr4hg7HyYOFHTwa5cqbvz/Mytt+p2/6FD1c97442uni5QybKi4rBwoWZEPAVVaZtj30fr1jB9OvTurbbee2/VPytszJwJDz0Ezz4LnTt7bc3JOessHYyuvlpTkfTt67VFJyV5kpodyzvvaL6bFSvg3HO9tiY2HnoIcnM19n/jxphErSo4na7W1WRZUXF4/vmYxKEqbXM8fW9mpu6kvO8+aNHC9UE7EOTkQJ8+MGSITmaCQKNGKvYdO0Jamq+zcPrUV+EyGzZolrmZM/VGCxJTpkCrVlpU5dAhV05RmbJ0seBa2oGoOAwdGvPMuCptc/r7AHRw6t9fB+28vPg/L8js26ciecMN6p4MEpdcovdkZiZs2eK1NSck+YR+927o3l0vqHbtvLamaqxYoQVULrlEV/8dpvyOzJRqhj2Hvonbr+54ndK9e78Xh2efjfm/nWi3aazJwxx9Ipk82fVB2/d8+60mNGvcWGtFBJFOnWDsWA2BLiz02poKSa40xV9/De3b6yxwtLuZGVxPU3zoEKSnw4UX6hqDw6zbcZCZOQW8ua6AklIXKw5VhW+/1bafc44WDq/kIvqx/vZYXDOupe8tKfk+nUV+vi70JQtlZTpZ2btXAyLOPDPhJjh6n/7ud5pNdskSqFHDmc88BTZN8bGUlKjAt2kDDz7otTXxU7u2RpisW6duKIdp1aQOjWrXoKTUYZdFvJSV6YJmWZm6bqoQKXXs00UsrhnHn0iipKaqK/HQIV/7eF2hVy/YulWLhnsg8o4zdqxupBo40Hcpl5ND6EV04au0FMaPD8xutlOSkaEpUqdN04o9DuOayyIeevaEbdscFQfP2xkdtHNy4PbbE3tur3jgAb12ly7Vp7MwUK2a1sPds8d/xc9jCbZ3+3B9w9TTT4tcfLHI4cPunqccrqUprogJEzQR1BtvOP7RiUqqFRMjR2rOn0jiLSc5VTsT8j28+6627/HH3TuHH5g4Ua/XKVO8tsSd+7SwUKRFC5EXX3T+s4+BRCY1i/dwVehnzBBJSxPZudO9c1RAQoVeRHfNpqToLtowEhWHadMSfuq1n38lGQ8vlPNHz5eMhxe6K/YuDtq+4IMPdDB79FGvLRERF+/TLVtE6tUTWbTInc+PEKvQh9t1s3IlDB4Mc+fqqn6YeeYZjce+9lp9dAwT77+vRd0feQT69Uv46V0JrzwRgwfDsGFw223qzgkT27bpXodbb4XHHvPaGndp3lwL2Q8YoBW0PCa8Qr9tG9xyi6Y2uPxyr61JDAsW6C69yy7TyJQwsG0bdOniqThUJSQzLrKydNBu3z48g3ZRkYaSXnyx7gxOBtq21TXB7t2hoMBTU8IZy/XllxrT+uij+jNZqFZNo3DS0/Wm2rjRv/l7YqGoCK64QuOsPRSH6IavyoZkxsWCBRo6e+mlGnro9xQdJyMaRlmzJqxa5bU1ieXWW2H7dt0Yt2JFhbWBE0GAVeAE/Oc/cPPNGp1RiaLRbvH222+TkZFB8+bNnS+QXRFnnqlFw/PztcxZUImKQ61aruwTqCzNzjZ8MH40/TtfzQUXXMC0ee+5686pVk2jcEAHO5+F61WK667TydeGDZ7vE0j4/QgaYXT11Rre7cIGx1gIl9CXlcFdd0G9epCoTjwJpaWlDBkyhEWLFrFp0yays7PZtGmT+ydOS9Msie+8o7VUg4iPxAFg+PDhdOnShc2bN7N+/Xp+2r6N+yGZNWpoGOn27cEdtH/+c53Fr1rlWFGfqlLR/bhz5073T2yM5mICzeXjwSbVcAn9I4/oTTF5si9cFh9//DHNmzenadOmVK9enb59+/LWW28l5uRt2mgR6qwsLULtQyryca/bcZDcrj9DouLgg5TDRUVFLF++nLvvvhuA6tWrc/0l57mTv+dYgjxoP/UUvPYavPWWuqE8pqL7cfXq1Yk5eWoq/O1vsHq1FrNPMN5PlZzilVfUj7tyZcK2H5+K3bt307hctE9aWlriLizQR8W8PE341by55oXxCet2HKTfxFXfZbTMvucqAN7/+UhGvj+be/o+xuCzGtLKYzsB8vPzqVu3LnfddRfr16+nVatWPPPMM7RqUicxKSGig3a/frpJ7pe/dP+c8TJ7tqYZycryzTpZRffjhx9+mDgDatXSbJdt22rQRO/eCTu199NeJ1i8GH77W13AqlvXa2u+Qyp4RIu1pq5jPPKI5sru2tVX2fVm5RRwtKQMAY6WlDErp4AvXs9m5PuvMabjPXxwXit/pFwASkpKyMnJYfDgwXzyySfUrFkzcf7dKH36aNTRvfdqLhU/8+mnugj5q19pqKhP8MX9mJYG8+ZpxtWPPkrYaYMv9Bs3arrXGTNcy89eVdLS0ti1a9d3rwsKCmjYsGHiDZk6VUMuW7f2TZbEY2+5I6vW0GXMMLJb3cSU1j38k3IB7ce0tDR+8pOfAJCZmUlOdKE0kfh00P4v9u3TGet118GECV5b819UdD+e44Vr8NJL1b18yy2a6ycBBNt1c+iQxqhmZWnMsc9o06YNW7ZsYfv27TRq1Ijp06czbdo0b4z56CMtlFCnDtx/v+drGMOKvqXZhj0aTCLCHZ8uZFX6RZjxExhx5KjzWSLjoH79+jRu3Ji8vDwyMjJYsmQJLVu29MaYqVPVHffjH+vCnt/CLl97Ta+zxYu9tuQ4KrofBw0a5I0xXbrAmDFaD+Ozz8Dl6ynYQp+aqrG5Bw54bclxnHbaaSxatIg77riDa665hrKyMjp27Eh+fj75+fkJt8cUF9O2dm1qHzhA7sGDvkjs9uOmtdl8CAqPCDvPrkfeD9PJyVlPpzTYs0EPv5CZmUn37t0pLi6mfv36DB8+nHnz5iXeEBEu+8EPaAzkFhVRduRI4m04CfXr1qW4Vi3WzJ0LKSlem3Mcx96PzZo1886YaO76s892/VTBz0f/+ef6qPjii9Cjh6N2hQYRDXP78ktdJPPRDRjdeHTGv4qYMWUUpw8ZTPpj7tYKCDRjx8KcObBsmU5y/MbRo+pauugiTcthqZipUzXD5apV0KBBlT8mefLRn3eehm/dfbdWZbccz9ixupaRne0rkYfvd53e06MV386eS/rE57Q/LcczZQq8/LIu5vlR5AGqV9cSne+9Z4X+RCxbpu7TBQviEvlKEUvmM7cPR7JXzp4t0rChyOefx/9ZYeKvfxVp0kRkzx6vLYmNjz8WOfdc/Wn5nqVLRerWFdm40WtLYmP7dr0f58zx2hJ/sXmzyI9+JLJ4sSMfR9Jlr+zVC0aN0phdn0SWeM6yZbrJJpEzh3hp00Znrb16qVvOAps3a8z1tGnqEgkC0SftX/wifFk4q8r+/apP48ZBx44JPXV4hB5g+HD9An/2M/UVJjO5uSoO2dm+2JVYKXr21HKP3brZQTsqDk8+mXBxiJvWrWHSJDtoA3zzja4h3nabrpclmHAJPcCf/wxnnaW7B32w0OwJ+/erSP7xj77aDVsphg2zg3ZUHPr31xxOQaRHDx20k/lJu6xM6zo3awa//70nJrgi9MaYLsaYPGPMVmNMYkMoUlL0EXfjRl2ETDaOHNGb6/bb4c47vbYmPqKD9qBByTdol5VpH3ooDo4xbBh06qQbhJJx0H7wQZ18vfKKZ2HNjgu9MSYFGA90BVoC/Ywxid1dUrOm5pSYNEkjFZKF0lIVhxYtdDNG0IkO2p99Bo8/7rU1iWXUKN0f4qE4OMqf/qS52O+5J7kG7QkTNEpqzhw4/XTPzHBjRn8lsFVE8kXkKDAd6OnCeU5O/fq6CDlypC5KJgOjRsFXX+liZhjEAXTQnjcPXn01eQbtCRN0ojJ7tqfi4CgpKRo7npubPIP2ggXa1oULPc/C6obQNwJ2lXtdEHnvvzDGDDLGrDXGrC2M7hBzmgsv1Blh7956gYWZ8eP1wgqTOESJDtojRsDSpV5b4y4+EgfHKT9oT57stTXu8skn6jqdPRuaNvXaGleEvqKp5HHPaiLykoi0FpHWdd3MOHnDDboo2a2b+snCyPz58MQTKg51/JEfxnFattQ01H36hHfQzsnxlTi4Qr16Opg98EB4B+1du7RQzIsvwlVXeW0N4I7QFwCNy71OA7ytcHznnbrq3aOHLlaGiZwcjcgIszhE6dBBi1l06wZffOG1Nc6ya5denz4SB9do2VLDfsM4aBcV6fV5//0aMeYT3BD6NUALY8z5xpjqQF9grgvnqRyPPaYZ/26/XRctw8DOnSoOf/kLRFLohp6BA2HAgHAN2ocPa/jhiBG++I964wAABwBJREFUEgdX6dBBKy2FadAuLtY8/Nde67tqYI4LvYiUAEOBd4Bc4G8i8pnT56k0xsDEibpYOWqU19bEz+HDepOMGKFha8nEo49q7YEwDNpRcWjfXmeBycSAAXr89KfBH7RFtChMaqrm+PFZMETws1dWloMHNdvl0KGazzuIFBfrDDAjA557zncXVUI4ehQ6d4bLL9d4+yAiouGG+/Zp+J0PiqAnHBF9Svv6a3jzTd8l3YuZceO0+NHy5br3I0EkT/bKylKnji5aPvGELmIGDREYPFgLTmRlJafIg2ZJnDULFi2C55/32pqq8eSTusYyfXpyijzo9fvyy/qE+sADXltTNaZP17WV+fMTKvKVITmvrvPP1xlU9+4qFK38UII6RsaN09CtZcuSVxyiRAftdu2gSRN1AQSF7GwVh5UrfSsOCSOa2rhdO70377vPa4ti58MPdefvkiXgRZnQWIklxaXbhyNpiqvCrFkijRqJ7Njhzfkry7RpwUo5nChWr9YUvmvXem1JbCxfrvZu2OC1Jf5i+3aRBg1E5s712pLYyMsTqVdP5N13PTOBpEtTXBVuvll3znbrpo+OfmbFCs3OOX9+cFIOJ4orr4SXXtKslzt3em3NyfnnP3XxdepUuPhir63xF+WLCK1b57U1J6ewUHXjiSfgxhu9tuaUJLfQA/z613D99ZCZqYucfiQvT8UhSPnIE02vXvCb3+gitV8H7cJCtS8g4uAJbdrooN2jB+zY4bU1FfPNNzqp6N1bB6UAkOROXnQxKCtLZ/ddu2oki99YuBD+8Ifg5SNPNMOHQ36+fk9XXum1Ncfz979D376BEQfP6NVLRb5TJ39e8xs36tNHgHL2JF945Yn49791xuzHNKrp6cFaaPSS0lJNfvavf3ltyfHUrq2FJ5I1UqqyzJypoad+44wztEbAGWd4bUnM4ZVW6C0WiyWg2Dh6i8VisQBW6C0WiyX0WKG3WCyWkGOF3mKxWEKOFXqLxWIJOVboLRaLJeRYobdYLJaQY4XeYrFYQo4vNkwZYwoBnya2iIlzgQNeG+EQti3+IyztANsWp2kiInVP9Ue+EPqgY4xZG8vutCBg2+I/wtIOsG3xCuu6sVgslpBjhd5isVhCjhV6Z3jJawMcxLbFf4SlHWDb4gnWR2+xWCwhx87oLRaLJeRYobdYLJaQY4U+TowxXYwxecaYrcaY0V7bEyvGmMbGmA+MMbnGmM+MMcMj759jjFlsjNkS+VnHa1tjxRiTYoz5xBgzP/L6fGPM6khb3jDGVPfaxlgwxtQ2xrxpjNkc6Z+rg9ovxpj7I9fXP4wx2caYM4LSL8aYV4wx+40x/yj3XoX9YJRnIzqwwRhzhXeWH48V+jgwxqQA44GuQEugnzGmpbdWxUwJMFJELgCuAoZEbB8NLBGRFsCSyOugMBzILff6j8D/RdpyEAhKsdZngLdF5H+AS9E2Ba5fjDGNgGFAaxG5CEgB+hKcfnkN6HLMeyfqh65Ai8gxCHghQTbGhBX6+LgS2Coi+SJyFJgO9PTYppgQkb0ikhP5/WtUTBqh9r8e+bPXgV7eWFg5jDFpQDfg5chrA3QA3oz8SSDaYoz5AdAemAQgIkdF5BAB7RcgFahhjEkFzgT2EpB+EZHlwFfHvH2ifugJTBZlFVDbGNMgMZaeGiv08dEI2FXudUHkvUBhjDkPuBxYDdQTkb2ggwHwI+8sqxRZwCigLPL6h8AhESmJvA5K3zQFCoFXI26ol40xNQlgv4jIbuB/gZ2owB8G1hHMfolyon7wtRZYoY8PU8F7gYpXNcacBcwEfi0iRV7bUxWMMd2B/SKyrvzbFfxpEPomFbgCeEFELgf+TQDcNBUR8V/3BM4HGgI1URfHsQShX06Fr683K/TxUQA0Lvc6DdjjkS2VxhhzGiryU0VkVuTtL6KPnJGf+72yrxK0A3oYYz5H3Wcd0Bl+7YjLAILTNwVAgYisjrx+ExX+IPZLR2C7iBSKSDEwC2hLMPslyon6wddaYIU+PtYALSJRBNXRhaa5HtsUExEf9iQgV0T+XO6f5gIDI78PBN5KtG2VRUR+KyJpInIe2gfvi0h/4AMgM/JnQWnLPmCXMSYj8tYNwCYC2C+oy+YqY8yZkest2pbA9Us5TtQPc4EBkeibq4DDURePLxARe8RxADcB/wS2AQ95bU8l7L4GfbTcAHwaOW5CfdtLgC2Rn+d4bWsl23U9MD/ye1PgY2ArMAM43Wv7YmzDZcDaSN/MAeoEtV+AMcBm4B/AX4HTg9IvQDa6tlCMztjvPlE/oK6b8REd2IhGGnnehuhhUyBYLBZLyLGuG4vFYgk5VugtFosl5Fiht1gslpBjhd5isVhCjhV6i8ViCTlW6C0WiyXkWKG3WCyWkPP/6xRRws/+rFYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "q_h.plot()" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "129.38095238095238" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_h.chi2 #chi-squared test statistic for the observed point pattern" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "19" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_h.df #degree of freedom" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.909272893094198e-18" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_h.chi2_pvalue # analytical pvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similar to the inference of rectangle tessellation, since the analytical p-value is much smaller than 0.05, we reject the null of CSR. The point pattern is not random." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hexagon quadrats & empirical sampling distribution" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "q_h_e = qs.QStatistic(pp_juv,shape= \"hexagon\",lh = 15, realizations = csr_process)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.001" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q_h_e.chi2_r_pvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Because 0.001 is smaller than 0.05, we reject the null." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 1 }