{ "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": [ "