{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Heller, Heller, Gorfine (HHG)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial, we explore\n", "\n", "- The theory behind the HHG test statistic and p-value\n", "- The features of the implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| | $d_{\\mathbf{x}}(x_i,\\cdot) \\leq d_{\\mathbf{x}}(x_i,x_j)$ | $d_{\\mathbf{x}}(x_i,\\cdot) > d_{\\mathbf{x}}(x_i,x_j)$ | |\n", "| :-------------------------------------------------------: | :---------------------------------------------------------------: | :---------------------------------------------------: | :----------------: |\n", "| $d_{\\mathbf{y}}(y_i,\\cdot) \\leq d_{\\mathbf{y}}(y_i, y_j)$ | $A_{11}(i, j)$ | $A_{12}(i, j)$ | $A_{1\\cdot}(i, j)$ |\n", "| $d_{\\mathbf{y}}(y_i,\\cdot) > d_{\\mathbf{y}}(y_i, y_j)$ | $A_{21}(i, j)$ | $A_{22}(i, j)$ | $A_{2\\cdot}(i, j)$ |\n", "| | $A_{\\cdot 1}(i, j)$ | $A_{\\cdot 2}(i, j)$ | $n-2$ |\n", "\n", "The following description is adapted from [[1]](https://arxiv.org/abs/1907.02088):\n", "\n", "HHG is a consistent multivariate test of associations based on the rank of the distances [[2]](https://academic.oup.com/biomet/article-abstract/100/2/503/202568). For every sample point $j \\neq i$, denote a point in the joint sample space as $(x_j, y_j)$. Let $d_{\\mathbf{x}} (x_i, x_j)$ be equivalent to the norm distance between samples $x_i$ and $x_j$ and $d_{\\mathbf{y}} (y_i, y_j)$ is similarly defined. The indicator function is denoted by $\\mathbb{I} \\{ \\cdot \\}$. The cross-classification between these two random variables can be formulated as in the above table, where\n", "\n", "$$A_{11} = \\sum_{k = 1, k \\neq i, j}^n \\mathbb{I} \\left\\{ d_{\\mathbf{x}} (x_i, x_k) \\leq d_{\\mathbf{x}} (x_i, x_j) \\right\\} \\mathbb{I} \\left\\{ d_{\\mathbf{y}} (y_i, y_k) \\leq d_{\\mathbf{y}} (y_i, y_j) \\right\\},$$\n", "\n", "and $A_{12}$, $A_{21}$, and $A_{22}$ are defined similarly. $A_{\\cdot 1}$, $A_{\\cdot 2}$, $A_{1 \\cdot}$, and $A_{2 \\cdot}$ are the sums of the columns and rows respectively.\n", "\n", "Once this table is generated, the Pearson's chi square test statistic can be calculated using\n", "$$S (i, j) = \\frac{(n - 2) {\\left( A_{12} A_{21} - A_{11} A_{22} \\right)}^2}{A_{1 \\cdot} A_{2 \\cdot} A_{\\cdot 1} A_{\\cdot 2}}.$$\n", "\n", "From here, the HHG test statistic is simply\n", "\n", "$$\\mathrm{HHG}_n = \\sum_{i = 1}^n \\sum_{j = 1, j \\neq i}^n S (i, j).$$\n", "\n", "The p-value is then calculated using a standard permutation test." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using HHG" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before delving straight into function calls, let's first import some useful functions, to ensure consistency in these examples, we set the seed:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt; plt.style.use('classic')\n", "import seaborn as sns; sns.set(style=\"white\")\n", "\n", "from mgcpy.independence_tests.hhg import HHG\n", "from mgcpy.benchmarks import simulations as sims\n", "\n", "np.random.seed(12345678)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, let's simulate some linear data:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAIlCAYAAABrdaqpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dfZxcZX338c8GSFKWJLjeGKCg+FB/NxIUEKyioFHKLdjerdaHqhVtS1Gr+EC0YgVFxWdAMEoVRRChiuhtrbY+IlEReWhVRIRfKwYFhKCsJmElWcju/cc5K8Owu5m9srMzs/N5v155ze6Z65z9zXDY+e51Xec6A+Pj40iSJM3Ugk4XIEmSepMhQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkCYiIgU7XIPWa7TtdgDTfRcQaYHFmPn6aNicBbwH+IDM3zVFpxSLiAcAbgL8A9gI2A9cB5wMfzsx76nZPAS4BjsjMr8xBXecCT8/MXWewzwBwIrAJeG+97SR66L+H1Cn2REjd4WPAE6g+jLtaRPwB8B3gucAHgCOBFwCXAe+nChITvk/1ur43x2XOxCLgrcCODdt65r+H1En2REhdIDNvBm7udB0tejawD/CYzPxRw/Z/j4jfAm+NiHdn5g8zcwNweUeq3AY99t9D6hhDhNQFmrvP6275vYCzgBOAhwM3Au/MzE807Lewfv5FwO7Az4EzM/P0puO/BHg51Yf/DsDP6nar6+f3AtYCxwF/C+wJnJyZp0xS7sRQwWQ9mWcDo8DG+rhPoWE4o67jHODxVL0W+wPrgJOArwCrgSOADcA5mfmmpvpenpkfnup9ay4mIhYArwWOAv4IGACuB96dmRc2HBfgLRHxlswcmOy4EbE/8HbgcVS9Ft8DTszMy5tqfB7w58AzgO2ALwOvyszbJnm/pJ7mcIbUvfaj6mZ/F/CnwC+AcyNi34Y2FwGrgI/UbT4LnBoR75poEBEvBT4OfB34M6qehBuBD9Qf8o3eXR/rKOBLU9T1FeAe4KsR8Y6IOLQe4iAzb8nMd2fmDVt5bZ8DPlnXcyPV8MEa4L+BZwEXA/8UEc/eynG25h31v/Oohl3+mirkXBARDwduBZ5ctz2bagjjfiLiycAVwFLgZVTvz47AtyLikKbmHwF+Q/U+H1+/xg9t4+uQupI9EVL3WgYckpnXAEREUvU0/F/gmoh4av3132bmOfU+X4+ITcCbI+JDdbf8I4AzMvOEiQNHxHeBO4CVVB/eE76YmR+crqjMvCYi/hL4MPBP9b+7I+IK4ELgrMwc3cprOzUz/7mu5e66hqsaeh4uAZ4DPJEqGJXaE3hLZp46sSEi1gL/BRyamedExJX1UzdP9CpM4j3ATcBhE68tIv6dqlfjFOCPG9p+MzNfWX/9jYh4LPDXETGQmd6sSPOKIULqXhsnAkRtYox+sH78k/rx3yKi8f/lL1D1YDwN+ERmvh4gIpYCj6QKFY+t2y5s+pnX0ILM/LeI+DJwaP1zDqH6IH0S8NKIeEpm3jHNIb7b8PW6+vH3H+CZeXdErAce0Eo909T51wAR8UDufe0r66ebX/ukImKQagjjvY3hKDM3R8RngNdHxE4Nu3y36RA3Uw0h7UDVCyLNG4YIqXv9rvGbzByLCLh3GPJ/1Y+/nmL/PwSIiIcC/wwcDmyhGjKY+KBrXhthHS3KzLuphh0urn/OUuCNVF34bwD+cZrdN0yybaTp+23+qz0iDgA+yL1XWlzHvUGp1XUhdq7bTjan4db6uaUN237X1GasfnT4WPOOJ7XUu35L9Zft44CDJvl3br0GwpeAh1D1Egxm5j7Aa0p/aERcFhEXNm/PzA2Z+UaqkLJP6fGnMBEotmvavmSqHSJiCdX8jTGq+SWDmbk/1byPmfht/fMnW3ti9/q56XpdpHnLECH1rjVUXfI7ZuZ/Tvyj+sv5ZGA3YBfgUcC5mXlZQ3f8kfVjye+AG4A/j4hHNT8RETsDDwJ+dL+9ts1Ez8WeTdubJzU22pvq9a/OzKszc0u9vfm1b7nfng0ycwS4EnhOROwwsT0iFlGtlXFFZrqehPqSwxnS3Ng1Iib76//WzLzfX/Ut+jLV5ZMXRcQ7gKupPjhPpupm/3E9br8WeFn9+CuqD943UP0FPTjpkaf3T8BTgMsi4kNUC0/9rv7ZrwGGgVOn3LtAZv4mIi6leh3XAbdQXYr60Gl2ux5YDxwfEXfVNR4JTEx6HKyPfXdE3AkcHBGH1q+n2RuBrwEXR8T7623HUQ0ZvWRbXpvUy+yJkObGQ6jWRWj+99rSA2bmGNVlnefUx/ka1Qf8RcDKhr+O/y/VuhAfo7q08hnA0VQh5NCCn3sT1foOH6VaD+EiqstHXwf8B/C4zJxqnsa2eDHVXI4zgU9ThZU3TFPnBqrXvhn4VP3vAKr37Cfce2knVOs/HET1njT3dpCZlwBPpbq09XzgXKo5HIdm5pptelVSDxsYH/eKI0mSNHP2REiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVGT7ThewNRHxEWC7zDx6mjYXAc9u2nxxZh7W1uIkSepjXRsiImIAeCtwDHD2VpqvAI4HPtGwbXObSpMkSXRpiIiIh1EFhxXAL7bSdiHwCODKzLxtDsqTJEl075yIJwA/A/YF1m6l7d5UYei6dhclSZLu1ZU9EZl5AXABQERsrfkKYBR4a0QcAdwFXAScnJmb2lmnJEn9rCtDxAztAwwACXyQqvfiNGBP4MWtHCAiFgEHAbcCW9pTpiRJXWM7YDfgqswsnkM4H0LECcApmTlcf39NRGwBPh0Rx2XmHS0c4yDgO22rUJKk7nQIcGnpzj0fIjJzDBhu2nxN/bgn0EqIuBXgggsuYNddd53F6iRJ6j633XYbL3zhC6H+/CvV8yEiIj4D7JCZz2zYfCDVJZ4/bfEwWwB23XVX9thjj1muUJKkrrVNQ/g9FyLqSzqHgOHMHAU+Sz10AXwB2B84hWqI487OVSpJ0vzWrZd4Tudgqu6XgwEy8zPAS4C/AX4MnAqcAby5Q/VJktQXur4nIjOf0vT9GqqrMRq3nQecN3dVSZKkXuyJkCRJXcAQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBUxREiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSpJ6yZWyc4fWbuG7tHQyv38TY2HinS+pb23e6AEmSZmL9xs0ce+olbBgZZengQlavWsnQssWdLqsv2RMhSeop64ZH2DAyCsCGkVHWDY90uKL+ZYiQJPWU5UODLB1cCMDSwYUsHxrscEX9y+EMSVJP2XnJIlavWsm64RGWDw2y85JFnS6pbxkiJEk9ZcGCAYaWLXYeRBdwOEOSJBUxREiSpCKGCEmSVMQQIUmSinT9xMqI+AiwXWYePU2bA4EzgP2BW4C3Z+Z5c1SiJEl9qWt7IiJiICLeBhyzlXa7AF8Fvg8cAHwAODsiDm9/lZIk9a+u7ImIiIcBZwMrgF9spfnRwHrg1Zk5BlwfEQcArwO+1tZCJUnqY93aE/EE4GfAvsDarbQ9BPh2HSAmrAGeGBHd+vokSep5XdkTkZkXABcARMTWmu8B/KBp2y+BHYEh4NezXZ8kSerenoiZ2BHY1LRtc/3ocmaSJLXJfAgRdwHNC6dPfO+t3SRJapP5ECJuAnZr2rY7cCfVhEtJktQG8yFEXAocGhEDDdtWAt9tmmwpSZJmUVdOrJxORCykmjA5nJmjVJeC/iPw4Yg4HTgMeAHw9M5VKUnS/NeLPREHA7fWj2TmOqrAsD/VVRqvBI7KzG92rEJJkvpA1/dEZOZTmr5fAww0bbsceNzcVSVJknqxJ0KSJHUBQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSQBsGRtneP0mrlt7B8PrNzE2Nt7pktTltu90AZKk7rB+42aOPfUSNoyMsnRwIatXrWRo2eJOl6UuZk+EJAmAdcMjbBgZBWDDyCjrhkc6XJG6nSFCkgTA8qFBlg4uBGDp4EKWDw12uCJ1O4czJEkA7LxkEatXrWTd8AjLhwbZecmiTpekLmeIkKQet2VsnPUbN9/nw3/BgoEZH2fBggGGli12HoRaZoiQpB7nhEh1inMiJKnHOSFSnWKIkKQe54RIdYrDGZLU45wQqU4xREhSj3NCpDrF4QxJklTEECFJkooYIiRJUhFDhCRJKmKIkCRJRQwRkiSpiCFCkiQVMURIkqQihghJklTEECFJkooYIiRJUpGuvHdGRGwHnAy8BFgCfAV4RWaum6L9RcCzmzZfnJmHtbNOSZprW8bGWb9x831utrVgwUCny1Kf6soQAZwEvBg4CrgDOBP4HPCkKdqvAI4HPtGwbXMb65Okjli/cTPHnnoJG0ZGWTq4kNWrVnrjLXVM14WIiFgIvBp4VWZ+vd72V8DaiDg4My+bpP0jgCsz87Y5L1iS5tC64RE2jIwCsGFklHXDI4YIdUw3zonYj2oIY83Ehsy8EbgROGSS9ntThaHr2l+aJM2uLWPjDK/fxHVr72B4/SbGxsanbb98aJClgwsBWDq4kOVDg3NRpjSpruuJAPaoH29p2v5LYM9J2q8ARoG3RsQRwF3ARcDJmbmpbVVK0iyY6fDEzksWsXrVyvvMiZA6pRtDxI7AWGbe3bR9MzDZ/1n7AANAAh8E9gVOowocL25jnZK0zWY6PLFgwQBDyxY7hKGu0I3DGXcBCyKiOeAsAkYmaX8CsGtmnpaZ12Tmv1DNqTgqIh7Y5lolaZs4PKFe1o09ETfVj7s1fA2wO/cf4iAzx4Dhps3X1I97Ul3dIUldyeEJ9bJuDBFXAxuBJwPnA0TEXsBewLebG0fEZ4AdMvOZDZsPpBr++Gmba5WkbeLwhHpZ14WIzNwcEWcCp0TEr4HbqdaJ+FZmXl5f0jkEDGfmKPBZ4NMRcRzwBWB/4BTglMy8szOvQlK/cjEo9ZOuCxG1E4AdqHoidqBesbJ+7mDgEmAlsCYzPxMRi4HXA++gCh1nAO+a66IlycWg1E+6MkRk5j3Aqvpf83NrqK7GaNx2HnDenBQnSdNwMSj1k268OkOSepZXW6ifdGVPhCT1Kq+2UD8xREjSLPJqC/UThzMkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISX1jy9g4w+s3cd3aOxhev4mxsfFOlyT1tO07XYAkzZX1Gzdz7KmXsGFklKWDC1m9aiVDyxZ3uiypZ9kTIalvrBseYcPIKAAbRkZZNzzS4Yqk3maIkNQ3lg8NsnRwIQBLBxeyfGiwwxVJvc3hDEl9Y+cli1i9aiXrhkdYPjTIzksWdbokqadN2RMREa+Yy0Ikqd0WLBhgaNli9n7oAxlatpgFCwY6XZLU06YbzjgjIr4REXvOWTWSJKlnTBcingQsB34cEX83R/VIkqQeMWWIyMzLgf2B9wEfjIj/iIg/nLPKJElSV5t2YmVm3gOcHBGfAs4AromIM4HfNbV7Z/tKlCRJ3ajVqzNuBq4CDgf+Ftjc8Nw4YIiQJKnPbDVERMSRwIeABwKvzcwPtb0qSZLU9aYMERHxIOADwHOAbwDHZObP56owSZLU3abribi+fjw6M8+Zi2IkSVLvmC5EfAt4eWbeNlfFSJKk3jFliMjMZ85lIZIkqbd4Ay5JklTEECFJkooYIiRJUhFDhCRJKjLdOhHfbPUgmfnU2SlHkiT1iul6Im5o+Hcr8BRgELgG+C9gADgE+El7S5QkSd1ouks8/37i64j4JPCezHxjY5uIOAF4XPvKkyRJ3arVORHPBD4+yfYLgcNmrxxJktQrWg0Rt1MNXTQ7Arhp9sqRJEm9otVbgZ8KnBkRjwe+TzUf4mDgecBL2lOaJEnqZi2FiMz8UERsAF4BvAAYB34IPCczv9DG+iRJUpdqtSeCzPwk8Mk21iJJknpIyyEiIh4BrAL+N/DXwJ8DP8nMNe0pTZIkdbOWJlZGxB8DVwMPp5oLsQjYB/h6RPxp+8qTJEndqtWrM94NvDszDwdGATLzFcB7gZPaU5okSepmrYaIA4BPTbL9bGDv2StHkiT1ilZDxAjwoEm2PxJYP3vlSJKkXtFqiPgU8P6I2Jvq8s7FEfFU4IPAZ9tVnCRJ6l6thog3AjcC1wI7Ud2E6+vAVfVzkiSpz7S62NQo8LyIeBOwH9Xkymsz84Z2FidJkrpXSyEiIn4GHJiZPwV+2rB9N+DqzJxsvoQkSZrHpgwREXEkcGD97V7A8RFxZ1OzR053jFIRsR1wMtV9OZYAXwFekZnrpmh/IHAGsD9wC/D2zDxvtuuSJEn3mi4ArAVOp7rZFsCzgS0Nz48DG4Fj21DXScCLgaOAO4Azgc8BT2puGBG7AF8F/gX4O+BPgLMj4rbM/FobapMkSUwTIjLzOqqeBiLiEuBZmfmbdhcUEQuBVwOvysyv19v+ClgbEQdn5mVNuxxNdZnpqzNzDLg+Ig4AXgcYIiRJapOWrs7IzJWTBYiIWBgRT5zlmvajGsJY0/Dzb6S6OuSQSdofAny7DhAT1gBPjIhWrz6RJEkz1OrEyscCHwX2ZfLgsd0s1rRH/XhL0/ZfAntO0f4Hk7TdERgCfj2LtUmSpFqrf6mfAdwFHEN1eefLqe6bsRl43izXtCMwlpl3N23fDCyeov2mSdoyRXtJkjQLWg0R+1PNUTgH+CGQmflG4A1UgWI23QUsiIjmXpJFVMtvT9Z+0SRtmaK9JEmaBa2GiAHgV/XX/0M1rAHwReAxs1zTTfXjbk3bd+f+QxwT7Sdreyfe10OSpLZpNUT8GDiy/vonwMRkyuXM7nwIgKupLh198sSGiNiLaq2Kb0/S/lLg0IgYaNi2Evhu02RLSZI0i1pdKOo9wIURsYXqZlxviYh/peqFuGQ2C8rMzRFxJnBKRPwauJ1qnYhvZebl9SWgQ8BwvRz32cA/Ah+OiNOBw4AXAE+fzbokSdJ9tXqJ5+eAxwNXZubPqXolRoH/AP6+DXWdAFwAnE8VUn5OtdgVwMHArfUj9SqWT6eat/ED4JXAUZn5zTbUJUmSagPj4+OdrqHj6uGStRdffDF77LHH1ppLktTTbr75Zp72tKcBPLRei6nIdPfOaHm1x8w8vLQASZLUm6abEzHZlRCSJEnA9PfO+Ju5LESSJPWWVpe9fsF0z2fmv8xOOZL60ZaxcdZv3My64RGWDw2y85JFLFgwsPUdJXVUq5d4nj/F9k3AzVS34ZakIus3bubYUy9hw8goSwcXsnrVSoaWuWq91O1aChGZeZ9LQSNiO6rbhP8z8JE21CWpj6wbHmHDyCgAG0ZGWTc8YoiQekDRrbIzc0tmXgccB7x9dkuS1G+WDw2ydHAhAEsHF7J8aLDDFUlqRavDGVO5h+o+FZJUbOcli1i9auV95kRI6n7bMrFyKdWtwa+Y1Yok9Z0FCwYYWrbYIQypx2zLxMq7ge8B/zB75UiSpF5RNLFSkiTJcCBJkoq0OifiscAHgRXA/WY8ZebCWa5LkiR1uVbnRHyM6tbfrwfual85ktrFVSElzbZWQ0QAB2Xmte0sRlL7uCqkpNnW6pyI7wMPbmchktprslUhJWlbtNoTcQzw+Yg4CPgZMNb4pDfgkrrfxKqQG0ZG2e+R/4tdHrAj1629w6ENScVaDRHPAv4IOGmS58bxBlxS12tcFXKXB+zIq09b49CGpG3Saoh4DXACcHpm/q6N9Uhqk8ZVIa9be4c3vJK0zVqdE7Ed8CkDhDQ/eMMrSbNhJpd4vgx4QxtrkTRHvOGVpNnQaohYBhwVEc8HbqC6b8bvZebhs12YpPbxhleSZkOrIWIH4FPtLESSJPWWVm/A9TftLkSSJPWWKUNERLwA+GxmjtZfT2U8M+2lkCSpz0zXE3E+8A3g9vrrqYzjUIckSX1nyhCRmQsm+1qSJAlan1j5exGxPfBoYF1m3jL7JUmSpF4wbQ9DRLwoIv4zIh5cf/8o4KfAVcDPI+JjEbHdHNQpSZK6zJQhIiKeC5wL/BiYuN3fJ4GlwNOBg4HHUy2JLUmS+sx0wxmvAk7IzHcBRMRjgP2Bt2fm1+ttJwJvB05td6GSJKm7TDec8Wjg8w3fH0Z1JcYXG7b9CHh4G+qSJEldbroQsQAYbfj+UGA98F8N2/4A2NSGuiQV2jI2zvD6TVy39g6G129ibGy80yVJmqemG864Fngi8LOIWAo8DfhSZjb+RvpLqjkTkrrE+o2bOfbUS9gwMsrSwYWsXrXSe2RIaovpQsSHgNUR8WiqMPEHwOkAEfEg4AXA8cBL212k1E+2jI2zfuPm+9xhc8GCgZb3Xzc8woaRqhNxw8go64ZHDBGS2mK6xabOi4jFwDHAFuB5mXl5/fRbgKOB92bmee0vU+of29qTsHxokKWDC3+///KhwTZWK6mfTbvYVGaeBZw1yVPvAt6cmXe0pSqpj21rT8LOSxaxetXK+/RkSFI7zHjFSoDMvHm2C5FU2daehAULBhhattghDEltVxQiJLWPPQmSeoUhQuoyc9mTsK2TOCX1N0OE1Me8HFTStvAW31Ifm2wSpyS1yhAh9bGJSZyAl4NKmjGHM6Q+5iROSdvCECH1MS8HlbQtHM6QJElFDBGSJKmIIUKaAW+zLUn3ck6ENAOuqyBJ97InQpoB11WQpHsZIqQZcF0FSbqXwxnSDLiugiTdyxAhzYDrKkjSvRzOkCRJRQwRkiSpiCFCkiQVMURIkqQihghJklTEECHNEpfEltRvvMRTmiUb7tzMq067hPV3uiS2pP5gT4S0jSZ6IG6+fSNv+bvH88gHP8AlsSX1BXsipG3UfFOuk45+PG89+3KXxJY07xkipG3UfFOuTaP3sHrVSpbt5JLYkua3rgsREfEg4IPA4cAocA7wpsy8Z5p9bgd2adp8Ymae3LZCpdrETbkmeiL+cJclPGCpcyEkzX9dFyKAzwHjwJOBPwTOBe4B3jRZ44hYThUgDgX+p+GpjW2tUqp5Uy5J/aqrQkREPAF4EvCwzFwLXB0RrwdWR8TbMnPzJLutoAoZV2Tm6ByWKwHelEtS/+q2qzMOAX5eB4gJa4AlwH5T7LMCuMEAIUnS3OqqnghgD+CWpm2/rB/3BK6YZJ8VwD0R8SXgwHr/0zPzk22rUpIkzW2IiIi9gLVTPL0ZOB/Y1LgxM++OiHFgqr7ifYAHAidSzZs4AjgnIrbPzHNmo25JknR/c90TcQuw9xTPjQHHAveZlRYROwADwFQr96wEFmbmxETKqyPiIcBxVFd2SJKkNpjTEJGZdwPXT/V8RNwEHNm0eff6sXmYY+KYm6l6MRpdAzy/sExJktSCbptYeSnwsIjYs2HbSqrLNX/Y3Dgito+ImyLitU1PHQhc274yJUlSt02s/B5wOXBhRLwSWA68Bzht4uqLiNgJ2Ckzb8vMeyLii8AJEXED8BPgL4AXAc/oyCuQJKlPdFVPRGaOA88E1gHfoZrTcDbwtoZmrwNubfj+tcCHgQ9Q9T68CHhuZn5tLmqWJKlfDYyPj3e6ho6buGrk4osvZo899uh0OZIktdXNN9/M0572NICHZuaNpcfpqp4ISZLUOwwRkiSpiCFCkiQVMURIkqQihghJklTEECFJkooYIiRJUhFDhCRJKmKIkCRJRQwRkiSpiCFCkiQVMURIkqQihghJklTEECFJkooYIiRJUhFDhCRJKmKIkCRJRQwRkiSpiCFCkiQVMURIkqQihghJklTEECFJkooYIiRJUhFDhCRJKmKIkCRJRQwRkiSpiCFCarBlbJzh9Zu4bu0dDK/fxNjYeKdLkqSutX2nC5C6yfqNmzn21EvYMDLK0sGFrF61kqFliztdliR1JXsipAbrhkfYMDIKwIaRUdYNj3S4IknqXoYIqcHyoUGWDi4EYOngQpYPDXa4IknqXg5nSA12XrKI1atWsm54hOVDg+y8ZFGnS5KkrmWIkBosWDDA0LLFzoOQpBY4nKE55xUQkjQ/2BOhOecVEJI0P9gToTnnFRCSND8YIjTnvAJCkuYHhzM057wCQpLmB0OE5pxXQEjS/OBwhiRJKmKIkCRJRQwRkiSpiCFCkiQVMURIkqQihghJklTEECFJkooYIiRJUhFDhCRJKmKIkCRJRQwRkiSpiCFCHbNlbJzh9Zu4bu0dDK/fxNjYeKdLkiTNgDfgUses37iZY0+9hA0joywdXMjqVSu9KZck9RB7ItQx64ZH2DAyCsCGkVHWDY90uCJJ0kwYIjQrSoYmlg8NsnRwIdr/za0AAA+OSURBVABLBxeyfGiw3WVKkmaRwxmaFSVDEzsvWcTqVStZNzzC8qFBdl6yaI6qlSTNBkOEZsVkQxNbCxELFgwwtGyx8yAkqUc5nKFZ4dCEJPUfeyI0KxyakKT+Y4jQrHBoQpL6j8MZkiSpiCFCkiQVMURIkqQihghJklSkaydWRsQi4ErgfZl5/lbavhB4M/Bg4Grg2My8qv1VSpLUv7qyJyIilgCfBx7dQtvDgI8DpwIHANcAX4uIXdpapCRJfa7rQkQdCn4ILG9xl9cDn8rMszLzOuClwDDw920qUZIk0YUhAjiSqmfh4K01jIgFwBOBNRPbMnMM+DZwSJvqkyRJdOGciMw8buLriNha852BQeCWpu2/BA6a3cokSVKjOQ0REbEXsHaKpzdn5kyXO9yxftzUfCzApRMlSWqjue6JuAXYe4rnxgqOd1f92HyjhkXASMHxJElSi+Y0RGTm3cD1s3jIYaqwsFvT9t25/xCHJEmaRd04sbJlmTkOXAY8eWJbPdnyUKrJlZIkqU26bmLl1kTETsBOmXlbvek04IsR8QPgm8BxwDLgYx0qUZKkvtCLPRGvA26d+CYzvwIcA6wCvg88Cjg8M3/dmfIkSeoPXd0TkZkDk2w7CTipads5wDlzU5UkSYLe7ImQJEldwBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIoYISZJUxBAhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFdm+0wVMJSIWAVcC78vM87fS9nZgl6bNJ2bmye2qT5KkfteVISIilgAXAo9uoe1yqgBxKPA/DU9tbE91kiQJujBERMRhwEeA37a4ywrgHuCKzBxtW2GSJOk+unFOxJHAx4GDW2y/ArjBACFJ0tzqup6IzDxu4uuIaGWXFcA9EfEl4EDgFuD0zPxkeyqUJEkwxyEiIvYC1k7x9ObMXFxw2H2ABwInAm8CjgDOiYjtM/OcFo+xHcBtt91W8OMlSeotDZ93223Lcea6J+IWYO8pnhsrPOZKYGFmTkykvDoiHgIcB7QaInYDeOELX1hYgiRJPWk34IbSnec0RGTm3cD1s3zMzcDmps3XAM+fwWGuAg4BbgW2zFJpkiR1q+2oAsRV23KQrpsTMRMRsT3V8Mhpmfn+hqcOBK5t9Th1ELl0lsuTJKmbFfdATOi5EBEROwE7ZeZtmXlPRHwROCEibgB+AvwF8CLgGZ2sU5Kk+a7nQgTwOuAtwED9/WuB3wAfoOqauR54bmZ+rTPlSZLUHwbGx8c7XYMkSepB3bjYlCRJ6gGGCEmSVMQQIUmSihgiJElSEUOEJEkq0ouXeM6aiFgEXAm8LzPP30rb24FdmjafmJknt6u+bjHD9+mFwJuBBwNXA8dm5jatiNYLIuJBwAeBw4FRqiXX35SZ90yzz7w/pyJiO+Bk4CXAEuArwCsyc90U7Q8EzgD2p1om/+2Zed7cVNtZBe/VRcCzmzZfnJmHtbPObhIRHwG2y8yjp2nTt+fUhBbfp6LzqW97IiJiCfB54NEttF1O9cv+UKq1KCb+vX+6/eaDGb5Ph1Hdxv1U4ACq5ce/FhHNH5Tz0eeAXYEnU30I/A3w1qka99E5dRLwYuAoqte6B9V7dT/1efJV4PtU588HgLMj4vA5qbTzTqLF96q2Ajie+54/z2lvid0hIgYi4m3AMVtp19fnVKvvU63ofOrLnoj6w+4jwG9b3GUFcA9wRWaOtq2wLlPwPr0e+FRmnlXv/1LgqcDfA+9sS5FdICKeADwJeFhmrqW6CdzrgdUR8bZ6WfVm8/6cioiFwKuBV2Xm1+ttfwWsjYiDM/Oypl2OBtYDr87MMeD6iDiAaoG5eb143Ezfq7r9I4ArM7Ovbj8cEQ8Dzqb6f+gXW2nez+dUy+/TtpxP/doTcSTVX8wHt9h+BXDDfP1lP42W36eIWAA8EVgzsa3+n/bbVDc3m88OAX5eB4gJa6i6pPebYp9+OKf2o3oP1kxsyMwbgRuZ/Jw4BPh2fd5MWAM8sT6/5rOZvld7U/0ReF37S+s6TwB+BuxLde+k6fTzOTWT96n4fOrLnojMPG7i64hoZZcVwD0R8SWqm3vdApyemZ9sT4XdYYbv087AINV70+iXwEGzW1nX2YPJXzfAnsAVk+zTD+fUHvXjZO/NnlO0/8EkbXcEhoBfz2p13WWm79UKqrk3b42II4C7gIuAkzNzU9uq7AKZeQFwAbT0e6lvz6kZvk/F59O8CxERsRdTp67Nmbm44LD7AA8ETgTeBBwBnBMR22fmOUWFdlgb3qcd68fmE24zUPKed42tvVfA+TS97sy8OyLGmfq1z7tzahI7AmOZeXfT9qnOiR2Z/PxhivbzyUzfq32o7h+UVBN69wVOowocL25jnb2mn8+pmSg+n+ZdiKBK8ntP8dzYFNu3ZiWwMDM31t9fHREPAY6jmoXfi2b7fbqrflzUtH0RMFJwvG6ytffqWJped0TsQPU/5VSvfT6eU83uAhbUwajxKpWpzom7mPz8YYr288lM36sTgFMyc7j+/pqI2AJ8OiKOy8w72lxvr+jnc2omis+neRci6iR//SwfczP3ptcJ1wDPn82fM5fa8D4NU/1PuVvT9t25fxdtT9naexURN1HNH2m0e/046Wufj+fUJG6qH3dr+BqmPiduYvLz506qyXHz2Yzeq3qMf7hp8zX1456AIaLSz+dUy7blfJrvE0u2WURsHxE3RcRrm546ELi2EzV1o8wcBy6jusQR+P1ky0OpJlfOZ5cCD4uIxrHrlcBG4IfNjfvonLqa6j1oPCf2AvZi8nPiUuDQiBho2LYS+G7TxLj5aEbvVUR8JiI+37T5QKpg+tO2Vdl7+vmcatm2nE/zridiNkTETsBOmXlbZt4TEV8EToiIG4CfAH8BvAh4Rifr7LTG96nedBrwxYj4AfBNqq75ZcDHOlTiXPkecDlwYUS8ElgOvAc4beLqi348pzJzc0ScCZwSEb8GbgfOBL6VmZfXl5UNAcP1+3Q28I/AhyPidOAw4AXA0zvzCuZOwXv1WequZuALVAspnULVJX1nZ15F53lOtWY2zyd7Iib3OuDWhu9fC3yYaqGSa6l+2T83M+f1dcYtuM/7lJlfoVrUZBXV4i6PAg7PzHk7Axp+3wvzTGAd8B2qOQ1nA29raNav59QJVDPEzwcuAX7OvaviHUz1nhwMUK/M+HSqX2A/AF4JHJWZ35zjmjtlJu/VZ7h3UbMfUy3wdgbVarH9zHOqNbN2Pg2Mj4+3rUpJkjR/2RMhSZKKGCIkSVIRQ4QkSSpiiJAkSUUMEZIkqYghQpIkFTFESJKkIq5YKfWJiDiKarGdfahuHPYj4AOZeWFDm3HgRZl5fptqOBfYIzMPa7H9o4CHZua/b8PP/BjwiMx8yiTPPZxqyenzM/NlTc/9KdXqff8nM79R+vOl+cyeCKkPRMQxVLf4PRN4DPDHwL8Dn4qIxlv97ka1BG63+AJwULsOnpk3UC3P/tKI+P1N1Or7oHwCeK8BQpqaPRFSf3gp8NHMPLdh208iIoBXU31g0nAflG4xsPUm2yYzz4qIPwPOjoh9gd8Cnwb+Bzix3T9f6mWGCKk/bAGeGBHLMrPxFsivAwYnvmkczqiHHu6mumXy0fUxTgf+H3AWcACQwNGZ+Z/N+092zOaiIuIvgeOBFcA41f0NXpOZV0XEGuDhwFsi4iWZuVdELALeSXUTpcG6/Rsy8/KGY76ifl3L61pb6XE9murWx6uB6+p69svMe1rYV+pbDmdI/eF9wOOAX0bEv0XE6yJiv8z8VWbeOM1+R9WPjwXeT3VTsX8F3lUfbxT4UElBEXEQ8BngXGBvqttgDwAfrZs8C7iR6mZAE0Ma51HdXv65VLcq/iZwSUQ8sj7mi6juJvtOqpsu3Qw8f2u11DdqOgb4K6reh5dm5tqS1yX1E0OE1Acy8yLgSVTzIA6lChU/iIjvR8Q+0+z6K+D19dyB99fb/iUzv5SZ11DdsXRFYVl3A/+QmR/KzBsz8yqqALFvXfMwVe/HnZn5q4h4BFV4eElmficz/zsz3wpcSnXnWKgmjp6fmR/NyvHAVS3W8y3gN8A9wGWFr0nqKw5nSH0iMy8DLouI7ah6Fv4MOBb4ckQ8IjNHJ9nthvpW52TmSDWFghsanr8LWFRYzw8j4rcR8Uaq28b/EbAfU/9xs3/9eEVdx4RFDTWsoJ7f0eBy4NEtlPRR4A6q34vnRcRTM3Oshf2kvmWIkOa5+kqDNwJvy8zbMnMLcCVwZUR8B/gq1Yfsf06y+92TbGv5gzUipvwdExErgS9TXYHxXeDjwCOBD0+xy0TIeQJVeGm0uX4c5/6TMScLR821vBR4JvBUqt+LX6fq3Xjf1vaV+pnDGdL8dxfVxMEXTPLcb6k+eG+fpZ91N7C04fs/mqbtPwBfy8znZeYHMvMSYC+AiJgIAuMN7a+tH5dn5k8n/gGvBf68fu6HwMFNP+fA6QqOiL2phmpOycxvZebFVBMsT46Ix0y3r9Tv7ImQ5rnM/HVEvBd4V0QsBT5HFSz2BU4GPpGZv5ilH/c94JiI+C6wHdWH8+Yp2v4KeEZEPB5YB/wp8Jr6uUXAJmAj8MiI2D0zfxoRFwJn1Vdg/Dfwt8DLgMPr/U4BPhsRV1L1cjyXai7IpZMVEBGLqS7nTO57Oefx9THPj4gDM3Oq1yD1NXsipD6QmSdQrRXxJ1QfqNcC76CaP3DMLP6ol1NdEnoF1aJVZ1FdITGZN1NdovlV4L+orsaYWPhq4mqM04AjgB9FxAKqHpX/oJrQ+eP6uWfVvQdk5r8CLwFeQbUi5+PrGqbyPqohlBc2zgnJzLuAFwH/m+pKFEmTGBgfH996K0mSpCb2REiSpCKGCEmSVMQQIUmSihgiJElSEUOEJEkqYoiQJElFDBGSJKmIIUKSJBX5/6b5MONp8NieAAAAAElFTkSuQmCC\n", "text/plain": [ "