notebooks/cda_4_gof.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "cda_4_gof.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyPuYxwfUiwphvEfVfKiJPca",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/DepartmentOfStatisticsPUE/cda-2022/blob/main/notebooks/cda_4_gof.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Setup environment"
      ],
      "metadata": {
        "id": "mkzPNRw1Sk4k"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Python libraries"
      ],
      "metadata": {
        "id": "WM8_gHwDTp8P"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "wxOH_c_TSKmM"
      },
      "outputs": [],
      "source": [
        "import scipy.stats as st\n",
        "from scipy.optimize import minimize\n",
        "import numpy as np\n",
        "import pandas as pd"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Setup R via Python"
      ],
      "metadata": {
        "id": "NX1jMDufT2W0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%load_ext rpy2.ipython"
      ],
      "metadata": {
        "id": "9yUv3YJ7SlnH"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "install.packages(\"vcd\")\n",
        "install.packages(\"fitdistrplus\")"
      ],
      "metadata": {
        "id": "tQ2Pq4mdT4vF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "library(vcd)\n",
        "library(fitdistrplus)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fc26Q3CLUGev",
        "outputId": "b5d65487-145c-46ff-e78b-4e5155b60503"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "R[write to console]: Loading required package: grid\n",
            "\n",
            "R[write to console]: Loading required package: MASS\n",
            "\n",
            "R[write to console]: Loading required package: survival\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Setup Julia via Python"
      ],
      "metadata": {
        "id": "QFZIqeDEUJ9Q"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%bash\n",
        "wget -q https://julialang-s3.julialang.org/bin/linux/x64/1.7/julia-1.7.2-linux-x86_64.tar.gz\n",
        "tar zxvf julia-1.7.2-linux-x86_64.tar.gz\n",
        "## python's module\n",
        "pip install julia"
      ],
      "metadata": {
        "id": "kHQMsiAWUKtS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import julia\n",
        "julia.install(julia = \"/content/julia-1.7.2/bin/julia\")\n",
        "from julia import Julia\n",
        "jl = Julia(runtime=\"/content/julia-1.7.2/bin/julia\",compiled_modules=False)\n",
        "%load_ext julia.magic"
      ],
      "metadata": {
        "id": "iXt4bjxlUQxp",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b32ab48f-767d-49fc-c16b-4a97ddf5a20b"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "\n",
            "Precompiling PyCall...\n",
            "Precompiling PyCall... DONE\n",
            "PyCall is installed and built successfully.\n",
            "\n",
            "PyCall is setup for non-default Julia runtime (executable) `/content/julia-1.7.2/bin/julia`.\n",
            "To use this Julia runtime, PyJulia has to be initialized first by\n",
            "    from julia import Julia\n",
            "    Julia(runtime='/content/julia-1.7.2/bin/julia')\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Initializing Julia interpreter. This may take some time...\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/julia/magic.py:166: UserWarning: Incompatible upstream libraries. Got ImportError: No module named 'prompt_toolkit.lexers'\n",
            "  warnings.warn(template.format(err))\n",
            "/usr/local/lib/python3.7/dist-packages/julia/magic.py:174: UserWarning: Incompatible upstream libraries. Got ImportError: cannot import name 'Completion' from 'IPython.core.completer' (/usr/local/lib/python3.7/dist-packages/IPython/core/completer.py)\n",
            "  warnings.warn(template.format(err))\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "using Pkg\n",
        "Pkg.add(\"Distributions\")\n",
        "Pkg.add(\"Optim\")\n",
        "Pkg.add(\"HypothesisTests\")\n",
        "Pkg.add(\"StatsBase\")\n",
        "using Distributions\n",
        "using Random\n",
        "using HypothesisTests\n",
        "using Optim\n",
        "using StatsBase\n",
        "using DelimitedFiles"
      ],
      "metadata": {
        "id": "jsOjnBvnUSV_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "Pkg.status()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "EbYj69bBlgp4",
        "outputId": "4efaec70-2f49-48d0-eef0-098a9d2980bc"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "      Status `~/.julia/environments/v1.7/Project.toml`\n",
            "  [31c24e10] Distributions v0.25.52\n",
            "  [09f84164] HypothesisTests v0.10.8\n",
            "  [429524aa] Optim v1.6.2\n",
            "  [438e738f] PyCall v1.93.1\n",
            "  [2913bbd2] StatsBase v0.33.16\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Solution in R"
      ],
      "metadata": {
        "id": "XxdM-ZJOYIn-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "set.seed(1)\n",
        "n <- 10000\n",
        "x <- rnbinom(n, mu = 3, size = 2)\n",
        "x_fitted <- fitdistr(x = x, densfun = \"negative binomial\")\n",
        "x_fitted"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fqYbxI4LYJdR",
        "outputId": "b7b555a6-e8c3-4f23-bc58-291b0c44a5cd"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "      size          mu    \n",
            "  2.01336939   3.00049449 \n",
            " (0.05148793) (0.02733434)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "cat(x, sep = \"\\n\", file = \"nb_sim.txt\")"
      ],
      "metadata": {
        "id": "yftkdXjFr-1S"
      },
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Small note regarding the `vcd` package. It gives correct answers but in a slightly different way.\n",
        "\n",
        "1. `goodfit` function returns a list with `fitted` (expected) frequencies under a given distribution using `pdf` of this distribution\n",
        "2. `summary(goodfit object)` returns a $\\chi^2$ and $G^2$ tests using `cdf` of a given distribution. \n",
        "\n",
        "Thus, if we use `object$fitted` created by `goodfit` to calculate $\\chi^2$ and $G^2$ statistics by hand we will get slightly different statistics values than the one returned by `summary.goodfit()` function. Thus, both are accepted results. \n",
        "\n",
        "In this material, we use the first approach in Python and Julia."
      ],
      "metadata": {
        "id": "NdkqM9lW0Uqa"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "gof_nb <- goodfit(x, \"nbinomial\")\n",
        "gof_po <- goodfit(x, \"poisson\")\n",
        "gof_nb"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wxcpCsPrYUMA",
        "outputId": "dc4ddfb6-3b7e-44d8-ffd1-549c2e7ca11c"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "Observed and fitted values for nbinomial distribution\n",
            "with parameters estimated by `ML' \n",
            "\n",
            " count observed       fitted pearson residual\n",
            "     0     1583 1592.9576699     -0.249491416\n",
            "     1     1969 1919.3226642      1.133924923\n",
            "     2     1699 1730.5759462     -0.759033671\n",
            "     3     1356 1385.4754570     -0.791883200\n",
            "     4     1025 1039.1753710     -0.439733965\n",
            "     5      748  747.9232225      0.002807408\n",
            "     6      565  523.1819782      1.828256376\n",
            "     7      357  358.4183143     -0.074916483\n",
            "     8      244  241.6616132      0.150422398\n",
            "     9      164  160.9035748      0.244105604\n",
            "    10      111  106.0489061      0.480781851\n",
            "    11       54   69.3104364     -1.839027630\n",
            "    12       45   44.9808274      0.002858694\n",
            "    13       24   29.0166287     -0.931297590\n",
            "    14       23   18.6216166      1.014623922\n",
            "    15       12   11.8967720      0.029928364\n",
            "    16       11    7.5704166      1.246468975\n",
            "    17        6    4.8004960      0.547467895\n",
            "    18        2    3.0345408     -0.593883338\n",
            "    19        0    1.9128455     -1.383056572\n",
            "    20        1    1.2027237     -0.184850942\n",
            "    21        0    0.7544899     -0.868613754\n",
            "    22        1    0.4723144     -0.226411689\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "summary(gof_nb)\n",
        "summary(gof_po)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "sBclqRTz1dO2",
        "outputId": "332f2bb7-2c4d-41ad-95c9-a964083041da"
      },
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "\t Goodness-of-fit test for nbinomial distribution\n",
            "\n",
            "                      X^2 df  P(> X^2)\n",
            "Likelihood Ratio 21.27568 18 0.2657337\n",
            "\n",
            "\t Goodness-of-fit test for poisson distribution\n",
            "\n",
            "                      X^2 df P(> X^2)\n",
            "Likelihood Ratio 5112.284 19        0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Note that we may use $\\chi^2$ test but the degrees of freedom and p-value is incorrect as it does not account for number of parameters. For the $\\chi^2$ we also use only the non-zero count values."
      ],
      "metadata": {
        "id": "o0fegsa1aloS"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "flag <- gof_nb$observed > 0\n",
        "\n",
        "chisq.test(x = gof_nb$observed[flag], \n",
        "           p = gof_nb$fitted[flag]/sum(gof_nb$fitted[flag])) "
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "j1SwzwM5oYZa",
        "outputId": "01df2d3d-6d6f-47cb-ae3b-fa27f2db7c7a"
      },
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "\tChi-squared test for given probabilities\n",
            "\n",
            "data:  gof_nb$observed[flag]\n",
            "X-squared = 14.509, df = 20, p-value = 0.8038\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "chisq.test(x = gof_nb$observed[flag], \n",
        "           p = dnbinom(x = gof_nb$count[flag], \n",
        "                       size = x_fitted$estimate[1],\n",
        "                       mu = x_fitted$estimate[2]), \n",
        "           rescale.p = TRUE) "
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UZhmh2UbZGtU",
        "outputId": "5e0076ab-bbaa-4753-ea38-b3e5345a692c"
      },
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "\tChi-squared test for given probabilities\n",
            "\n",
            "data:  gof_nb$observed[flag]\n",
            "X-squared = 14.509, df = 20, p-value = 0.8038\n",
            "\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We need to correct p-value"
      ],
      "metadata": {
        "id": "ptlaw25W1ztl"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "pchisq(14.509, 18, lower.tail = FALSE)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qiyHBCTSxLLo",
        "outputId": "33575d9d-cd5d-4a1e-d2b2-9d1dd009b37c"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1] 0.6953586\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%R\n",
        "rootogram(gof_nb)\n",
        "rootogram(gof_po)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 977
        },
        "id": "bP1dbwrfYqE7",
        "outputId": "8e9e69e5-197c-4831-a0b3-5cdea0a700f4"
      },
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3dd0DUZh8H8N+x9xAQFVBQBBwgCiooKg4UreJCW/cWF85ibfWt2taiddXWurXYKrW4Rx1YFdxVq1a0AiKCbNl7HHf3/nEtRW4QlLvk4Pv5K/fkSfLLPfAl5HIJTyQSEQAAcI8a2wUAAIB0CGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYD+V0EBVVSwXQQAwH8afUAXFVFQEPF4ZGxM2trUty/9/DPbNQEAEBHxRCIR2zWwp7CQhg2jyMia7Z98QuvXs1EQAMB/GndAr1xJX38tfdb9++TurtxqAADe0rhPcchKZyI6ckSJdQAASKHBdgEsyM3N/fTTT3lC4U7ZfR4fP76zsNDW1vbTTz9VXmUAANU0xiPojIyMpKSk4aNGyenT1M5u5MiRFy5cUFpVAAA1NMaAJiIDAwNra+tCLy+ZPXr2tLa2VldXV2JRAABvaaQBLZY1bZqsWQX9+yuzEgAASY06oIu7dcucNatGY3mrViJ19ZbLlqkVFbFSFQCAWGP8kLAKTyg0vH6db2mZNX26VlKSUFe3tF27Qm9vg1u3Wi5bZjdjhpGpKds1AkDj1agD2vTXX3ViYl5/+22NExqF3t6vt22zWbz4Ow0NSkuj5s3ZqhAAGrPGe4pDPT/fcseOou7dpZ5uLvTySty1y6KsjPr1o5QU5ZcHANB4A9ry22/VCgvTV6yQ1aHY3X2FiwulppKXF8XHK7M2AABqtAFtm5treuJE9sSJZfb2crpFGRvT1atUWEh9+9KLF0orDwCAGmlAi0TTHj8WmJpmzplTe2c3N7p+nfh86tWLoqIUXxwAwD8aY0Abnz7tlJmZvnSpwMCA0QLt29O1a6ShQd7edP/+P435+fTsGZWVKa5OAGjkGl9AFxZabNkSa2aWN2xYHZZydKSbN8nEhAYOpG++od69ycSEOnYkXV3y9aVHjxRWLgA0Xo0voNes0cjOPti5M/F4dVvQ1pauXSNtbfrkE7px47/2S5eoSxe6d69+ywQAaGQB/fw5ff99nr9/XJMm77K4kRFlZEif9ckn71MXAICkRhbQS5eSgUHmokXvuHh4uMxZERH06tU7rhYAQJrGFNAnTtDFi/Tll4K6Hz7v3bvXx8dn96pVcvosGzNm4MCBN6qf/QAAeA+NJqBLS2nZMnJxISaX1kl48ODBzJkzB4weLafP5MDAHj16REdHv2uJAABvaTT34tiwgRIS6No1etdbPBsYGAjk3D+aSMPVVTc9/d1WDgAgqXEcQb9+TRs30rhx5O39PqupaNUqd8QIqbNSP/tMhLv7A0C9ahwBvXgxEdH69e+/ptRVq/L8/Go05o4cmTNu3PuvHACgukZwiuPKFTp5koKDqWXL91+ZSFs7ed26zBkz9B4/1szIqLC2tvzuO50XL0gkqvOF1QAAcjX0gK6spMWLyd6eliypx7WWt25d3rq1eJpXWWn1+edGEREFffvW4yYAABruKY6KCiKibdvo6VPato20tRW0nbzhw8tbt2763XckFCpoEwDQOCkpoCMiIng8XtUlaBcvXnR2djYzM/P19U2v3ysfrl8nPz/i8Uhbm3r0oJUrydeXhgypz028TaSm9mb+fJ24OJPz5xW3FQBohJQR0OXl5UuWLLG0tBS/zM/PnzRp0p49ezIyMtzd3QMDA+ttSz/9RH360Nmz/7y8c4fKyykjgzIz620T0uT7+JQ6Ozf9/nt1gUChGwKARkUZAR0cHOzn52dubi5+GR4e7ubm5unpqaGhERQUdPbs2fLy8nrYTHo6TZkipf3RI5L92JT6weNlzJunlZrqUnUzUgCA96bwDwljY2OPHz9+//7948ePV7U4ODiIp42NjU1MTBITE6taqmRnZyckJNRoLC8vt7Ozay71Ka4nTsgs4sAB2rmTtLTecR8YKPLyKu7atceVKykDBihuKwDQqCg8oOfOnbt161YdHZ2qlpKSEl1d3aqXenp6JSUlkgtGRUVdvHixRmNMTEyHDh2++uqr6o2hoaFPnjwZevWqnO/5bVi0KNfYeM6cOba2tu+yGwxkBAa2njy5w9WrNH++gjYBAI2KYgP64MGDzZs3H/D2QaW+vn5qamrVy6KiIgNpTzbx9vb2lvjiX1hYWFZWVo3Gw4cPf/TRR2YxMST7DEO7zp1/jYx8/vy54gK6pHPnl+3adbp0ifLzydhYQVsBgMZDseegT506denSpWbNmjVr1iwmJqZXr17nzp1zcnKK+vfhfikpKcXFxe8fmu3atdPr10/W3JJOnVr36FH1KaXiRA4apFVSQps2KXpDANAYKDagT548mZmZmZ6enp6e7ujoeOPGjaFDh/r4+Pz9999Xr16trKxct26dv7+/hkY9HMgX9e5d3KWL1FmSX85WkMxmzeLd3WnrVpn39QcAYIyFL6oYGhqGhoYGBgZaWlomJSVt2bKlXlYrUldP+frrks6da7RnTp+eM3ZsvWyCiT+HD6eKCgoOVtoWAaChUt5XvZ8+fVo17ePj8+zZs3rfRIWV1av9+40uX9Z99kwjJ6e8devC3r3LHB3rfUNy5DdtSlOn0s6dtGgR2dkpc9MA0MA0tHtxiDQ184cMyVfkVwdrt2YNHTpEX31F+/ezWQYAqLiGey8OFrVoQXPn0sGD9Pw526UAgApDQCvGp5+Svj6tXs12HQCgwhDQimFuTkuW0LFjdO8e26UAgKpCQCvMxx+ThQUOogHgnSGgFcbAgJYvp4sXKSKC7VIAQCUhoBVp/nyysaFVq9iuAwBUEgJakXR0aNUqunWLfvuN7VIAQPUgoBVsxgxydKRPP8UDsQCgrhraF1U4R12dVq+m8eNpyhTS0CANDXJxoY8+IgsLtisDAK5DQCueSEREdOjQfy0LF9KFC+Try1ZFAKAScIpDwR4+pAkTpLQPHkwSz4sBAKgOAa1gcu7Vt327EusAANWDUxyK8uTJEz6f7/D8uaGMDoXXr8f++aeenl67du2UWhkAqAgcQStESUnJ0KFDd+zYYfjwoaw+qcnJO3bsGDx4sDILAwAVgoBWCJFIZGdnt2TJkiIvmU+ytezcecmSJU2bNlVmYQCgQhDQipUn+1INpT2ICwBUFAJasfKGDcsfNEiyPeejj4p69lR+PQCgQvAhoYKpqSV9802xu7vx5cv69+4RkVBHp9LCInXlSrYrAwCuQ0Arnppazkcf5Xz0EQmFxOM1CQtr8dVXeg8flsh4BjkAgBhOcSiRmhrxeLkjRlSam1uEhLBdDQBwHQJa2UTa2jkffWQYEaEdF8d2LQDAaQhoFmR/9JFQV9f84EG2CwEATkNAs0BgbJw7apTJuXOa6els1wIA3IWAZkfWlClEZBYaynYhAMBdCGh28Js1yx80qElYmKFAwHYtAMBRCGjWZE2frlZSMvLNG7YLAQCOQkCzpszBoahHj/EZGVRWxnYtAMBFCGg2ZU6b1oTPp8OH2S4EALgIAc2m4u7dn+rr08aNeKQsAEhCQLMstFkziomhM2fYLgQAOAcBzbKrpqZkb0/BwWwXAgCcg4BmmZDHoyVL6N49unmT7VoAgFsQ0BwwfTo1a0YbN7JdBwBwCwKaA3R0aO5cOnuWnj1juxQA4BAENDcsWED6+rRlC9t1AACHIKC5oUkTmj6dfv6ZkpLYLgUAuAIBzRkff0xEtH0723UAAFcgoDnDxobGjqVduygvj+1SAIATENBcsnw5FRbS7t1s1wEAnICA5hIXFxo4kLZuxe2TAIAQ0JwTFEQZGbh9EgAQAppz+ven7t1x+yQAIAQ0Fy1ZQjExdPYs23UAAMsQ0Nzj74/bJwEAEWmwXQBIUFenwEBatIh8fEggoKZNydubpk4lHR22KwMApUJAc09mJh09SkT0++//tPz6K4WG0rFj1LQpi3UBgJLhFAf3zJ4t5dajN27QrFlsVAMArEFAc8yrV3TqlPRZZ85QXJxyqwEANuEUB1ckJCTcv3/f6v79HrL73Nm5M9nDw83NrXXr1sqrDABYgoDmikOHDsXGxg4tL5fTJ+X16/M5OU+fPl27dq3SCgMAtiCguUIkEvXp06dz06YUFiarT6cJE4rz8+Pj45VZGACwBeeguaXczq6oh/STHEUeHuX29kquBwBYhIDmnNTVq0tcXKS2K78YAGARAppzKlq0SNy5801AQJGHBxGVOjsTUcbixRXW1myXBgBKhXPQXCQwMnqzYEHVS7spU5ocOZI1ebJIU5PFqgBAyXAErQKyZs7UTE83Pn+e7UIAQKkQ0CqgsFev0nbtLPbvxz1IARoVBLRqyJo6VfvVK6OICLYLAQDlQUCrhgJf34pWrSz27WO7EABQHgS0ahCpqWVNnqwbFdXs+XO2awEAJUFAq4zcESMqmzZ1/u03tgsBACVReEAfOXLEwcHBxMSkX79+cf/ejO3ixYvOzs5mZma+vr7p6emKrqFhEGlpZU+Y0OLp0xZpaWzXAgDKoNiAjo2NXbBgwfHjx7Ozs93d3efNm0dE+fn5kyZN2rNnT0ZGhru7e2BgoEJraEiyx46t0NPzunWL7UIAQBkUG9BaWlqHDh1ydnZWV1cfNWpUTEwMEYWHh7u5uXl6empoaAQFBZ09e7Zc7i3coIrQwCCmf//2z59TbCzbtQCAwin2m4S2tra2trZEVFBQsHv3bj8/PyKKjY11cHAQdzA2NjYxMUlMTKxqqZKenv706dMajVFRUSYmJgqtmeP+9vFxunBBc+NG2ruX7VoAQLGU8VXvoKCgTZs29erV69SpU0RUUlKiq6tbNVdPT6+kpERyqRcvXvxe9VC+f8XExHTs2FGh1XJcmZHRI1fXbgcP0uefk40N2+UAgAIpI6A3bty4du3anTt39u3b9/Hjx/r6+qmpqVVzi4qKDAwMJJfq1atXr169ajSGhYVlZWUptlzOu92jR7fHj2nbNtq0ie1aAECBFHsO+q+//rp69SoR6enpLVq06OnTpxkZGU5OTlFRUeIOKSkpxcXF4tMgwFCuiQl9+CHt2kXZ2WzXAgAKpNiATktLmzx58suXL4no8OHDFhYWlpaWPj4+f//999WrVysrK9etW+fv76+hgZvq1dHKlVRaSjt2sF0HACiQYgPa19d30aJF/fv3NzU1/e6778LCwng8nqGhYWhoaGBgoKWlZVJS0pYtWxRaQ8PUrh0NGULbtlFREdulAICiKPzQNSgoKCgoqEajj4/Ps2fPFL3pBm7VKvLwoAMHaOFCtksBAIXAV71VVvfu1KsXbdxIFRVslwIACoGAVmUrVlByMv3yC9t1AIBCIKBV2ZAh1KULrV+PG/kDNEgIaBX38ccUHU1nzrBdBwDUPwS0ihs7ltq2peBgtusAgPqHgFZx6uq0dCndu0fXrrFdCgDUMwS06ps2jVq0oPXr2a4DAOoZAlr1aWvTwoUUHk5//sl2KQBQn/Ad6wZh7lwKDqZ586hXLyovJxcXGjmSzM3ZLgsA3gsCukGIiqL8fLp3j+7d+6dl9mw6cYJGjmS1LAB4LzjFofrS0sjLS0r7qFH05InSqwGAeoOAVn0//CBz1saNSqwDAOoZTnGosNTU1LKysuYREboyOpTFxqbGx2toaLRs2VKplQFAfcARtArr0aPHJ598Ei/7AbI69+598sknnp6eBQUFyiwMAOoFAlqFNW3adPXq1VZdu8rqUOThsXr1aicnJyFu1gGgghDQKi9/6FBZswr691dmJQBQvxDQKq+wV69caZfTFfbunTtmjPLrAYD6goBuCFLWrk393/+K3dzELystLERqaulBQSJ1dXYLA4D3gas4GgQeL2fs2JyxY3l8PgmF6gUFDoMHmx84kPLFF2xXBgDvDkfQDYpIU1OkrV1pYZEzdqzJ6dPaCQlsVwQA7w4B3TBlzpol0tW12LWL7UIA4N0hoBsmgalp9kcfmVy4oCP7KmkA4DgEdIOVNWOGQF+/KQ6iAVQWArrBEhgaZk+caPT77474GiGAakJAN2TZkycLjIwm46NCANVUt4BOSEi4f//+vXv3Xr16paCCoB4JDAyypk71yMrSePCA7VoAoM4YBfTLly/nzp1raWnZuXPnCRMmTJw4sUuXLpaWlnPmzHn58qWiS4T3kT1hQq6Wls7XX7NdCADUWe0BvXnzZh8fn1atWl2/fj03Nzc2NjY2NjY3N/fGjRu2trY+Pj5btmxRQqHwboS6ukdattS4do0iItiuBQDqpvaAfvXq1dOnT1esWOHo6Fi93cHBYcWKFU+fPo2Pj1dYeVAPzlhbC1u0oFWr2C4EAOqm9oDevn27np6eq6vr+vXrJbNYT09v+/btiqkN6keFmlr50qV06xaFh7NdCwDUAdMPCb/66qsXL154eHh069Zt8+bNr1+/VmhZUL/KJ08mOztauZJEIrZrAQCmmAb00KFD9+/fn5aWtmnTpsTExO7du/fs2XPHjh14VIdq0NSklSvpwQM6d47tUgCAqbpdZldUVPTy5cu4uLiSkhIzM7MnT564uLhcuHBBQcVBfZo6lRwdaeVKwtNVAFQE09uNHjt2LDQ09OLFi927dx8/fvzhw4dNTU2J6O+//x48eHBiYqIii4T6oK5Oq1bRpEl04gT5+7NdDQDUjukR9IYNG7y8vOLi4q5duzZr1ixxOhNR+/btx48fr7DyoF6NH0/OzrRqFQkEbJcCALVjGtA3btwoKyszMzMjoqSkpA0bNpSXl4tnBQcHK6o6qF9qavS//1FMDB05wnYpAFA7pgE9c+bMu3fvCgQCIjIwMHjw4MHs2bMVWRgohr8/de5Ma9ZQZSXbpQBALZieg46IiEhMTFRXVyciU1PTX375pWXLloosDBSDx6O1a8nPjw4epBkz2K4GAORhegStrq6ekpJS9TI6OlodDyRVUcOGkYcHffEF/XuSCgC4iekR9Lp169zc3Hr27GliYpKZmXnz5s19+/YptDJQoLVradAg2r+f5s1juxQAkIlpQE+cOLFnz56XL1/Oy8uzsLDYu3dvixYtFFoZKNDAgdSnD335JbVuTa9fk5UVubuTpSXbZQHAW5gGNBHZ2dnhg8GGo2dPioykwYP/a1m6lDZsII06/EgAgEIx/W08f/78l19+mZqaKqh2CW1ycrJiqgIFCwsjyTtEb9lCfD599x0bBQGAFEwDetmyZV9++WXbtm3x2aDKE4noww+lz/r+e1q8mFq3Vm5BACAd04C2sbHxx/eDG4aYGHlzr1xBQANwBNPL7Lp163blyhWFlgKKduLEifbt2wfI/UO786uv3NzcgoKClFYVAMjC9Aj63LlzwcHBhoaGBgYGVY04B61a8vPzJ0yYMMbN7a3PBt82bPZsOze3o0ePKrMwAJCKaUCHhIRo4PP9BqHC2rq4Sxf9hw+lzi3q2ZPKypRcEgBIxTRzXV1diYjP52dlZTVv3lyRJYHCvVmwwG76dMn2jMWLK83NCf8YAXAD03PQb968GTlypL6+fvv27YkoMDDw9u3biiwMFKi4a9dX+/eXdO5cvbGySZOsSZPYKgkAJDEN6I8++sjLyysnJ8fY2JiIJk+evGjRIkUWBopV3K1b/E8/xYSHv9q378XZswm7dmnk5Jj98gvbdQHAf5ie4khMTFy2bFnVy65du+bl5SmmJFAefvPm/ObNiajc1rawd++mO3fmDRnCdlEA8A+mR9AGBgZPnjypehkdHa2jo6OYkoAdaStW8Ph8y++/Z7sQAPhHHe5m17dvXzc3t6ysrBEjRty+ffvHH39UaGWgZBU2NtkTJ5qHhBj36cN2LQBAxDyghw4d+uTJk99++23AgAHNmzffvXu3JW5+1uC8mT3b5Ny59j/8wPPyYrsWAKjL3eysrKxwN7uGTaivn75wofWqVZ64khKAA5gGtJOTU40WkUgUI/+uDqCC8vz89H/5Zcz9+1RQQEZGbJcD0KjV4ZuEVdMFBQXHjh2TjGxoCHi85wEBnosW0ddf0/r1bFcD0KgxDWgPD4/qLwcOHPjBBx8sXbpUASUBywratr1lb++1dStNn04ODmyXA9B4Mb3Mrobi4uLExMT6LQW443jXrqSrS9WufAcA5XuXc9BCoTA1NTUgIEAxJQH7CnR0aNUqCgqi8+cJX10BYMm7nINWV1dv3ry5tbW1QioCjli0iPbto0WLqH9/0tZmuxqAxojpKQ6Parp27co8nc+cOdOuXTsTExNvb+/Y2Fhx48WLF52dnc3MzHx9fdPT09+lcFA0TU36/nuKiyN8txCAJUwDulWrVvayyVoqOTl58uTJ+/bty8nJ8fLymjNnDhHl5+dPmjRpz549GRkZ7u7ugYGB9bMrUO98fGjoUPriC0pLY7sUgMaI6SmOhQsXPnjwYPLkyRYWFhkZGXv27PH09PT19a11wb179/bs2ZOIRo8eLT5PEh4e7ubm5unpSURBQUGWlpbl5eXa+Ceam779ljp0oJUr6cABtksBaHSYBvTBgwcfP36spvbPEffAgQM7deq0YsUK+UtZW1uPGTNGPH358mUvLy8iio2Ndfj34i1jY2MTE5PExEQHicu5kpOTo6OjazRGRUUZGhoyrBnqQZs2tHgxffMNBQRQ9+5sVwPQuDAN6Nzc3MzMzKr7b7x58yY3N5f5Zi5durRr167r168TUUlJia6ubtUsPT29kpISyUVSUlL+/PPPGo3x8fFt2rRhvl2oB6tW0c8/0/z5dO8eqb3jdZkA8A6YBnRAQEC7du369OljbGycl5cXGRnJ/Ib9oaGhX3zxRXh4uPijRX19/dTU1Kq5RUVF1R9EW6V79+7dJQ7ZwsLCsrKyGG4X6oeBAQUH05Qp9PPPNGUK29UANCJMD4hWrVp18+bNQYMGOTg4DBo0KDIycs2aNUwWPH369Pr16yMiIqo+S3RycoqKihJPp6SkFBcX29ra1rlwUKZJk8jLiz75hPLzSSjEU2UBlKMOd7Nr375927Zt6/TQ2Nzc3AULFty4caNZs2ZVjT4+PgEBAVevXu3du/e6dev8/f3xvHCu4/Fo2zZydydXV0pIICLq3Zt69aLPPiM9PZZrA2i4FPvQ2FOnTqWkpDg5Oen8Kzs729DQMDQ0NDAw0NLSMikpacuWLe+1B6Acv/1GItE/6UxE16/TunU0dCgVF7NZFUCDptiHxk6bNk0oFJZVY2ZmRkQ+Pj7Pnj3Lzs4+e/asuAU4LSqKPv9cSvu1axQcrPRqABoLPDQWGAgNlTkrIkJ5ZQA0MnhoLMizadOmgICAP8PCZPa4dWv+rFkLFizIyclRYl0AjQLTgBY/NHbgwIHih8b27t17Pe7m3gicOXNm5MiRzezs5PQZNmpUUlJSSkqK0qoCaCTw0FiQR11d3dramufpSVeuSO1Q5OFhbWMj9Up2AHhPTAN6wIABv//+Ox4a2zjlDRvW4quvpM7KxldXABSG6SmOpk2bHj16VCgUKrQa4Cahnl7ciRMlnTrVaC/p3LnQy4uVkgAaA6ZH0NHR0dOmTZs6daqJiQmPxxM3JicnK6ww4Jaytm0T9u83Cg/XiY5WLy4ua91aOzGxSViY4Y0bhb16sV0dQMPENKD379+vqamp0FKA44Ta2nnDhtGwYeKXvIoKvUePrNaseXHiBLuFATRUtZ/i6NixIxF17ty5Y8eOGzZssLa27vgvxZcH3CXS0kr++mv13NzmuJ4HQDFqD+jqN2W+fPlyGW6UA/8qc3LKnDXL5Ny5bq9fs10LQAOE2/vCe8mcPbu0Y8fp9+9rZGezXQtAQ4OAhvciUldPXrdOu7LS6osv2K4FoKFh9CFh1dUaQqEwLS2tsrJS/JL5s72hAStv3fq4i8u4a9fo0CGaOJHtcgAajtoDWiAQ2NjYVL3s0qVL1bRIJFJIUaBqfnNy8lNT01+wgPr0oWo/LQDwPmo/xcGXTQn1gUoQ8XjJX31FAgHNmEH4sw1QT2oP6KVLl1ZWVmpIQ0RlZWXMH04IDViFlRWtX0+XL9O+fWzXAtBA1B7QNjY2HTt23Lx5c1xcXPX2ly9fbt68uWPHjjgTDf+YN48GDaIlS+jlS7ZLAWgIag/ooKCg8+fPR0dHe3h4mJqaOjo6Ojo6NmnSpFu3bs+fPz9//nxQUJASCgUVwOPR/v2kpUVTpxJu2wLw3hhdxeHg4PD555/v2bPn5cuXWVlZRGRmZpaTk9O9e3cFlweqxsqKNm2iGTPo++8J574A3g/Te3EMHjz46dOn9vb29vb2RFRUVNStW7fc3FxF1gaqafp0OneOVqwgHx9q357tagBUWO2nOEJCQpo1a/b333/rVGNsbNy1a1cl1AcqadcuMjKiGTNIIGC7FAAVVvsR9NSpU6dMmTJq1Ki9e/dWNWpqaoof7w0gRdOmtG0bjRtHixaRmRlFR5O5OXXrRmPHkq4u28UBqAxGpzh4PF5hYaG5ubmiq4GG48MPac0a+uGH/1p27KA9e+jHH8nBgb2yAFQJnqgCirFrF8XE1Gy8fZumT6eKCjYKAlA9TANa/EQVQ0NDKysr638ptDJQbYcOSW+/dYvOnFFuKQCqCk9UgfpUUVERFRWlkZvb6fZtWX0yzp5NtrOzsLBo2bKlMmsDUDlMA7pz584KrQMahhs3bixfvrxv27Y1ny9bTdSDB6wStrcAACAASURBVKE//JCSknLp0iXlVQaggpgGdKtWreQcQdf4Fjg0WiKRqEePHtNmzqRff5XVp8OIEQvHjFm+fLkyCwNQRUwDeuHChQ8ePJg8ebKFhUVGRsaePXs8PT19fX0VWhyoKJG6es6HHzaRkdH5Pj5KrgdARTEN6IMHDz5+/FhN7Z8PFQcOHNipU6cVK1YorDBQbRkLF2q/eKH/8GGN9txRoyrNzenfZz4AgBxMr+LIzc3NzMysevnmzRt8zxvkEBgZJe7cmTlrVrG7u7iloF+/EhcXk3PndJ88Ybc2AFXB9Ag6ICCgXbt2ffr0MTY2zsvLi4yMxG2gQT6hnl7GwoXVW9QLCtqMG9dq8eIYWRfhAUA1TAN61apVo0aNun79ek5Ojqmp6RdffOHi4qLQyqDhERgZJX7/fZsJE2yDgjQtLNguB4DrmJ7iKC0tzc3NnTNnzqJFi4qKisLDw4uKihRaGTRI5a1bJ3/1lX5U1KLnz9muBYDrmAb0zJkzT506RUTz588PDw9/+vTpzJkzFVkYNFgF/ftnTJs2MDWV9uxhuxYATmN6iuP27dsvX74sKys7efJkQkKCiYlJ27ZtFVoZNGBp8+cnnDnTfcECcnKi3r3ZLgeAo5geQaurq/N4vCtXrri5uZmamhJReXm5IguDBk1NLdjZmeztaexYSk5muxoAjmIa0J6enoMGDQoICFiwYAERrV692tXVVZGFQQNXoqFBJ09SWRkNH06lpWyXA8BFdbhZ0pkzZ5o3b96zZ08iMjc3P3DggCILg0bA0ZEOHqRRoygggH76ie1qADiHaUBraWn5+/tXvVz49vWtAO9o+HD63/9o7Vrq1o0WLGC7GgBuYXqKA0BRVq8mf39asoQiItguBYBbENDANh6PDhwgR0fy96dXr4iIKispO5vtsgDYh4AGDjA0pBMnSCCgAQOoXz/S1CRzc+LxaMYMyshguzgA1iCggRscHGjiRIqPp2vX/ms8cIBGjKCUFPbKAmATAhq4ISWFtm+X0n73Ln3yidKrAeAEBDRww4kTMmcdPkxlZUosBYArmF5mB6Age/bs+fPPP8fcuTNAdp+V06ZlGRl9/PHHuMEANCo4ggaWnT59evDgwXYdO8rp09/Pr7KyEo++hMYGAQ3ss7a2Vpd9y6QSV9emzs7GxsbKLAmACxDQwAlFPXoUd+8udRYeMguNFgIauEFNLSk4uLhbtxrNQh0ds7AwjWrPwwRoPPAhIXBFpYXFq717DW7f1nv2TD07u7x16yIvL42cHNvZs+2mTzccNIjtAgGUDQENXKKmVuTlVeTlVdVQYW2dsHOnbUBAwIkTf/frx2JpAMqHUxzAdSWdO7/+7rumubnun39OhYVslwOgPAhoUAFFHh4Hhg0zio8nX1/C04qh0UBAg2qItrV9vHw53btHI0fii4XQSCCgQWVk9OxJoaF07RqNGEF4JCY0AghoUCljxtDevXT5Mo0bR5WVbFcDoFgIaFA106bRtm108iSNG0cCAdvVACgQAhpU0IIFtHUrHTtGM2fSjz/SBx8Qj0fe3jRxIt27x3ZxAPUG10GDalq8mHJy6MsvKSTkn5bISCKiw4dp716aOZO9ygDqDY6gQWWZm0tvnzWLcN87aBAUHtCVlZXLly9XU1PLysqqarx48aKzs7OZmZmvr296erqia4CG6fhxmbMOHlRiHQCKovCA9vf319PTU1P7b0P5+fmTJk3as2dPRkaGu7t7YGCgomuAhqSkpGTTpk3frF9P16/L6hN76tSGDRt+/fVXZRYGUO8UHtBr1qxZs2ZN9Zbw8HA3NzdPT08NDY2goKCzZ8+W45pWYCwpKen48eMampryOmloaGpqfvvtt8oqCkAhFP4hoaura42W2NhYBwcH8bSxsbGJiUliYmJVS5VXr149ePCgRuPdu3ctLS0VVCqoChsbm4GDBhWHhenLuGbD2Nt74MCBZ8+eVXJhAPWLhas4SkpKdHV1q17q6emVlJRIdsvOzo6Pj6/R+ObNGzxZA8RyRo+WFdDFnTsruRgARWAhoPX19VNTU6teFhUVGRgYSHZzd3d3d3ev0RgWFlb9w0ZozPKHDNGNiTE/cKBGu8DIqNXCha+3bmWlKoB6xEJAOzk5Xbp0STydkpJSXFxsa2ur/DKgAUhfsqTA29v44kXthAShjk6Zk1Pu6NE8Pr/VvHm2AQED7e3ZLhDgvbAQ0D4+PgEBAVevXu3du/e6dev8/f01NPB9GXhHJZ07l0ic0Ij/+eeWS5eu+OMPWrGCgoOJx2OlNoD3pNirOLKzs3V0dHR0dAQCgbW1tY6OTkZGhqGhYWhoaGBgoKWlZVJS0pYtWxRaAzRCAiOjhN27z7ZoQRs20Ecf4fakoKIUe+hqZmZWJu13w8fH59mzZwrdNDRyInX1rQ4OfsuX09KllJxMp06RhQXbRQHUDb7qDQ3aokX066/06BF5elJMDNvVANQNTv5CQ+fvT9bWNHw49ehBJ06QhwcdOUK3b1NyMrVtSz4+NGQITlIDNyGgoRHw8KC7d+mDD8jHh+zsKDb2v1nbttGoURQaStra7NUHIB1OcUDjYGdHN2+SoeFb6Sx24gStWMFGTQC1QEBDo/HiBeXkSJ/17bdUXKzcagBqh1Mc0PAlJia+ePHC5swZR9l97v34Y4GTU5cuXZo0aaK8ygDkQkBDwxccHFxWVjYiIUFOQP9x48aF8+dHjhw5a9Ys5VUGIBcCGhoFf39/5+Tkfx6LJY3P/Pmld++KRCJlVgUgH85BQ2NR1K1bSadOUmcJdXS0EhKUWw5A7RDQ0FiItLTSly2TbC9r21ZgZmY3fXqfCxfUBALlFwYgC05xQCNS0rlz7G+/Wezbp/X6tf6ffxZ361basWPmjBkidfVmmzZ1P3bMLiGB+vcnicdHALACAQ2NS0XLlilffCHZnrp69XVdXd8TJ8jVlYKDadEi5dcGUANOcQD8I7Zjx6Nr1lDfvrR4MQ0eTDWeN19aShUVLJUGjRQCGuA/pUZGdO4c7d5N169Tp0507hzx+bR+PfXpQ3p6pK1N3t60fTsJhWxXCo0CTnEAvI3Ho9mzycuLJk4kPz9q2ZISE/+bGxlJkZF07x799BN7JUJjgSNoAGnat6c//qChQ99K5yo//0ynTyu9Jmh0ENAAMmhqUnm5zLmhoUosBRopnOIAqEkgECxatIjP5+8OD5fVJ/bmzc0BASYmJhs2bFBmbdCo4AgaoKaysrL79++PHDkyx9lZVh8zC4uRI0dWPZ8eQBEQ0ABS6OjoWFtb89q0kdVBVyRq2aSJlpaWMquCxgYBDSBT1tSpsmbpRUW1HTZsSFYW4f5KoDAIaACZytu0Sfz+e8n2VyEhLw8d4jdrtubVK/LwoDt3lF8bNAb4kBBAnkJv7+c3bhhduaLz4oVIQ6PMwaGgf3+hvj4Rxf/887ERI9a8fk09e5K/P23aRC1b/rPYX39RZCTFx5OdHfXuTZ07s7kPoLIQ0AC1EJiY5I4eLWWGmtp5c/M1167Rxo20YQOdP08ff0zLltHixXTgwFs9J0+mvXsJJ6yhjnCKA+D96OvTmjUUE0OjRtHatWRlVTOdieinn0janU4B5ENAA9SHli3pp5/o9GkqLJTeYft2evNGuTWBykNAA9SD8vLy+Pj49KwsOX0yTp2Kj48vl/PtRIC3IaAB6sHWrVvHjh27/7vv5PQJ2b593Lhx69evV1pVoOrwISFAPaisrJwzZ05fAwMaN05Wnyk+Ps07dYqLi1NmYaDSENAA9aa0ffsSV1e9x48lZ4k0NZtt2fKBrW1kp04kEJC6uvLLA5WDUxwA9UdNLWPJEqlzErdvT163TrOsbNTp0+TkRNu2UVnZf7Pv3aMZM6hfP+rXj8aPJ9ziA4gIR9AA9au4S5e4o0ctt283jIwUtxT27v1m/vzS9u2J6JqVldqFCxPi42nxYlq/ngICaMkS+uknWrjwrbX88gstW0abNim/fuAUBDRAPStzckrcvp2EQo2sLIG5uUjtv/9TRTxejIMDHTlCN2/Shg20di1t2kTFxVLWsnkzeXnRiBHKqxu4B6c4ABRDTa2yadPq6fwWLy86e5YePiQrK5lr2LdPQaWBqsARNAA7jh49mpubO5THayGrx2+/7dmzR1NTc8qUKWqygh4aNIw6ADs+//zzoqKiSoFATp+ioqKtW7fm5+crrSrgFBxBA7DD0NBw4MCB+g8ekOwro+fu369TWcnLziZT05rzHj6k8HCKjaUWLahHD/L1JRxlNzgIaAA25Y4da/bLL1Jn5Ywfr/fo0bzYWOrYkfz8aNIkGjyYNDRIKKQFC2jnzrd6DxhAR46QmZkyigZlwZ9cADaV2dunfPGFZPubuXNTP/00LixspodH2Zw5FBFBfn7UqhUtWkSLF9dMZyL6/XeaMkUZFYMS4QgagGW5I0eWduxoevSo9qtXPKGwvGXLvGHDSrp0Ec9N0NcvXb1a55tv6NIl+vln2rmT+HzpK/rtN7p/n7p2VV7poGAIaAD2lbVtm/bZZ/J6aGmV+fikduigO2xY88mTZfXKPncu38ysSZMmJiYm9V8lKB0CGkA1hIWFbd++3cfQcJ3sPpFhYfsfPGjWrNn+/ftrzisvpzNn6OFDKiigDh3Iz4+srRVZL9QDBDSAahAKhcOHDx/dqxddvSqrz6jo6IElJdeTkujmTerZk3i8f2Y8f04zZ9Lt2/91nT+fdu2igAAFVw3vBR8SAqiSyiZNCr29Zc1NX7Kk3MzM99kz6tWLbG0pIIDOnqXcXGrf/q10Fpszh8LDFVotvCcENICKSVu+XGp74nffZU2f/iA4ePGHH9LBg9S9Ox05Qn5+1KyZzHXhfkzchoAGUDEVNjbRV67kVruPUmHfvi9DQwv79hW/LNLRocmTKSws8+nT9H37+HJu93H5cvzLl69evRKJRFLmJiXRnj20ZAkFB1N4OAmF9bwnUBucgwZQPZVNm6Z8+WXKl1+qlZYKdXVldfMbO9ba2nopn+8pe1WHp0//NSXlyKlTHTt2fGvGt99SjXtb9+tHP/9MLWTeOwTqHY6gAVSYnHQmIh0dndWrV9sPGiSnz/+uX/8rPr7Nhx/SihV09izl5RERhYXVTGciunqVJk7EcbQyIaABGrj8wYNlzcqaPDlh797THTsKjIxo61by8yNzc3J3l5LOYteu0fnzcjeWT9ev0/Hj9PQpyb0PFDCBgAZo4MocHd/MmSPZXuzm9mbevCIPjxMuLvG7d1N2tujixbw5c8qFQkpNlbW2vHPn4uPjpdxgj8+noCAyMaE+fcjfn5ydSe4VgcAEAhqg4Xszf37Sxo3F3btXtWRNnZq4Y4dQX/+/TgYGkdraXS9dmmFoKGdVOaGhJ4YP/2TIEKqsfGvGvHk1rwm5c4f696c7d+pjDxopfEgI0Cjk+/rm+/ryBAJeRYWsM9cCgcDX13fO9On0751AJNny+R8/fUpEpKdHbduSmxu5uZGxscznv6xeLe9q6zt36Ndf6flz0tcnV1eaMUPeI2YaHwQ0QCMiUlcXyf1ckYhEmpp5w4aZnD0rdW7ciRPCiooTq1Z9NngwPXxI587Rzz/LW93ly5SXR1LvDfLpp7R+/X8vT56k1avp0iUaOLC2/WgscIoDAGpKX7y4xMVFSntQUIWNTbmt7bVmzeibb+j33xeMHz/K1fWW5PMEqjnUocO8Vq0eHzv21seGR468lc5VBg2i9HSZ6+LzKTSUFiyg4cMpKIguXWK8TyoJR9AAUFNl06avv/uu2datJqdPi1uKu3XLHjeuYMCAGj0TEhJW79nTMjSUvv1W1tompqZOJKIxY0hbm9q2JScncnSkkydlbj4khFaskNKenk4ffkjXr//XsmkTjRlDhw6RllYtu1ReTmpqpKlZSzeOQUADgBSVZmbJX32V/NVXGtnZAiMjkdxoK/TyspQR0EWenkmbNl3+/ntvS0sXbW1RdHTp3bu6J0/yZF+El3X69GM3N9MmTdzc3N6aMW3aW+ksdvQo2djQ5s0ydqOStm6lc+f+WbBPH/rwQ5oz57/bSHEbAhoA5Klk8BitMkfHXH9/02PHJGe9WbBAYGSU0Lx58oABLkOGlJaU/PDDD1RREbRqlay1md+923fIkGJzc3JxIVtbatWKbG2ptJQuXpS+wJYt9MUXVP2KFDGBgMaMoVOn/muJjKTISPrjDwoJkbkzfD4dOUKRkRQfT7a21Ls3jR9f+xG6YiCgAaAepK5cKTAyMj9woKqlxNU1IzCwxrlsPT29oKAgIqLwcCmHw2K+vurOzkaJiZSQQKdO0Zs3tW798dat2V27unt4GBsb/9caEvJWOlc5eJBGjqThw6XMys0lf/+3Lt/+8UcKCaFjx8jcvNYy6h0CGhqFS5cumdV2JBgVFTVI7reiQQ6Rhkb6kiWZs2bpxMaq5+aW29qW29nJe9D46NEyA3rbNnJwqHq1cNYstcTEkUlJfaKjZa3M9X//ExGVGhmRnR1ZW5OlJVlZkbQj+n8cPiw9oOfMkfLlmshImjVL3klzhUFAQ8M3e/bsBw8e1NrNyspqsOxvRQMTAgODYtnXUL9l3jy6eZOOHq3ZvmdP9XQmIr6GxpDAQNuUFJo7V9bK0pcte/nkiWlZWXsTE0pPz4+IMCwtVZN925D8336727mzoEmTIZMnk5kZmZuTuTnx+RQWJn2BU6fo5Utq04bRrtUfBDQ0fF26dOnCMDVAaTQ06NdfafBgOnqULlwgb29q0YIWLqRqX3esrqhbt5JOnfT++ktyVmHfvllTp946f97Y2Lh9QAAR5b9+nV1R0WLsWJ1Hj6SuzaiycuDz57zycuZfRo86eDCjd28nJydrJT4qDAENACzh8WjaNJo2jUlfkZZWxpIldlOnSs5K/fTTGi0tW7YkInJ1JRkBzVu0iL75hoqLKSuL3rz59YcfMp498ygt7fbsmawC/vrzz7MxMQ4ODl9++SWTgusFOwF98eLFoKCg1NTUrl27hoSENJPzxAfOKCoqKigokN+ntLS0ev+///671tVWXwRUV7NmzXbs2CHlUa1vq6iomDFjhnJKaniK3dxenD1rsWePZlqa/oMHxV27ljk5vZk1SyDrazLLltGPP0qfFRhIRKSvT/r61KpVXNu2TXv2NG3WjPz8ZG3dY9o0YUnJixcv3n9HmGMhoPPz8ydNmnTmzJmuXbuuWbMmMDDwqOR5KI7x8vI6fvw4k54DBw4kIh0dnSFDhjx8+LDW/uPHj3/f4oADZs6cOXPmTLarYIGmpuann36qrq4uv1tCQsKQIUPef3PltrbJX3/NtHeHDnT6tJQPA69eJRsbKSu3syvy8DC4e1dyVnHXrmWOjrKOxxWHhYAODw93c3Pz9PQkoqCgIEtLy/Lycm1t7fdZ571797Rqu1Dx9evX77z+SZMmTZo0iXl/dXX1uv4fZGdnt3///iNHjsjvVllZOWXKlDqtGUBxjhw5UlJSUms3dXV1IyOjd94Kw/9f37rGTszPjzIz6fhxiooiDQ1ycSF/f5JdSernnztI+0OS+r//sfLdFhYCOjY21uHfT2mNjY1NTEwSExMd3v7clohevHjx+PHjGo3Xrl3Lz8+3sLCo3ujo6BjO4OHEampqiYmJR48eTUlJSUpKYrJIZmam0o7utbW1582bx7CzuKrU1NSNGzfW2vnFixenTp3S19e/f/9+bm6ubm33ysnOzn716pV4EykpKbNmzap1EwkJCW5ubtHR0VFRUXFxcbW+twKBICMjQ7yJtLS033//vda/r7GxsTo6OkVFRampqcnJyQodvvLy8pcvXzJ5b5OTk9/5JyQnJ4fJXrx58+b06dPi4Xv+/HlSUpL8/oWFhWlpaeKqMjMzmWzi9evX4eHhMTExDIdPKBRWDZ8SiESi4ODgWrvxeLzRo0dLr6pJE+rT559pabfvePr0qbGxcWFhIRGprV1rFBGh9fq19suX5fb2FdbWBX37Cl+8oBcv4uPj09PTa2yCx+MNHz5cUzFfImchoEtKSqpnhJ6entS/wKWlpbm5uTUazczM9PX1a7Q7OTk5OTkx3Hpubm5BQQGfz6/1DzIR8fl8yRq4Y8iQIUxOYfft27e8vLyioqKkpKS0tLTWHS8qKqqoqBDveEBAQHl5OZNixOPC5/MfPnwYGxsrv7NIJDI3Nxdvgs/nh4eHM/kfuU2bNuLhq6ioUOjwiUSiAQMGCBg8EOSDDz54558QgUDAcC/y8vIqKiosLS1fvnyZkJBQ6yLOzs5V7y2TTYjfz9zc3KKiovLy8loXEQqFlZWVSvvVqNNFOO9WVWlpqYaGRtWO5/Xu/dZsoZAKCohI/BtUYxNqampChT0GjIWA1tfXT632vIaioiIDAwPJbi4uLi7S7qf1/qKjo589e+bv719rz8jIyNmzZyuihnpR19o0NTVTU1OHS70+v5rk5OSioqJ33nEmBzvV9erVKyUlhUlPDw8PAwODuLi4kJCQZ7I/ba9iZGT0znsREBDwbgsyt2/fPiY/hDdv3pwyZYqp3NvFyfLLL78w2URcXJy/v7+Li8vVq1f/+OMPJgFtZWXF5V+NusrMzDQwMOjatav8bk+ePElJSVHmjrMQ0E5OTpf+/S8jJSWluLjY1tZW+WUAR7Rr165du3bM+9vb21+5coXJMYuJ1HsQg2yenp5r165l0rMhpTMRubm5HThw4NatW7X2nDhxohLqqcJCQPv4+AQEBFy9erV3797r1q3z9/fX0FB2GRkZGXelfVZbQ0VFhRKKgbpqGH/RDQwMlsh6Nms1OTk57/wROp/PZ/JznpaWJp7Q1dUdIHFD0cbA19fX19eX7SqkYCGgDQ0NQ0NDAwMD09PTe/ToESLntlKKYWNj88EHH2RnZ9fak/mndgB1dVXxD1SdP38+k4uXBg0a1DD+5jU87HxRxcfHh8k5RAXR19f/535aAA3auHHj2C4B3gu+6t2IWFpa7t69OyIiQn43Pp/ft29fpVQEAPLwRCIR2zUAAIAUeGgsAABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEc1wIC+ePGis7OzmZmZr69venp6rf0rKyuXL1+upqaWlZXFcBNnzpxp166diYmJt7d3bGxsrf2PHDni4OBgYmLSr1+/uLg4hlshooiICB6PFx0dXWvP8vJyHo+n86+xY8fWusjr16/79u1rYGDQqVOnx48f19o/JCREpxoej5ednS1/kbCwsA4dOtjb2/v4+MTHx9e6CSI6ceJEu3btjIyMhg8fnp+fL6ub1FGTM/RS+8sfeqlz5Q+91EXkjL6cAmQNveQi8ode6ibkD73kIvKHXuom5A+91EXkDL3Ut13+b7rUReSPuNRF6vrLXs9EDUteXp65ufnt27f5fP7KlSv9/f1rXWT48OGrV69WV1fPzMxksomkpCRjY+ObN28KBIKVK1f27dtXfv+YmBgzM7MnT55UVlYGBQX5+Pgw2hORqKyszNXV1dLS8vnz57V2TktLMzc3Z7hmsV69em3cuLGiouLHH3+cMmVKnZa9fPlyv3795PdJTk42MTGJj48XiUSbN2/u379/ratNSEgwNjZ++PBhWVnZ+PHj586dK6un5KjJH3qpoyx/6CXn1jr0kovIH31ZBcgZeslF5A+91E3IH3r5b4vk0Ev2r3XoJReRM/RS33b5wy1rpOTsmtRF6vrLXu8aWkCHhYUNGjRIPJ2Xl6etrV1WViZ/kUePHolEojoFdFhYmHj64cOHVlZW8vu/evXqwoUL4uk7d+60bNmSyVZEItHq1as///zzDh06MAno6OjoNm3aMFyzSCSKj4+3trYWCATMF6nC5/OdnZ2fPHkiv1tkZGTHjh3F01FRUZaWlrWu+aeffvrggw/E08+ePTMzM5PVU3LU5A+91FGWP/SSc2sdeslF5I++rALkDL3kIvKHXrJ/rUMv522ROvSS/WsdeslF5Ay91Ldd/nDLGik5uyZ1kbr+ste7hnaKIzY21sHBQTxtbGxsYmKSmJgofxFXV9c6bcLa2nrMmDHi6cuXL3t5ecnvb2tr6+vrS0QFBQW7d+/28/NjspXY2Njjx49/+umnDKvKy8srKSnp169f06ZNfXx8av1f7K+//mrbtm1AQICtrW3//v2ZnEWpcvjw4Q4dOjg7O8vv1qlTp8zMzIcPH4pEolOnTvn4+DDfBBGZmJhkZ2fn5eVJnSs5avKHXuooyx96ybm1Dr3kIvJHX2oB8odechH5Qy/Zv9ahl/O2SB16yf61Dr38d77G0Et92+UPt6yRkrNdqYvU9Ze9/in5D4KiffbZZ8uXL696aWdnJ/6bWSvmR9BVLl68aGdnl5SUxKTzxx9/TES9evXKzs5m0r9fv36XL18WiUQMj6CfPXs2e/bs58+fl5aWrlixwtnZWX7/H3/8UVdX98qVK0KhcPPmzbX2r659+/bi371aHTp0SENDw9TU1Nra+uXLl7X2T0xMNDIyun//Pp/PX7p0qZqaWlpampz+1UeNydBLHWX5Qy91rvyhl1xE/ujX6M9k6KsvwmToq/dnOPRSd1zO0Nfoz2Toqy/CZOirv+0Mf9OljpT8EZe6SJ1+2etRQwvodevWzZ8/v+qlhYXFixcvmCxY14A+fPiwo6Mjw5WLFRcXb9q0ycXFRSgUyu8ZEhIyYcIE8TTDgK6uoqJCS0srJSVFIaNmgQAACcFJREFUTp8TJ064urqKpysrK7W0tLKyspis/P79++3atWPS88mTJ7a2tgkJCSKR6OjRo46OjkzOqISFhXXs2NHJyemHH35QU1MrKSmR07n6qDEZ+noJ6FqHXuoK5Yx+9f4Mh15WzbKGvnp/hkMvuQn5Q1+9P8Ohr7EJ+UNf421nMtyyRkrOiEtd5B1+2etLQwvo48eP9+7dWzydnJysp6fH5/OZLFingD516pSzs7P8g7sqjx8/vnLliniaz+fXelQoEolGjBhhbm5uaWlpaWmpoaFhbm5+9uxZ+YukpqY+e/ZMPF1WVlbr7jx69KhVq1bi6crKSk1Nzby8PCa78/nnnwcFBTHpuXXr1vHjx1e91NLSqtMByB9//OHk5CS/T/XdZDL07x/QTIa++iJMRr96f4ZDX30RJkNfvZHh0EuuR/7QV+/PcOhlvfOSQy/5ttc63HJGStZ2pS5Sp1/2etfQArqgoMDc3PzKlSt8Pn/u3LmTJ09muCDzgM7JybG2tn716hXDNV+4cMHKyiouLk4kEoWEhFhaWtZ6BF0dwyPo8+fPt2zZ8tWrV5WVlatWrerevXuti7i4uBw4cEAoFH777bfu7u4M6xk2bNiPP/7IpGd4eLitra346Cw8PNzMzKzWP5ZpaWn29vbx8fFFRUUDBw5cv369/P7VR43J0L9nQDMc+uqLMBl9WQUwPIJmMvQ1NsFk6CWrkj/01fszHPoa16LIGnqpb7v84ZY/UlLfcKmL1PWXvd41tIAWiUTh4eHt27dv0qTJ0KFDa/23PSsrS1tbW1tbm4jEE+np6fIXOXDgAI/H066m1q188803rVq1MjEx6dKlS2RkZJ12h/kpjuDgYCsrKzMzs0GDBjH5kYqLi+vSpYuJiUnPnj2jo6MZ1tO5c+eqyxJqtW7dOnt7e3t7e3d3d4Y7vnXrVnNzc1NT0zlz5sgKdFmjJmvopfaXP/RS58ofelkrlDX6tf7sSQ69rEVkDb2s/nKGXk5VUodeVn85Qy9rEVlDL+ttl/ObLnUR+W+41EXe4Ze9fvFEIpGyP5cEAAAGGtpldgAADQYCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENDwjiorK3k83qJFi6paTp06NWLEiHdb2+PHj+3t7euptP/MnTvXzMzs1KlTVS1CoXDz5s1OTk7a2trNmjWbNWvWmzdv6n27APUCAQ3vTldX9+jRo3///TfbhcgUFhZ2586d6n82Fi9evHv37u3bt2dmZl69ejUnJ6dfv36VlZUsFgkgCwIa3p1IJFq9evWSJUtqtN+9e9fJyanG9OPHj11dXZcvX+7t7d2hQ4erV6+OHj26U6dOgYGB4p48Hm/ZsmUtW7Zs27bttWvXxI3nzp1zcXFp06bNgAEDMjMziejRo0ddunSZNGmSj49P9Y0ePXq0Y8eO7dq18/b2jo6OJqKRI0fm5eX5+fmdO3dO3CcpKWnnzp0nT54cMGCAkZFR+/btw8LC+Hx+SEgIER06dKh169bW1tYTJ04sLy+XbJG6X9XrqVGbZPHiN2HlypUDBw50cnK6dOmSeG2Sm5ZctqKiYtKkSW3atLGzs5swYUJpaWl9jSNwlzIfIQ4NCZ/P19DQEAgEnTp1On36tEgkOnny5PDhw0Ui0Z07dxwdHcXdqqajoqLU1NRu3bolEonGjx/foUOHsrKy0tJSAwODtLS0R48eaWhoHDhwQCQShYSEtG3bViQSJScnN2nS5K+//hKJRJs2bRoxYoR4Pfr6+r/88kv1YhITE42MjGJiYkQi0fbt27t37y5uNzY2TkpKqup2+PDhDh061NiRVatWjR49+tWrV2ZmZvHx8ZWVlUOHDg0ODpZskbVfVfVUn5ZVvJqaWnh4uEgkOnr0qKenp0gkktyQ1GWPHj3q4+MjFAoFAsGyZctu3rxZPwMJHIaAhnfE5/PV1dVFItG1a9fs7e3Ly8trDWgzMzNx46pVq+bOnSuetre3f/jw4aNHj/T19QUCgUgkKi8v5/F4mZmZu3fv9vHxEXcrLCzU1NSsqKiIiorS1dUV96yyb9++YcOGiadLS0vV1NTy8/NFEgG9Y8eO/v3719iRHTt2eHt779mzx8/PT9xSXFxcVlYm2SJrv6rqqT4tq3hjY2Nx4+PHj1u2bCkSiSQ3JHXZGzduWFlZnTt3rrS0tA7jBKpMg+0jeFB53t7eLi4uW7dudXR0lN/TwMBAPKGurq6np1c1LRAINDQ0zMzM1NTUiEhLS0tPTy83NzcvL+/27du2trbinvr6+llZWURkamoq7lklMzPT3NxcPK2jo6Onp/fmzRsjI6MaBdjY2KSmptZoTEtLa9q0aW5urqmpqbhFXJtkiyzV66mallV89TdBIBBI3ZDUZb28vLZu3bpx48bx48ePHj36+++/19fXl1MVNAA4Bw31YNOmTZs3b666HEJdXV0oFIqnCwoKGK4kNzdXJBIRUUVFRXFxsZmZWYsWLQYMGJDwr9zc3ObNmxMRj8ersaylpWV2drZ4uqysrKSkxNLSUnITXbt2TUxMfPz4cVWLUCg8duzYBx98YG5uLj7VK64kKSlJskXWflWvp2paVvGSJDcka9kxY8ZERETEx8e/fv16x44dtb+noOIQ0FAP7OzsZs6c+c0334hfNm/ePDU1tbi4mIhOnDjBcCXFxcXHjx8nol9//bVDhw5NmjTx8fG5c+dObGwsEd2/f7/q40RJPj4+169fj4uLI6Jdu3b17NnT0NBQspulpeXHH3/s7+8fERFRVlaWkpIyfvx4PT29iRMnDho06NatW8+ePRMIBLNnzz58+LBkS532i3nxkhuSuux33323Zs0akUjUpEkTGxsbyb9S0PDgFAfUj88+++zgwYPiaWtr61mzZnl6etrY2PTv37/qkgw5+Hx+hw4d7t69u2rVKiLav38/EVlaWh44cMDf37+0tNTY2Pjbb7+Vtbi1tfWBAwdGjRpVUVFhY2MjvipDqjVr1hgbGwcEBCQkJJiYmPj5+V26dElNTc3Kymrv3r2+vr6VlZUDBgxYsmSJtra2ZAvz/WJevNRNSy47bty4adOmtWrVSl1dvWvXrnPmzKn1XQVVxxP/UwkAAFyDUxwAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgqP8D88uPJ8fi5F8AAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3dd0AT5/8H8E9C2DOCIkNFRQQRxS0VNyBaVwU6XK2tin5bHLVaa221Vn/ar1qrtbaOKrZ1FPdCxbqw1fZrW6mKIOJA9p4JBAj5/RGLSNahJHdJ3q+/Lk+eu/skD7w5nlzueDKZjAAAgHv4bBcAAADKIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwlIDtAvRWWhplZJCXF7Vsyah/WRnduUMCAfn4kLW1losDAEOAI+im++EH4vHIw4MCA6lVKxo0iK5cUdc/N5dee43s7SkggPr0IRsbmj6dSkt1VS4A6CueTCZjuwa98vXXNGeOkvZLl2jwYCXt5eU0ciT99lvj9mHD6MwZMjVt/goBwFDgCLopCguVpzMRffKJ8vavvlKSzkR04QJFRzdXXQBgkDAH3RRxcSqfunKFLl6kdu2eaRSL6dgxlascP04zZjRbbQBgcDDFwcjGjRt//PHHybm58zIymmub/wiFH/TqZWpqGhsb21zbBABDgikORu7cufPZZ59FzJyppk/OBx8cmzTpfGQkxcRIfvxxmY9P+vr1kg4dVPXv2KXLhg0bCgoKtFAvABgCBHQTiPr1U/Ns4ZQpyX5+D3r1ooiI2ldeudSyZWlISHlgoKr+FWq3BgCAgG6CalfXojfeUPpU5sqVMr6SN7Ng2jSl/cX+/kWvvdacxQGAwUFAN032okWFEyc2asxasqR43Dil/WudnO4dOSLq1euZRqEw44sv6qystFUlABgEnMXRNDKBoHDyZMd9+0rGjq309q728BB36ya1s1OzisTT8+GuXZbJyeYpKWRqavnPP45795oUF5Orq87KBgB9hIBuMoeTJ4nHy50zp6ZVK6br8HiVPj6VPj5EVDZokMPJky137Hi8YYMWqwQA/Ycpjiazj42t6NOnCen8rDobm8KJE+3On7e4d695CwMAA4OAbhrLW7fMHz0qHT36RTZSOHWq1Nq65fbtzVUVABgkBHTTOJw6JTM3Lxs+/EU2IrW1LXrtNbuzZ80fPWqmugDAACGgm4BXV2d/9mz54MFSW9sX3FTBW2/JLCxwEA0AaiCgm8D62jVBQUHJyy+/+KakDg5F4eH2J0+2qap68a0BgEFCQDeBw6lTUjs7NV8ObJKCadNkpqZTc3KaZWsAYHgQ0EyZSCR2Fy6UhoTIzMyaZYO1Tk7FEya8XFBAaWnNskEAMDAIaKac//iDLxK94PkbjRS8/baMiNaubcZtAoDBQEAz5XbpUk3r1qIePZpxmzWtW59ycqIdOygzsxk3CwCGAQHNiLVE4vT33yWjRpGyKyK9iF2urlRXR+vXN+9mAcAAIKAZ6XP/Pr+2tnnnN+Syzcxo4kTaupXy8pp94wCg1xDQjPS9d6+ibduqTp20svWlS6m6mr78UisbBwC9hYBmIC2tQ25uxrBh2tq+pydFRNDmzYS7qwBAAwhoBn76iYiyBg7U4i4++YQqK2nTJi3uAgD0DQKagX377rm4VDo7a3EXPj70yiu0aROVlGhxLwCgVxDQmty4QYmJ/9PS7HNDn3xCZWW0ebPWdwQAegIBrcmePWRm9reHh9Z31L07jR5NGzZQebnW9wUA+gABrVZdHe3fT6NGiSwsdLG75cupuJi++04X+wIAzkNAq3XxImVm0qRJOtpdz54UHEzr15NYrKM9AgCHIaDV2rOH7OyoOa4vytTy5ZSbS7hONAAgoNWpqqIjRyg8nCwtdbfTgAAaMoS++IIqK3W3UwDgJAS0aidPUkmJ7uY36n3yCWVnU3S0rvcLABwjYLsADtuzh1xdafBgXe932DAKDKT/+z/q0YMyM8nDg3x9STefUgIAl+AIWoXiYjp9mt54g0xMWNh7796UkUEBARQeTr17k6UlbdvGQhkAwCoEtAoHDpBEwsL8BhHt2EFffdW4MTISp98BGBsEtAp79pC3NzXr5fkZqaqiGTOUPzV7Nk6/AzAqCGhl0tPp119p8mQWdn31qrpnf/tNV3UAAPsQ0Mrs20cyGU2cqMt9fv/994GBgf83b56aPmvef3/gwIHfYa4DwDjoKKAvXbrE4/GSk5PlD8+cOePn5+fo6BgaGpqTk6ObGprgp5/opZeofXtd7jMrK2vatGlvzJ2rps/rc+ZMnz49OztbZ1UBAIt0EdASiWT+/PnO/16us7S0dMqUKdu2bcvNze3du3dUVJQOamiCO3fo1i12Ph4kqvT3V/OsWPdz4gDAHl0E9OrVq8eOHevk5CR/GBcX16tXr4CAAIFAsHDhwhMnTkgkEh2UwdQPP5BAQGFhrOy8ztw8a9kypU9lLV1ah7OhAYyJ1gM6JSXl0KFDH330UcMWLy8v+bK9vb2Dg0NaWpq2y2BKJqOff6bQUGrViq0SisLDMxUyOuuTT4pee42VegCALVr/JuHs2bM3bNhg0eDQTywWWza4uoWVlZVY2dljly9fPn36dKPG4uLi4cOHv/rqq1qqluLj6dEjWr1aW9tnpjg8vHTMGPN798wyMx1iY63/+KN43Dh2SwIA3dNuQO/evdvFxSUoKKhho7W1dVZWVv3DiooKGxsbxXV9fX0V28+dO5eRkaGNUp/Ys4esrWnMGC3ugpk6c/PKrl0ru3aV2traXrhgc/Vq+dChbBcFADql3YA+evTor7/+2rp1ayIqLCwcOHDgrl27vL29z549K++QmZkpEok8lN2vxMnJqX7aut79+/cLtHfr6+pqOnyYJkwga2tt7aLpRP361QqF9nFxCGgAY6PdgD5y5Ej9cteuXQ8ePOjt7V1eXh4ZGXnhwoVBgwatWrUqPDxcIGD1mk03blBMDN25QyUlVFhIo0ezWYwCmYlJ+fDh9qdP8yWSOnNztssBAN1h4Ysqtra2e/fujYqKcnZ2Tk9P//LLL3Vfw1MrVlDPnrRmDR0/TvHxRESvvfZkgTNKQ0L4IpGN+i8ZAoDB0d2h6+3bt+uXg4ODExMTdbZrlQ4eJKXntA0eTAUF5Oio84KUE/XtW9uihd3Zs2WY5QAwJsb9VW8135netUuHdWggMzEpGzbM7tIlPqdOGAcALTPGgJZKpRcuXPjl3Dk6f15Vn9zY2F9++SUhIUGXhalRNmIEXySywcWSAIyJMQZ0amrq/Pnzz8TGqumTm54eGxsbGRmps6rUE/XpUysU2v179gsAGANjDGiZTObt7f3WO++IevZU1cdl6NC3337byspKl4WpITMxKQsKsrt0yaS6mu1aAEBHjDGg65UPGaLqqRIOfFelkbIRI/hisdutW2wXAgA6YtQBXTh5stKMzl6yRKLba40yIerTp7ZFC4/r19kuBAB0xKgDWmZq+njjxuwlS6S2tvKWsuHDH+7aVfjGG+wWppSMzy8LCmpz44ZZTQ3btQCALhh1QBORjM8vfP11ksmK3njj9q1bj7/6StS7N9tFqVQaEiKQSDqmprJdCADogrEHNBGZZWSYVFRUdu7MdiGaifv0qbSz871zh+1CAEAXENBkkZxMRFXe3mwXopmMz3/cu7d3SgqJRGzXAgBah4Amy7t3ZSYmVZ6ebBfCyKM+fUxrakjhStkAYHgQ0GSRlCTp0EGmJxeKy/X2LrexoQMH2C4EALQOAU0Wycl6Mb8hJ+Pxkry96eRJzHIAGDxjD2hBcbFpXp4eBTQRJXbpQmIxqf2qOgAYAGMPaIukJCKq1KuATmvXjlxdMcsBYPCMPqDlp3Dowzl29WQ8Hr3yCp06hVkOAMNm7AFtefdujaur1N6e7UKaKCKCxGI6dYrtOgBAi4w9oC2SkvRrfuOJgQMxywFg8Iw6oPlVVWZpafo1v/EEn08TJlBsLFVUsF0KAGiLUQe0RUoKr66uyseH7UKeC2Y5AAydcQe0/BQOfTyCJqLAQHJzwywHgAEz7oC+e1dqZ1fj4sJ2Ic9FPstx6hSVl7NdCgBohVEHtGVycpWPD/F4bBfyvCIiqKoKsxwAhsp4A5pXV2d+756+zm/IDRiAWQ4AA2a8AW328CG/qkq/vuTdGJ9PYWEUG0tlZWyXAgDNz3gD2lL+HUI9PYWjHmY5AAyX8Qa0RXKyzMxM4uHBdiEvZsAAatsWsxwABsmoA7qqUyeZQMB2IS+Gx6MJE+j0acxyABge4w1oy7t39X5+Q04+y3HyJNt1AEAzM9KAdhSLTYqL9fsUjnoBAZjlADBIRhrQHiUlpCc3itWMx6OwMDp9mkpL2S4FAJqTkQZ0+5IS4vOrvLzYLqSZRESQRIJZDgADY6QB3a64WNK2bZ2VFduFNJP+/aldO8xyABgYIw1oj5ISA5nfkJPPcpw5g1kOAEOi5yeZPReT8vKWIlGeIQU0EXXuTBIJBQRQy5bk7Exvvkkvv8x2TQDwQozxCNo8KYmnv1cZVerLLykykogoKYni4+nAARo9mubNY7ssAHghxhjQT24UazBH0H/+SQsWKGnfuJEOH9Z5NQDQbIwxoM2TkootLGqdnNgupJls367yqZ07dVgHADQzY5yDNk9OTnZwMIBX/vPPP5eWlo6+fNlVVY9Tp7Zt2yYQCN566y0+3xj/GAPoNeP7pa2uNr9//5FQyHYdzWD58uUVFRW1UqmaPhUVFRs2bCjDlToA9JDxBfTt27yamocODmzX0QxsbW1DQkKs+/RR1aFs6NCQkJBWrVrpsioAaC7GF9AJCUSUZhBH0HJFr76q8qnXX9dlJQDQvIwxoOusrHKsrdmuo9lIPD0zVq5UbM+dM6fipZd0Xw8ANBcD+KisiW7ckPj4yPT3RrHKlIwbV9m1q/DQIYsHD0zKyixv3cqbPTt/xgy26wKAF6L1I+j9+/d7eXk5ODgMGzYsNTVV3njmzBk/Pz9HR8fQ0NCcnBxt1/CUTEY3bxrOGdANSDp2zFm06NF3393fs6fG1dUyKYntigDgRWk3oFNSUt57771Dhw4VFhb27t37P//5DxGVlpZOmTJl27Ztubm5vXv3joqK0moNz7h/n8rKJIZxnX5VeLzSoCCbq1dNKirYLgUAXoh2A9rMzOynn37y8/MzMTGZMGHC3bt3iSguLq5Xr14BAQECgWDhwoUnTpyQSCRaLeOpGzfIAG4Uq0nZiBG86mrbS5fYLgQAXoh2A9rDwyM0NJSIysrKtm7dOnbsWCJKSUnx+vdCzPb29g4ODmlpaVot46mEBBIIJB076mh3LBH7+dW4utrHxbFdCAC8EF18SLhw4cJ169YNHDjw6NGjRCQWiy0tLeuftbKyEovFimsdOnTou+++a9SYm5sbFBT0/KUkJFCXLjJz8+ffgl7g8UqDghz37zepqJDa2LBdDQA8J10E9Nq1az/77LNvv/126NChCQkJ1tbWWVlZ9c9WVFTYKAuRsLCwsLCwRo0xMTEFBQXPX8qNGxQc/Pyr64+yESOcfvjB9tKlktGj2a4FAJ6Tdqc4/vnnnwsXLhCRlZXV3Llzb9++nZub6+3tfevWLXmHzMxMkUjk4eGh1TKeyMuj7Gzy99fFvtiGWQ4AA6DdgM7Ozp46der9+/eJaM+ePS1btnR2dg4ODr5z586FCxdqa2tXrVoVHh4uEOjkdOwbN4iIevTQxb5Yx+OVDh9u89tvOJcDQH9pN6BDQ0Pnzp07fPhwoVC4adOmmJgYHo9na2u7d+/eqKgoZ2fn9PT0L7/8Uqs1PJWQQDwede+uo92x7cm5HJcvs10IADwnrR+6Lly4cOHChY0ag4ODExMTtb3rxhISqF07EgopN1fXu2aDuFu3GldXu7NnybC+NglgPIzpWhw3bhjL/IYcj1c6fLjtb79Z1dayXQoAPA+jCWiRiO7dM5JPCOvJZzkCCgvZLgQAnofRBPTNm1RXZ2wBLe7WrcbFZXBeHtuFAMDzMJqANqpTOOrxeGVBQX0LC3nl5WyXAgBNZjQBnZBAjo7Upg3bdehaaUiIaV2d6dmzbBcCAE1mTAFtZPMbcuLu3fMsLEyPHmW7EABoMuMIaKmUEhONbn5Djse70rKl4JdfCPeNBdA3xhHQSUkkFhvnETQRXW7ViieR0KlTbBcCAE1jHAGdkEBERhvQd+zt69zd6cABtgsBgKYxmoC2sKDOndmugx0yoprRo+n0acxyAOgX4wjoGzfIz490c0kmTqoZP56qqjDLAaBfjCOgb9400k8I/1Xbty+1aYNZDgD9YgQB/fgxFRQY7QT0EzwevfIKZjkA9IsRBLRxf0L4VEQEVVVRbCzbdQAAU0YQ0DduEJ9Pfn5s18G2AQMwywGgX4wgoBMSyMuLcO9U+SxHbCzhuhwAesI4AhrzG3LyWQ6cywGgJww9oEtKKC0NAf3ESy+RmxtmOQD0haEH9I0bJJMZ+Tl2T/H5NGECnT5NuJMsgD4w9ICWn8JhNDeK1SwigiorMcsBoBeMIKBdXcnZme06OGPAAMxyAOgLQw9oY7tRrEaY5QDQHwYd0BIJJSfjE8LGIiJILMY3VgC4z6AD+vZtqqlBQDc2YAC5umKWA4D7DDqg5TeKRUA3Ip/liI3FLAcAxxl0QCckkJ0ddejAdh3cg1kOAH1g6AHdvTvxDfo1Pp/AQMxyAHCf4YZXXR3dvIn5DeXqZzlEIrZLAQCVDDeg79+n8nIEtEqY5QDgPMMNaPknhDgJWpX6WY7cXJJK2a4GAJQwxICurqZ//qGjR0kgoC5d2K6Gq9LSyMqKDhyg1q1JIKDQULp0ie2aAOAZhhXQUimtXEnm5uTvT/v2UW0thYTQ9etsl8U9SUnUoQOlpj5tOXuWhg6lmBj2agKAxpoW0I8ePbp+/fr//ve/hw8faqmgFzJ/Pn3yyTMt8fHUty/98w9LBXHVggXK2197DTctBOAORgF9//792bNnOzs79+jRY9KkSZMnT+7Zs6ezs/OsWbPu37+v7RKZun2bvv5a+VMff6zbUrgtO5tOn1b57NmzOiwFANTRHNDr168PDg5u165dfHx8cXFxSkpKSkpKcXHxlStXPDw8goODv/zySx0UqtnJkyqfOnUK35ojor/++isyMvKL2bPV9Dm8bl1kZOT+/ft1VhUAqKI5oB8+fHj79u3Fixd37ty5YbuXl9fixYtv37794MEDrZXHyIcffhgUFHTk22/V9HkzNHTgwIF//fWXzqrioJs3b9ra2gaOGKGmj3fv3gEBAefOndNZVQCgiuaA3rx5s5WVlb+//5o1axSz2MrKavPmzdqpjanbt2+vWbNmwNixavosXru2T58+eXl5OquKm4RCoUNAgJoOZkOHtm7dWmf1AIAaTD8kXLly5b179/r379+3b9/169c/fvxYq2U9h3LVuVPRr5/U1laXxXCZTCDIUfEhYVlQUJW3t47rAQBVmAb06NGjv//+++zs7HXr1qWlpfXr12/AgAFbtmwp48yH/pXdupWOGqX0qfxZs3RcDMcVvPlmwVtvNWosHzgw8/PP2SgHAJRr2ml2FRUV9+/fT01NFYvFjo6ON2/e7Nat22k1pwToVuaKFUXh4Y0a0zZvFvXuzUo93MXj5SxYcO/YseyPPiqcOLF4/HgiKp4wQWpjw3ZlAPCUgGG/gwcP7t2798yZM/369Zs4ceKePXuEQiER3blzZ+TIkWlpadoskqk6c/OsZcvyZ8+2uHNHUFoq8fCo9PGRmZmxXRdHSTp0kMivxVpXZ/O//7U4eLAsKIjtogDgKaYB/cUXX7zxxhubN292dXVt2N6lS5eJEydqobDnV9OqVU2rVmxXoVf4/OJXXmm1ZYtZenp1mzZsVwMATzCd4rhy5UpVVZWjoyMRpaenf/HFFxKJRP7U6tWrtVUd6EpRWJiMzxcePcp2IQDwFNOAnj59+u+//y6VSonIxsbmzz//nDlzJpMVjx8/7uPj4+DgMGTIkJSUFHnjmTNn/Pz8HB0dQ0NDc3Jynq90aEa1LVtWDBokPHyYV1vLdi0A8ATTgL506dKRI0esrKyISCgU7tu3j8l3GTIyMqZOnbpjx46ioqLAwMBZs2YRUWlp6ZQpU7Zt25abm9u7d++oqKgXeQHQXIrCwgQFBbaXL7NdCAA8wTSgTUxMMjMz6x8mJyebmJgwWXH79u0DBgzg8/lhYWHyI+i4uLhevXoFBAQIBIKFCxeeOHGifrYEWFQ+cGCNq6vw4EG2CwGAJ5h+SLhq1apevXoNGDDAwcEhPz//119/3bFjh8a13N3dIyIi5Mvnzp0LDAwkopSUFC8vL3mjvb29g4NDWlpafQuwhs8vHjeu1datlrm5bJcCAETMA3ry5MkDBgw4d+5cSUlJy5Ytt2/f3uh0DvXOnj373XffxcfHE5FYLLa0tKx/ysrKSiwWK67yxx9/XFK4hPzNmzc9PT2Z7xeapDgsrOXWre5xcWRhwXYtAMA4oImoffv2DD8YbGTv3r0rVqyIi4tzd3cnImtr66ysrPpnKyoqbJR9P8LV1bVXr16NGisqKqytrZ+jBmCixtm5YtAg9zNnTNRe2AQAdINpQMfGxn7++edZWVnSBvevy8jI0LjisWPH1qxZc+nSpfpL8Hh7e5/996LDmZmZIpHIw8NDccU2bdq0UTgnt6ioqKCggGHN8ByKwsLaXbrUjcHIAoC2MQ3oBQsWfP755506dWL42aBccXHxe++9d+XKlYYXSAsODo6MjLxw4cKgQYNWrVoVHh4uEDThQB60qmLQoKpWrQbdvct2IQDAOKDbtGkTrnCZC42OHj2amZnp3eACaZmZmY6Ojnv37o2KisrJyXnppZeio6ObulnQHhmfnxEc7Ld3L6WlUbt2bJcDYNSYnmbXt2/f8+fPN3Xr06ZNq6urq2pA/l3E4ODgxMTEwsLCEydOyFuAOzJCQ2VEtHMn24UAGDumAX3y5MmQkBAHBwf3BrRaGbClysnptrs7ff894VuFAKxiOsURHR2NmWLjEd+5c7dffqHYWMLpHADsYZq5/v7+RFRTU1NQUODi4qLNkoB9/7i7k5sbbd+OgAZgEdMpjry8vFdeecXa2rpLly5EFBUVdfXqVW0WBmyq4/Pp7bcpNpa4caVvAOPENKBff/31wMDAoqIie3t7Ipo6dercuXO1WRiwbcYM4vFo1y626wAwXkwDOi0tbcGCBfVf+evTp09JSYnWqgIOaNOGQkLo+++pwVeTAECXmAa0jY3NzZs36x8mJydb4HINBm/mTMrIIM7ccxLA2DAN6FWrVg0dOjQkJKSgoGD8+PGDBg1as2aNVisD9o0e/eSjQgBgA9OzOEaPHn3z5s1Tp04FBQW5uLhs3brV2dlZq5UB+wQCmjaNVq+mx4+pbVu2qwEwOk04tdnNze35rmYHemzGDFq9mqKj6dNP2S4FwOgwDeiG19OQk8lkd3FJHYPXti0FB9OOHfTxx9SU62QBwItrwjcJ65fLysoOHjyoGNlgmGbMoLAwOnOGXn6Z7VIAjAvTgO7fv3/DhyEhIS+//PL777+vhZKAY8aOJVdX2r4dAQ2gY0zP4mhEJBKl4TtmRkIgoLfeolOnqMFdgwFAB5gGtHcDXl5ezs7OI0aM0GplwCHTp1NdHS5ACqBjzzMHbWJi4uLigsuNGpH27SkoiLZvpyVL8FEhgM485xw0GJ0ZMygiguLiaORItksBMBZMA7pdu3ampqaqnk1NTW2meoCrxo8nFxfavh0BDaAzTAN6zpw5f/7559SpU1u2bJmbm7tt27aAgIDQ0FCtFgccIhDQlCm0fj0tWEAyGXXpQmPGEL5NCqBNTAN69+7dCQkJfP6TDxVDQkK6d+++ePFirRUGHJOSQufPk1RKX375tHHXLnrrLdZKAjB0TM/iKC4uzs/Pr3+Yl5dXXFysnZKAeyorado0+uuvxu3TptGFC2wUBGAUmB5BR0ZG+vj4DB482N7evqSk5PLly7hgvxHZs4dU3UBn7VoaNky31QAYC6YBvXTp0gkTJsTHxxcVFQmFwhUrVnTr1k2rlQEXFBQUlJWVtTx92lZVjzNnHty/b2pm1qZNG10WBmAMmnA1uy5dunTq1Ak3jTUqY8aMcXd3X/C//6k5y/LjRYv+vHnz2LFj8vtVAkBzwU1jQR0LC4tly5Z5Bger6fPxZ5/17dtXijtjATQ33DQWNCtVfe5zAc7iANAa3DQWNKv08cl/+23FdnGPHvmRkbqvB8BI4KaxwEju/PkZK1eKevZs2FgcHi799282ADQ7ph8Sym8a26tXL/lNY69evbpr1y6tVgZcUzJuXMm4cbzqar5YLLWz83zttZbfflsycqRM9TUAAOBF4Kax0DQyMzOpmRkR5c6Z0+4//xEePFj0xhtsFwVgmJgGdFBQ0C+//IKbxkK98oEDK/r1a/XttyVjxrBdC4BhYjoH3apVqwMHDtTV1Wm1GtAvufPmCUpKnH74ge1CAAwT04BOTk6eNm2ara2tm5ub+7+0WhlwX2XXrqUhIU67dztUVbFdC4ABYjrF8f3336u5HjQYrdx58+wuXHjl1i22CwEwQJqPoLt27UpEPXr06Nq16xdffOHu7t71X9ovD7iu2t29ODx8aGqq+YMHbNcCYGg0B3RycnL98rlz56rwzyw8K2/27GqBwHnzZrYLATA0TOegAVSpFQpjfXzsf/lF5SVJAeC5IKChGZzy9q5t2ZJwhx2AZsXoQ8KMjAz5Ql1dXXZ2diN/6LEAACAASURBVG1trfwhTuQAOYlAkBsZ6bZyJZ04QTgtGqCZaA5oqVTa8FrsPRtcjUEmk2mlKNBDxeHhbocO0aJFNHIkCZpwnXEAUEXzFEeNajqoD/SFjM+nlSspOZl272a7FgADoTmg33///draWoEyRFRVVYULQ8MTEybQgAH06ackFrNdCoAh0BzQbdq06dq16/r161NTUxu2379/f/369V27dsVMNDy1Zg1lZdGmTWzXAWAINAf0woULY2Njk5OT+/fvLxQKO3fu3Llz5xYtWvTt2zcpKSk2NnbhwoU6KBT0Q2AgjR1La9ZQYSHbpQDoPUan2Xl5eX366af5+fnXr1/fvXv37t27//jjj9jY2B07dnh5ealft7a2dtGiRXw+v6CgoL7xzJkzfn5+jo6OoaGhOTk5L/QKgGu++IJEIlq1iu06APQe00/bR44cefv2bU9PT09PTyKqqKjo27dvcXGxxhXDw8P9/f35/Kd/CUpLS6dMmXL8+PE+ffosX748KirqwIEDz1c9Z/H5/OTk5AULFmjsmZWVpYN6dMrbm958k775ht59lzp2ZLsaAD2mOaCjo6MXL16cl5fX8B5XNTU1w4cPZ7KD5cuX+/v7r1y5sr4lLi6uV69eAQEBRLRw4UJnZ2eJRGJubt704rnL0tLy3r17TG50LTDIM9JWrKB9+2jZMvrpJ7ZLAdBjmtPhrbfeevPNNydMmLB9+/b6RlNTU/ntvTXy9/dv1JKSklI/MWJvb+/g4JCWlqZxqkTv2NnZsV0Ce1xdae5cWrOG5s+nXr3YrgZAXzE6fOPxeOXl5U5OTs2yS7FYbGlpWf/QyspKrOysrEOHDn333XeNGnNzc4OCgpqlDNCujz6i77+nhQvpwgUqKyNrazIxYbsmAD3D9P9r+R1VwsLCGs4mPx9ra+uGE68VFRU2yu4MHRYWFhYW1qgxJiam4YeNwF22tjRnDi1dSjzek5aQEJo7l0aNYrUsAH3Cwh1VvL29b/17fffMzEyRSOTh4fF8mwLuys2lkyefaYmLo5dfpm++YakgAP3Dwh1VgoODIyMjL1y4MGjQoFWrVoWHhxvmB2VGbvFi+v13Je3vvUdjx1KDq7sAgCpMk7FHjx7PsfXCwkI3Nzcikkql8iPutLQ0Z2fnvXv3RkVF5eTkvPTSS9HR0c+xZeC06mpSM6yHDtG8eborBkBvMQ3odu3aqTmCbvQt8HqOjo5K78ASHBycmJjIcNdGgsfjnT59WuPbkpGR4ePjo5uSnsOjR49Wr14trKhYo7rPxe+/35+U1KNHj1mzZumuMgA9xDSg58yZ8+eff06dOrVly5a5ubnbtm0LCAgIDQ3VanFG5Z133mH4R4vLAZ2amlpVVRU8bhzt3auqTzsfn7Fjx27cuBEBDaAe04DevXt3QkJC/SkcISEh3bt3X4w7aDQfFxcXFxcXtqtoBnZ2ds5duoh79LC6cUNpB5OBA+UTXwCgHtOzOIqLi/Pz8+sf5uXlMfmeNxitoldfVdou9vevGDxYx8UA6CmmR9CRkZE+Pj6DBw+2t7cvKSm5fPkyLgMNapSMHm2ak+O8ceMzrTyeDF9XAWCMaUAvXbp0woQJ8fHxRUVFQqFwxYoV3bp102ploO/yp08vDQmxvXLF/NGjWnv7yq5dTUpK3D/5xPmrrzLmzGG7OgA9wDSgKysri4uLZ82aJRKJtmzZEhcX16FDB6XfAASoV922beGkSQ1bLJOSnHbtErdrx1ZJAHqE6Rz09OnTjx49SkTvvvtuXFzc7du3p0+frs3CwDDlLFok6tfP/f/+r3NpKdu1AHAd0yPoq1ev3r9/v6qq6siRI48ePXJwcOjUqZNWKwODJDMxebx+fYfXX//sn38oK4tcXdmuCIC7mB5Bm5iY8Hi88+fP9+rVSygUEpFEItFmYWCwpPb2DzdutJBKKTyc8FMEoBrTgA4ICBgxYkRkZOR7771HRMuWLVO80DMAQ1UdOqz19aXff6fISLZrAeCuJlws6fjx4y4uLgMGDCAiJyennTt3arMwMHC/tWpFS5fS559T79703ntslwPARUwD2szMLDw8vP7hHJwmBS/us88oKYnmzydfXxo6lO1qADjnRa++D/D8eDzatYu8vSk8nO7fZ7saAM5BQAOrbGzo+HHi82nCBBKJ2K4GgFtwpXymcnJyTDR9TbmwsLBFixa6qcdwtG9P+/bRyJE0ZQodOvT0FlkARg8BzcigQYOOHTumsRuPx2s4Uw9MBQXRmjX0wQe0fDm1aUPXr1NeHnl705gx9NJLbBcHwBoENCOTJk2a9OxXlqGZLVhA167RihXPNK5ZQ//5D23ejMNqME6YgwZuqKmhBvd6f2rLFmp0STwAo4GABm44eZKuXVP+VEwMyWS6rQaAEzDFASxLTU199OiR54EDHqp6XLt2JSampmXLgIAAS0tLHZYGwDIcQQPL5s2b9+OPP969dUtNn5OHDy9evPjy5cs6qwqACxDQwDKZTBYVFeUXFqamz9TFiwMDA2WY6AAjg4AGTigdPlzVU9Xt2vGkUl0WA8ARCGjgBKlQmL5+vWJ7jYuLWXp6hylTnEpKdF8VALvwIaEe69ix4/z58zV2MzMzs7Cw0EE9L6g0JKR6//6WO3faxcURkahvX1GPHgXvvGN186b7hx++v2dPcuvWNHIk22UC6A4CWo/t27eP7RKaWaWv72OF4+iKfv3u//yz+dSp/qtXU2Ehff01mZmxUh6AjmGKA/RAjbPz5oiIBxERtH07DRhADx+yXRGALiCgQT/U8fl3p02jI0coNZX69KHTp595Oi8PF8MDw4OABr0ybhwlJFCHDvTyyzR3LhUUUFQU8Xjk7Ew2NjRoEG3bhq8dgsHAHDTom3btKD6ePvyQNm2i3buptPTpU1eu0JUrdOcOffUVe/UBNBscQYMesrCgjRtpzJhn0rnexo30++86rwmg+SGgQW+dOKHyqZ9/1mEdANqCKQ7QM0VFRe+++66stna/6j7XY2LWZmZ6e3uvaHSBaQC9giNo0DP5+fllZWVvvfOOmj5tOnR4++23L126pKuiALQCAQ36x9ra2t3dvXzQIFUdnH//vf+WLZ5VVbqsCqDZIaBBX+WrOIiu8vYuDguz++WXndevU2AgnTjxzIl3xcW0ahWNHk1DhlB4OG3dSjU1OqoYoIkwBw36StyzZ/oXX7T58MNnGnv0yPz8c0m7drnz5p2fMuW9x49p7Fjy9KT33qPp0yk9nXx8ntnKoUMUE0OHD5O9vU6rB2AAAQ16rHTUqIqAALv4ePP796XW1lXe3hWBgTITEyKS2tgccnd/Ly6Ojh6lDRto3jxatoysrZVs5cIFmjePdu3SdfUAmiCgQb9JhcLiceNUPm1qShERFBFBf/1FS5ZQXJzybtHRtGEDOThoqUiA54OABsN3+/btpAcPPN3ceqjuc+Gbbwq9vF566SU3N7fGz126RLt2UXo6EVHHjjRzJvXpo71qAeohoMHwbd682dTU1CY7W00fp5iY321sRJmZb82b98wTn35Kn3/+9OHFi7RjB23cSHPmaKdYgKcQ0GAURowY4d2jB505o6qD39273SQSunqVvv6agoJowAAaMoTu3HkmnevNnUsvvUS9eyvfVnU1/fADXb5MmZnk6kqDB9PUqWRu3kwvBYwIAhqMhbh3b3HPnlZ//634VMmYMZkrVvy1a5d3Rkbf6mqKiaFt24hI+YeKclu3Kg/owkIKD6eG35HZs4d++okOHKBWrVRu7eFDio6mmzepupr8/GjSJPLz0/ySsrLo3j1yc6MOHYjP7JRZmYyqq/HXQl+wcx70mTNn/Pz8HB0dQ0NDc3JyWKkBjI2Mz8/85BPFdlGfPtmLF8sEghw3t4QRI+jECVle3qlPPvk9PFzNNabzz537YcOGQ4cONX5i5kxS/AZjfDzNmKGysoMHqUMHWrGCjh6l2Fj64gvq1k3DBfni42ngQHJzoyFDqFMnMjGhrVvV9Sei48dpxAji88nCggYPpk8/pcpKdf0lElq3joKDicejoUNp6lS6eVPDLlJT6d13adgw4vFo1CjavFnzOeYHDlB4OA0ZQsOG0YwZmneRl0cff0wjRhCPRy+/TBs2kMavI/3yC02aREOH0rBh9NZbdPWqhv7l5bRmDY0aRTwejRxJK1cqvyaXTrBwBF1aWjplypTjx4/36dNn+fLlUVFRBw4c0H0ZYIQknp5J8fGOe/ZY3b7Nq6mps7KqCAgoCg+XPXsPLXF19fz9+ydOnNj/4EFVm2qZljb1/fdFfD75+pKHB7VvT+3akbU1HT6sfIXjxyk1lTw9G7ffu0cREUr6z59P/v40ZIiSp65cocGDGzfOmkWlpbRokfK9r15NS5Y8fRgfT/Hx9OuvdOKE8v8SKipo/Hg6f/7JQ/mfnB9/pMOH6ZVXlO/il18oOPjpw9On6fRpOnGCDh9Wvou6Opo2jX744WmLfH7/xx9p8mTlu7h9+5l/LGJjKTaWjh6lw4fJ0VH5KkuW0OrVz7Ts3k3r1tGCBcr7Z2TQa689DfEzZ+jMGYqNpb17ycND+SraxEJAx8XF9erVKyAggIgWLlzo7OwskUjM8T8X6IRUKMx77z2N3VxcXMLDw8t+/dWuPqSeVREQUB4YGP/DD6+1b0+PHtVdvMivqFC/zfvvvpvr5+fZu3crLy9ycCChkOzt6dtvVa6waZPygFb2fwAR0Ycf0tSp1Lp14/bExGfSud7Fi7R6Na1cqeSpZctI6QufMIHy88nJqXF7Wdkz6VwvLo4++YS+/FLJU9u2PZPO9aZMoQEDqH37xu21tTR7tpL+8fE0Zw7t2aPkqdjYxuks98EHFBhI/fopeSoyUskh9rVrNHOmynM0tYmFgE5JSfHy8pIv29vbOzg4pKWl1bcAcEfRq6+qCujc+fMrfXzWX7r02rFjRLR61ap7f/wxobR0bHy8qq11jIvrqPBLLjMx4ala4ciRB2++KTM379C5M8/Wlng8cnCg/Hy6fFllxfv20VtvPX1oZkbW1rR3r8r+Sq8nJZUqj1S5gwdp1izFUlX237CBvviCTE0bt+9XfTnC6Gj67LPGjVeu0K+/Ku+/dy9t2qTkIHr7dpW72LZNSUCnplJsrPL+585RYiL5+qrcoHawENBisdjS0rL+oZWVlVgsVux28+bN3xUuux4fH19eXi4QPFP248ePjxw5Yqr4E/CsW7du1dbWpqen5+TkPHjw4KDq/17rZWVlbZN/WGQQrl69WlJSUqNpWrC4uPju3bvyF56VlcXkjUpNTT148OC1a9eSkpISExM1riKVSjMyMuS7aOrw5ebmNnX4kpKSBAKBk+Jxn8IurKysiEgikWRnZ8t3YTtkiK1CipWMGiVOTKTExPz8fPkurl+/bufgcLNt2xzVAV04aVJmRUVFbu5LPXpQZWXswYPtWrVqn5VlJZWqrOmHH4hI+Z8Ipd5/n95/n3l3+u23alNTIjI1N+dZWtbJZFVVVXyZzEL1GrVz5lQuXmxiZiZ/r6qqqmpra81FIjXjJ27Vqs7ExMrSkm9iIqurE4nFRGRTWKhyF2vWVH3zDZ/Pl++iurq6urratKpKzX/ZlR06SE1NLSws5OEgqqiQEVmWlZmo6C/98cfKY8d4PJ61tTUR1dTUSCQSQXW1mheePGtWvodHx44dXV1dG7bz+fwpU6ZoaQ6AhYC2trbOysqqf1hRUWFjY6PYzcLCQigUNmq0tbWtq6tr1G5qampra2v27DSiInNzcxsbG6FQKBaLTU1N7ezsNJZqamqqWIP+srKykkgkGl94TU2NmZmZ/IUzfKPMzMzs7OyEQqGtra25ubnGVaRSaf1729Thq6yslO9OY1X1u5Cvq3EVS0tLKysroVAokUgEAsGT/mPG1Pj6Wl+/LigqIh6vxslJ1L+/wN1dvi0TExP5LiwsLKysrKw6dLD08rJISVHcuKRTp4qePenBgwqBQDh0KBHdvnJl0H/+4xATY/3HH0rrqfL2zp8x46edO+fMnm1tafnHxYslJSW9nJycoqNVvQRRr15FrVo9ePBgcJ8+JJOdi4vz6tzZKSXFOi9P1Sq5Pj5ZWVmdOnWysbEpKioqSE9vJRQKlb0EObGVVZ2jY1FhoZunJxElJyfb2tpa1dbaqP6wTmZuXiYWm9na2traymQyEz6fZDJ71QFdIxCIra0FAoE8PTOzsgR8vq1AoPqUGuKbmpbV1ZkQCe3siMiESCaTWVdUqAq4Wj7f1Nqab2Jia2NDRHl5eXyZzFogUPMjYs7nFxQUVFZW+j57HM3j8fgMT6FpOhYC2tvb++zZs/LlzMxMkUjkoWz23cvLS3HeQyaTFRQURDz7oUp0dHRQUJCFhZo/fkREGRkZAwcOHDlyZHJy8rVr10JCQjSWeuLEiQilH+Dop4qKiqysLI0vPCMjIycnR/7Ct2zZwuSN+uuvv0JCQvz8/IRCYXFxscZVamtrz507J99FU4fv7t27v/32W5OG7/z584GBgUp/zBpV5eDgEBERIRKJtm/f/swuFP+pJyKi/fv3y3dx9+7d1q1b9+/f38zX1ys0VLFnyldfVbdt+/fff6empspXWbt2bUhIiJWzcwcVAZ0+bVppSEhsbOy4iROFQqGIz8/MzBw4dmxX1QGdsmLFQxOT2z//HLFzJxFtGTp0ztdft/j5Z1elE81EFf36PdqxY82aNa8vWtStW7fz58/fPnr09Zkzu3brpmoXGUuWFI4Y8cEHH3x67hwRjR49etmyZa0vX26zcKGqVRLj4r6PjvYNCho1atTT1mHD6OJF5St89BEtXVr/aPbs2aNGjepaXNz+zTdV7SLpxIkTV66QvX1EZOTT1vBwUjzNRu7NNxue97Jq1apWrVoNbNtW6djJ+X/wwbXi4nv37ukyE1g4zS44OPjOnTsXLlyora1dtWpVeHh4oykLAP1V7eaW9Ntv+TNnivr1I6KKfv3yp09PunKlum1bpf3F3bvnKvvQsnjChNKRI5WswOOlr12rdFNFr7+udC8lY8aoqrZwyhSlu8hRkbbi7t3Lhg9XbC8fOlTcvbvSVfL+8x+Z0gPMhnPljUybptgm8vcX9eqltHvx+PFSW1slTyj9UFFO2V/caje38qFDlXYvDwyUKH5uqX0sBLStre3evXujoqKcnZ3T09O/VPNxBIAektrZ5UZFPdyx4/atW4927MidO1eq9jJM+ZGRj7ZtKxs2TNSnDxGVDxqU/t//Zn72GfGUf3xYGhqaqXArr8JJk7KfvfJqvTorq9RDhxQDNPvDD8sVT9cjIqKCyZOLw8MV27OWLq1TNtlaZ26e/fHHSkvNV3UC+JQpyr8uf+wYKV4OhYj4/MzPPhMrHNpX9O+v6oXT8OH0f/+npP2776iH8uuyZC5bJu7Zs1GjuEePTMUPLXWCnUPX4ODgxMREVnYNwEEVAQEVAQHM+xe/8kppSIjVrVtmaWm1rVqJu3atbdlSTf8qL69HO3bYnT1rmZzMr6iQdOxYNnx4dZs2Klfg8zOXLSsNDbW9eNH80aNae/vKrl1LVB2oEhFRpY9PUnx8i4MHLRMTTcrLq93dywcPLhs6VNWfGeLxntya/eef6f59Mjen7t1p5kzq0EGxb0VFRVlZGQmFJWvXOh89anvrlkAkqhYKS/v1ywsNldXVUVlZZWWlveJFvT/6iIYOpehounePTEyoc2eaPp2UHeyLxeKysjIyNS1Zt875+HG7GzdMS0pqHBzKunfPHTu2zsKCysqUns6gVZhbANBLddbWFf37U//+TPtbWJSMG1ei5tKsCir69atQerKwClKhUOXxsipBQRQUpL5LQEDAkUan8bm4PFlITaXNm+ubhyqdoGDwLnXr1m3Tpk0nT5582mRqSvK/eXfu0J079c2zVHwaoSUIaONy/fp1JqfZ6aYYACamTp06depUre5izJgxY1TP1LMIAW1ERowYoTGdicjV1bW3quu0ATyrXbt2UVFRPFXzGP8qLy83pBOidAYBbURcXV1nzpzJdhVgUL755hu2SzBkuKs3AABHIaABADgKAQ0AwFGYgwaAp1JTU+M0XVdTqubqTtCsENAA8ETPnj1fUXU9/md99NFH2i4GCAENAPWEQiHO8+EUzEEDAHAUAhoAgKMwxQGgN+7fv694m6FGCgoK6pctLS1nMLg4Rm5urrXS+7oC2xDQAPrhpZdeSk5OTk1N1djztddeky/EqrrDHugJBDSAfvDy8lqzZg3bVYBOYQ4aAICjENAAAByFKQ5Qp2PHjvPnz9fYraampkWLFjqoB8CoIKBBnR07drBdgt54cs8ktUQikW6KAcOAgAZoBr169dq1a9evv/6qsefkyZPlC/n5+Z8xuBVpUlISX+ldscEIIKABmsHIkSNHjhzZpFWuXr1aWVmpsZtAIFByL1QwDghoAHa41N/5FEAF/OsEAMBRCGgAAI5CQAMAcJQxzkGbmZldvnz58ePHGnvizhEAwCJjDOgOHTrcvXu3rq5OY08rKysd1APa1rNnz40bNzLpuWTJEvlCenr6zp07NfZveOk4gGZnjAFNRDhvyajMnDmzSTcKsba23rVrV01NjcaeERERL1AXgAZGGtAA6g0ePJjtEgDwISEAAFfhCBqak4uLS3JyMpPrKwUEBOigHgC9hoCG5uTr63vu3Dm2qwAwEJjiAADgKAQ0AABHIaABADgKc9DAMoFA8NFHH2m85PHDhw9HjRolX87Nzf399981blkikTRDfQDsQUADy37++Wcml0UmIgcHByJyc3MbOHBgamqqxv7Tpk170eIAWIWABpZZWFhYWFgw729jY7Ny5Urt1QPAHZiDBgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwlNYDura2dtGiRXw+v+G9J86cOePn5+fo6BgaGpqTk6PtGgAA9JHWAzo8PNzKyqrh98RKS0unTJmybdu23Nzc3r17R0VFabsGAAB9pPUvqixfvtzf37/hNwvi4uJ69eolvxzwwoULnZ2dJRKJubm5tisBANAvWj+C9vf3b9SSkpLi5eUlX7a3t3dwcEhLS9N2GQAAeoeFr3qLxWJLS8v6h1ZWVmKxWLHboUOHvvvuu0aNRUVF48eP1259AADc0PwBvWnTphUrVhDRp59+OmfOHMUO1tbWWVlZ9Q8rKipsbGwUu4WFhYWFhTV7eQAA+qL5A/qdd9559dVXicjW1lZpB29v77Nnz8qXMzMzRSKRh4dHs5cBAKDvmn8O2traunXr1q1bt7a2tlbaITg4+M6dOxcuXKitrV21alV4eLhAgIvqAQA0pt0PCQsLC+UXk5RKpe7u7hYWFrm5uba2tnv37o2KinJ2dk5PT//yyy+1WgMAgJ7S7qGro6NjVVWVYntwcHBiYqJWdw0AoO/wVW8AAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHCUgO0CmkHLli3fe+89Ho+nvltBQUFYWJhuSgIAeHE8mUzGdg0AAKAEpjgAADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEcZYECfOXPGz8/P0dExNDQ0JydHY//a2tpFixbx+fyCggKGuzh+/LiPj4+Dg8OQIUNSUlI09t+/f7+Xl5eDg8OwYcNSU1MZ7oWILl26xOPxkpOTNfaUSCQ8Hs/iX6+++qrGVR4/fjx06FAbG5vu3bsnJCRo7B8dHW3RAI/HKywsVL9KTEyMr6+vp6dncHDwgwcPNO6CiA4fPuzj42NnZzdu3LjS0lJV3ZSOmpqhV9pf/dArfVb90CtdRc3oqylA1dArrqJ+6JXuQv3QK66ifuiV7kL90CtdRc3QK33b1f+mK11F/YgrXaWpv+zNTGZYSkpKnJycrl69WlNT8/HHH4eHh2tcZdy4ccuWLTMxMcnPz2eyi/T0dHt7+19//VUqlX788cdDhw5V3//u3buOjo43b96sra1duHBhcHAwo1cik1VVVfn7+zs7OyclJWnsnJ2d7eTkxHDLcgMHDly7dm11dfWuXbvefPPNJq177ty5YcOGqe+TkZHh4ODw4MEDmUy2fv364cOHa9zso0eP7O3t//7776qqqokTJ86ePVtVT8VRUz/0SkdZ/dArPqtx6BVXUT/6qgpQM/SKq6gfeqW7UD/06t8WxaFX7K9x6BVXUTP0St929cOtaqTUvDSlqzT1l73ZGVpAx8TEjBgxQr5cUlJibm5eVVWlfpUbN27IZLImBXRMTIx8+e+//3Zzc1Pf/+HDh6dPn5YvX7t2rW3btkz2IpPJli1b9umnn/r6+jIJ6OTk5I4dOzLcskwme/Dggbu7u1QqZb5KvZqaGj8/v5s3b6rvdvny5a5du8qXb9265ezsrHHLP/zww8svvyxfTkxMdHR0VNVTcdTUD73SUVY/9IrPahx6xVXUj76qAtQMveIq6odesb/GoVfztigdesX+GodecRU1Q6/0bVc/3KpGSs1LU7pKU3/Zm52hTXGkpKR4eXnJl+3t7R0cHNLS0tSv4u/v36RduLu7R0REyJfPnTsXGBiovr+Hh0doaCgRlZWVbd26dezYsUz2kpKScujQoY8++ohhVSUlJWKxeNiwYa1atQoODtb4v9g///zTqVOnyMhIDw+P4cOHM5lFqbdnzx5fX18/Pz/13bp3756fn//333/LZLKjZyX0bgAACz9JREFUR48GBwcz3wUROTg4FBYWlpSUKH1WcdTUD73SUVY/9IrPahx6xVXUj77SAtQPveIq6odesb/GoVfztigdesX+Gode/TvfaOiVvu3qh1vVSKnZr9JVmvrL3vx0/AdB25YsWbJo0aL6h+3bt5f/zdSI+RF0vTNnzrRv3z49PZ1J5w8++ICIBg4cWFhYyKT/sGHDzp07J5PJGB5BJyYmzpw5MykpqbKycvHixX5+fur779q1y9LS8vz583V1devXr9fYv6EuXbrIf/c0+umnnwQCgVAodHd3v3//vsb+aWlpdnZ2169fr6mpef/99/l8fnZ2tpr+DUeNydArHWX1Q6/0WfVDr7iK+tFv1J/J0DdchcnQN+zPcOiVvnA1Q9+oP5Ohb7gKk6Fv+LYz/E1XOlLqR1zpKk36ZW9GhhbQq1atevfdd+sftmzZ8t69e0xWbGpA79mzp3Pnzgw3LicSidatW9etW7e6ujr1PaOjoydNmiRfZhjQDVVXV5uZmWVmZqrpc/jwYX9/f/lybW2tmZlZQUEBk41fv37dx8eHSc+bN296eHg8evRIJpMdOHCgc+fOTGZUYmJiunbt6u3t/c033/D5fLFYrKZzw1FjMvTNEtAah17pBtWMfsP+DIdeVc2qhr5hf4ZDr7gL9UPfsD/DoW+0C/VD3+htZzLcqkZKzYgrXeU5ftmbi6EF9KFDhwYNGiRfzsjIsLKyqqmpYbJikwL66NGjfn5+6g/u6iUkJJw/f16+XFNTo/GoUCaTjR8/3snJydnZ2dnZWSAQODk5nThxQv0qWVlZiYmJ8uWqqiqNL+fGjRvt2rWTL9fW1pqampaUlDB5OZ9++unChQuZ9NywYcPEiRPrH5qZmTXpAOSPP/7w9vZW36fhy2Qy9C8e0EyGvuEqTEa/YX+GQ99wFSZD37CR4dArbkf90Dfsz3DoVb3zikOv+LZrHG41I6Vqv0pXadIve7MztIAuKytzcnI6f/58TU3N7Nmzp06dynBF5gFdVFTk7u7+8OFDhls+ffq0m5tbamqqTCaLjo52dnbWeATdEMMj6NjY2LZt2z58+LC2tnbp0qX9+vXTuEq3bt127txZV1f31Vdf9e7dm2E9Y8aM2bVrF5OecXFxHh4e8qOzuLg4R0dHjX8ss7OzPT09Hzx4UFFRERISsmbNGvX9G44ak6F/wYBmOPQNV2Ey+qoKYHgEzWToG+2CydArVqV+6Bv2Zzj0jc5FUTX0St929cOtfqSUvuFKV2nqL3uzM7SAlslkcXFxXbp0adGixejRozX+215QUGBubm5ubk5E8oWcnBz1q+zcuZPH45k3oHEv//3vf9u1a+fg4NCzZ8/Lly836eUwn+JYvXq1m5ubo6PjiBEjmPxIpaam9uzZ08HBYcCAAcnJyQzr6dGjR/1pCRqtWrXK09PT09Ozd+/eDF/4hg0bnJychELhrFmzVAW6qlFTNfRK+6sfeqXPqh96VRtUNfoaf/YUh17VKqqGXlV/NUOvpiqlQ6+qv5qhV7WKqqFX9bar+U1Xuor6N1zpKs/xy968eDKZTNefSwIAAAOGdpodAIDBQEADAHAUAhoAgKMQ0AAAHIWABgDgKAQ0AABHIaABADgKAQ0AwFEIaAAAjkJAAwBwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGp5TbW0tj8ebO3dufcvRo0fHjx//fFtLSEjw9PRsptKemj17tqOj49GjR+tb6urq1q9f7+3tbW5u3rp16xkzZuTl5TX7fgGaBQIanp+lpeWBAwfu3LnDdiEqxcTEXLt2reGfjXnz5m3dunXz5s35+fkXLlwoKioaNmxYbW0ti0UCqIKAhucnk8mWLVs2f/78Ru2///67t7d3o+WEhAR/f/9FixYNGTLE19f3woULYWFh3bt3j4qKkvfk8XgLFixo27Ztp06dLl68KG88efJkt27dOnbsGBQUlJ+fT0Q3btzo2bPnlClTgoODG+70wIEDXbt29fHxGTJkSHJyMhG98sorJSUlY8eOPXnypLxPenr6t99+e+TIkaCgIDs7uy5dusTExNTU1ERHRxPRTz/91KFDB3d398mTJ0skEsUWpa+rYT2NalMsXv4mfPzxxyEhId7e3mfPnpVvTXHXiutWV1dPmTKlY8eO7du3nzRpUmVlZXONI3CXLm8hDoakpqZGIBBIpdLu3bsfO3ZMJpMdOXJk3LhxMpns2rVrnTt3lnerX7516xafz//tt99kMtnEiRN9fX2rqqoqKyttbGyys7Nv3LghEAh27twpk8mio6M7deokk8kyMjJatGjxzz//yGSydevWjR8/Xr4da2vrffv2NSwmLS3Nzs7u7t27Mpls8+bN/fr1k7fb29unp6fXd9uzZ4+vr2+jF7J06dKwsLCHDx86Ojo+ePCgtrZ29OjRq1evVmxR9brq62m4rKp4Pp8fFxcnk8kOHDgQEBAgk8kUd6R03QMHDgQHB9fV1Uml0gULFvz666/NM5DAYQhoeE41NTUmJiYymezixYuenp4SiURjQDs6Osobly5dOnv2bPmyp6fn33//fePGDWtra6lUKpPJJBIJj8fLz8/funVrcHCwvFt5ebmpqWl1dfWtW7csLS3lPevt2LFjzJgx8uXKyko+n19aWipTCOgtW7YMHz680QvZsmXLkCFDtm3bNnbsWHmLSCSqqqpSbFH1uurrabisqnh7e3t5Y0JCQtu2bWUymeKOlK575coVNze3kydPVlZWNmGcQJ8J2D6CB703ZMiQbt26bdiwoXPnzup72tjYyBdMTEysrKzql6VSqUAgcHR05PP5RGRmZmZlZVVcXFxSUnL16lUPDw95T2tr64KCAiISCoXynvXy8/OdnJzkyxYWFlZWVnl5eXZ2do0KaNOmTVZWVqPG7OzsVq1aFRcXC4VCeYu8NsUWVRrWU7+sqviGb4JUKlW6I6XrBgYGbtiwYe3atRMnTgwLC/v666+tra3VVAUGAHPQ0AzWrVu3fv36+tMhTExM6urq5MtlZWUMN1JcXCyTyYiourpaJBI5Ojq6uroGBQU9+ldxcbGLiwsR8Xi8Rus6OzsXFhbKl6uqqsRisbOzs+Iu+vTpk5aWlpCQUN9SV1d38ODBl19+2cnJST7VK68kPT1dsUXV62pYT/2yquIVKe5I1boRERGXLl168ODB48ePt2zZovk9BT2HgIZm0L59++nTp//3v/+VP3RxccnKyhKJRER0+PBhhhsRiUSHDh0iop9//tnX17dFixbBwcHXrl1LSUkhouvXr9d/nKgoODg4Pj4+NTWViL777rsBAwbY2toqdnN2dv7ggw/Cw8MvXbpUVVWVmZk5ceJEKyuryZMnjxgx4rfffktMTJRKpTNnztyzZ49iS5NeF/PiFXekdN1NmzYtX75cJpO1aNGiTZs2in+lwPBgigOax5IlS3bv3i1fdnd3nzFjRkBAQJs2bYYPH15/SoYaNTU1vr6+v//++9KlS4no+++/JyJnZ+edO3eGh4dXVlba29t/9dVXqlZ3d3ffuXPnhAkTqqur27RpIz8rQ6nly5fb29tHRkY+evTIwcFh7NixZ8+e5fP5bm5u27dvDw0Nra2tDQoKmj9/vrm5uWIL89fFvHilu1Zc94033pg2bVq7du1MTEz69Okza9Ysje8q6Due/J9KAADgGkxxAABwFAIaAICjENAAAByFgAYA4CgENAAARyGgAQA4CgENAMBRCGgAAI5CQAMAcBQCGgCAoxDQAAAchYAGAOAoBDQAAEchoAEAOAoBDQDAUQhoAACOQkADAHAUAhoAgKP+HxqPwl++pSs1AAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Solution in Python"
      ],
      "metadata": {
        "id": "YsaYPWjIa-0W"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Python uses a different specification of negative binomial than we specified in R so we need to rewrite it as follows: \n",
        "\n",
        "$$\n",
        "p = \\frac{\\text{size}}{\\text{size} + \\mu}.\n",
        "$$"
      ],
      "metadata": {
        "id": "j_wQp-cqd3Ah"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "np.random.seed(1)\n",
        "N = 1000\n",
        "x = st.nbinom(n = 2, p = 2/(2+3)).rvs(N)\n",
        "np.mean(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tgric81PbBfu",
        "outputId": "282d2ce9-abf5-4c1b-c2ae-54fb6c8ca33a"
      },
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "3.09"
            ]
          },
          "metadata": {},
          "execution_count": 17
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "However, these pseudo-random numbers are not the same as in R. So, to make our example comparable with R, we read data generated from R that was saved to text file."
      ],
      "metadata": {
        "id": "RT1yAjlGs8OV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "x = np.loadtxt(\"nb_sim.txt\", dtype = np.int64)\n",
        "np.mean(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qRx74zKitNKD",
        "outputId": "2954b6e9-c149-4be0-e19c-603c9dbf2833"
      },
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "3.0006"
            ]
          },
          "metadata": {},
          "execution_count": 18
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Optimization using pmf function"
      ],
      "metadata": {
        "id": "feo9K8D8azao"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def pdf_nbinom(par, x):\n",
        "  pdfnbinom = st.nbinom(par[0],par[1]).logpmf(x)\n",
        "  return -np.sum(pdfnbinom)\n",
        "\n",
        "res = minimize(fun=pdf_nbinom, x0=[2, 0.5], args = (x), method = \"Nelder-Mead\")\n",
        "res"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FB_ouajwazwU",
        "outputId": "f01772e4-1d41-4657-d4d8-ffc83e4b8d86"
      },
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              " final_simplex: (array([[2.01357331, 0.40157875],\n",
              "       [2.0136646 , 0.40158661],\n",
              "       [2.01348695, 0.40156478]]), array([22061.74799012, 22061.74799024, 22061.74799191]))\n",
              "           fun: 22061.747990122178\n",
              "       message: 'Optimization terminated successfully.'\n",
              "          nfev: 76\n",
              "           nit: 40\n",
              "        status: 0\n",
              "       success: True\n",
              "             x: array([2.01357331, 0.40157875])"
            ]
          },
          "metadata": {},
          "execution_count": 19
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Goodness of fit -- step by step**"
      ],
      "metadata": {
        "id": "hbEQsnOkfTwY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "$G^2$ and $\\chi^2$ GoF can be calculated using `scipy.stats.power_divergence` function"
      ],
      "metadata": {
        "id": "wkzALQIDfV2p"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "x_uniq_vals, x_uniq_counts = np.unique(x, return_counts=True)\n",
        "## we simply use pdf(NB(2.01357331, 0.40157875), x) \n",
        "est_pdf = st.nbinom(res.x[0],res.x[1]).pmf(x_uniq_vals) \n",
        "## probs are not summing to 1 - rescale\n",
        "est_pdf = est_pdf/np.sum(est_pdf) "
      ],
      "metadata": {
        "id": "1T6GJ-MwX2P4"
      },
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "$G^2$ test"
      ],
      "metadata": {
        "id": "qFII4Sscuxim"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "st.power_divergence(x_uniq_counts, \n",
        "                    sum(x_uniq_counts)*est_pdf, \n",
        "                    lambda_ = 0, ddof = 2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Tr4KF7KubVHT",
        "outputId": "20f7a815-088c-4bc5-de83-967382bb5c3d"
      },
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Power_divergenceResult(statistic=14.378707061013243, pvalue=0.7040911371933183)"
            ]
          },
          "metadata": {},
          "execution_count": 21
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "$\\chi^2$ test"
      ],
      "metadata": {
        "id": "P-eMeMEDu1v5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "st.power_divergence(x_uniq_counts, \n",
        "                    sum(x_uniq_counts)*est_pdf, \n",
        "                    lambda_ = 1, ddof = 2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nhwbN077cCOh",
        "outputId": "7ab5a661-8707-4c5c-b3fd-39ad401627a9"
      },
      "execution_count": 22,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Power_divergenceResult(statistic=14.510966027945067, pvalue=0.695226354877517)"
            ]
          },
          "metadata": {},
          "execution_count": 22
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Rootograms** are not available in python. There are some scripts that recreate this plot using `matplotlib`, see\n",
        "\n",
        "https://stackoverflow.com/questions/38252879/how-to-plot-a-hanging-rootogram-in-python\n",
        "\n"
      ],
      "metadata": {
        "id": "Z-YxCNlDdfRD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Solution in Julia"
      ],
      "metadata": {
        "id": "IHcp32zRkCe-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "Random.seed!(123)\n",
        "n = 1000\n",
        "x = rand(NegativeBinomial(2, 2/(2+3)), n)\n",
        "mean(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AUiU6w4SkD9q",
        "outputId": "c551df69-c15f-42a7-e45f-b671f767081f"
      },
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "3.026"
            ]
          },
          "metadata": {},
          "execution_count": 23
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "However, these pseudo-random numbers are not the same as in R. So, to make our example comparable with R, we read data generated from R that was saved to text file."
      ],
      "metadata": {
        "id": "Zd6spYTothGW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "x = readdlm(\"nb_sim.txt\", Int)\n",
        "x = vec(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "RhBWVlggthw8",
        "outputId": "1665b76b-bd3e-4fb1-f812-d3cdcd946548"
      },
      "execution_count": 24,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([3, 1, 1, ..., 3, 3, 3], dtype=int64)"
            ]
          },
          "metadata": {},
          "execution_count": 24
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "function llnb(par, data)\n",
        "  ll = logpdf.(NegativeBinomial(par[1], par[2]), data)\n",
        "  return -sum(ll)\n",
        "end\n",
        "\n",
        "res = optimize(par -> llnb(par, x), [2, 0.5])\n",
        "res.minimizer"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SSLwLDL1kcaa",
        "outputId": "29df48cc-757b-4547-a66f-85cfc0fc2637"
      },
      "execution_count": 25,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([2.01360006, 0.40157958])"
            ]
          },
          "metadata": {},
          "execution_count": 25
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Goodness of fit statistics"
      ],
      "metadata": {
        "id": "egyxA_bwmeom"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "x_uniq_dict = sort(countmap(x))\n",
        "x_uniq_vals = Int.(keys(x_uniq_dict))\n",
        "x_uniq_counts = Int.(values(x_uniq_dict))\n",
        "params = res.minimizer\n",
        "est_pdf = pdf.(NegativeBinomial(params[1], params[2]), x_uniq_vals) ## not sum to 1 - rescale\n",
        "est_pdf = est_pdf ./ sum(est_pdf)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "b9TSRl6Vmh-L",
        "outputId": "5cd97f08-aee3-405d-d9a3-c0b52cc22a82"
      },
      "execution_count": 26,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([1.59332469e-01, 1.91992341e-01, 1.73119475e-01, 1.38600619e-01,\n",
              "       1.03958803e-01, 7.48224991e-02, 5.23393542e-02, 3.58562103e-02,\n",
              "       2.41757016e-02, 1.60965655e-02, 1.06088651e-02, 6.93355259e-03,\n",
              "       4.49964681e-03, 2.90262677e-03, 1.86274929e-03, 1.19003169e-03,\n",
              "       7.57253287e-04, 4.80174579e-04, 3.03527061e-04, 1.20296310e-04,\n",
              "       4.72387904e-05])"
            ]
          },
          "metadata": {},
          "execution_count": 26
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "$G^2$ test"
      ],
      "metadata": {
        "id": "1fHv42w590Hc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "PowerDivergenceTest(x_uniq_counts, lambda = 0.0, theta0 = est_pdf) ## throws an error on colab"
      ],
      "metadata": {
        "id": "WqTuP8Ss3m69"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "```julia\n",
        "Multinomial Likelihood Ratio Test\n",
        "---------------------------------\n",
        "Population details:\n",
        "    parameter of interest:   Multinomial Probabilities\n",
        "    value under h_0:         [0.159332, 0.191992, 0.173119, 0.138601, 0.103959, 0.0748225, 0.0523394, 0.0358562, 0.0241757, 0.0160966  …  0.00693355, 0.00449965, 0.00290263, 0.00186275, 0.00119003, 0.000757253, 0.000480175, 0.000303527, 0.000120296, 4.72388e-5]\n",
        "    point estimate:          [0.1583, 0.1969, 0.1699, 0.1356, 0.1025, 0.0748, 0.0565, 0.0357, 0.0244, 0.0164  …  0.0054, 0.0045, 0.0024, 0.0023, 0.0012, 0.0011, 0.0006, 0.0002, 0.0001, 0.0001]\n",
        "    95% confidence interval: [(0.1491, 0.1676), (0.1877, 0.2062), (0.1607, 0.1792), (0.1264, 0.1449), (0.0933, 0.1118), (0.0656, 0.0841), (0.0473, 0.0658), (0.0265, 0.045), (0.0152, 0.0337), (0.0072, 0.0257)  …  (0.0, 0.0147), (0.0, 0.0138), (0.0, 0.0117), (0.0, 0.0116), (0.0, 0.0105), (0.0, 0.0104), (0.0, 0.0099), (0.0, 0.0095), (0.0, 0.0094), (0.0, 0.0094)]\n",
        "\n",
        "Test summary:\n",
        "    outcome with 95% confidence: fail to reject h_0\n",
        "    one-sided p-value:           0.8108\n",
        "\n",
        "Details:\n",
        "    Sample size:        10000\n",
        "    statistic:          14.378553064834735\n",
        "    degrees of freedom: 20\n",
        "    residuals:          [-0.258657, 1.12004, -0.77377, -0.805987, -0.452445, -0.00822523, 1.81864, -0.082495, 0.144257, 0.239165  …  -1.84171, 0.000526521, -0.932932, 1.0131, 0.0288963, 1.24553, 0.546826, -0.594231, -0.185051, 0.767654]\n",
        "    std. residuals:     [-0.282106, 1.24602, -0.850924, -0.868412, -0.477971, -0.00855137, 1.86819, -0.0840149, 0.146033, 0.241114  …  -1.84813, 0.000527709, -0.934289, 1.01405, 0.0289135, 1.246, 0.546958, -0.594321, -0.185062, 0.767672]\n",
        "```"
      ],
      "metadata": {
        "id": "fmGbs30p-Eg3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Correct p-value for the correct number of dof."
      ],
      "metadata": {
        "id": "98K2lcXZ-JCl"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "1-cdf(Chisq(18), 14.378553064834735)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7rIFaWj3-OZJ",
        "outputId": "1bd7b185-0194-445d-c2fa-5a1e495eb981"
      },
      "execution_count": 29,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.704101421610712"
            ]
          },
          "metadata": {},
          "execution_count": 29
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "$\\chi^2$ test"
      ],
      "metadata": {
        "id": "Azip379k93NP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "PowerDivergenceTest(x_uniq_counts, lambda = 1.0, theta0 = est_pdf) ## throws an error on colab"
      ],
      "metadata": {
        "id": "gD8P8r509xi1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "```julia\n",
        "Pearson's Chi-square Test\n",
        "-------------------------\n",
        "Population details:\n",
        "    parameter of interest:   Multinomial Probabilities\n",
        "    value under h_0:         [0.159332, 0.191992, 0.173119, 0.138601, 0.103959, 0.0748225, 0.0523394, 0.0358562, 0.0241757, 0.0160966  …  0.00693355, 0.00449965, 0.00290263, 0.00186275, 0.00119003, 0.000757253, 0.000480175, 0.000303527, 0.000120296, 4.72388e-5]\n",
        "    point estimate:          [0.1583, 0.1969, 0.1699, 0.1356, 0.1025, 0.0748, 0.0565, 0.0357, 0.0244, 0.0164  …  0.0054, 0.0045, 0.0024, 0.0023, 0.0012, 0.0011, 0.0006, 0.0002, 0.0001, 0.0001]\n",
        "    95% confidence interval: [(0.1491, 0.1676), (0.1877, 0.2062), (0.1607, 0.1792), (0.1264, 0.1449), (0.0933, 0.1118), (0.0656, 0.0841), (0.0473, 0.0658), (0.0265, 0.045), (0.0152, 0.0337), (0.0072, 0.0257)  …  (0.0, 0.0147), (0.0, 0.0138), (0.0, 0.0117), (0.0, 0.0116), (0.0, 0.0105), (0.0, 0.0104), (0.0, 0.0099), (0.0, 0.0095), (0.0, 0.0094), (0.0, 0.0094)]\n",
        "\n",
        "Test summary:\n",
        "    outcome with 95% confidence: fail to reject h_0\n",
        "    one-sided p-value:           0.8037\n",
        "\n",
        "Details:\n",
        "    Sample size:        10000\n",
        "    statistic:          14.510593683125506\n",
        "    degrees of freedom: 20\n",
        "    residuals:          [-0.258657, 1.12004, -0.77377, -0.805987, -0.452445, -0.00822523, 1.81864, -0.082495, 0.144257, 0.239165  …  -1.84171, 0.000526521, -0.932932, 1.0131, 0.0288963, 1.24553, 0.546826, -0.594231, -0.185051, 0.767654]\n",
        "    std. residuals:     [-0.282106, 1.24602, -0.850924, -0.868412, -0.477971, -0.00855137, 1.86819, -0.0840149, 0.146033, 0.241114  …  -1.84813, 0.000527709, -0.934289, 1.01405, 0.0289135, 1.246, 0.546958, -0.594321, -0.185062, 0.767672]\n",
        "```"
      ],
      "metadata": {
        "id": "x5U_rhP2zFak"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Note that p-value is calculated wrongly as correct df is 18 not 20. So correct p-value is: "
      ],
      "metadata": {
        "id": "AWn0Qdrfqmsi"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%%julia\n",
        "1-cdf(Chisq(18), 14.510593683125506)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PMnRSVQ8rune",
        "outputId": "8ad0ce83-489e-46b9-ec25-2e08b67059a5"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.6952513989555996"
            ]
          },
          "metadata": {},
          "execution_count": 72
        }
      ]
    }
  ]
}