{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Mixed Effects Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: The R code and the results in this notebook has been converted to markdown so that R is not required to build the documents. The R results in the notebook were computed using R 3.5.1 and lme4 1.1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%load_ext rpy2.ipython\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%R library(lme4)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "array(['lme4', 'Matrix', 'tools', 'stats', 'graphics', 'grDevices',\n",
    "       'utils', 'datasets', 'methods', 'base'], dtype='<U9')\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing R lmer to statsmodels MixedLM\n",
    "=======================================\n",
    "\n",
    "The statsmodels imputation of linear mixed models (MixedLM) closely follows the approach outlined in Lindstrom and Bates (JASA 1988).  This is also the approach followed in the  R package LME4.  Other packages such as Stata, SAS, etc. should also be consistent with this approach, as the basic techniques in this area are mostly mature.\n",
    "\n",
    "Here we show how linear mixed models can be fit using the MixedLM procedure in statsmodels.  Results from R (LME4) are included for comparison.\n",
    "\n",
    "Here are our import statements:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Growth curves of pigs\n",
    "\n",
    "These are longitudinal data from a factorial experiment. The outcome variable is the weight of each pig, and the only predictor variable we will use here is \"time\".  First we fit a model that expresses the mean weight as a linear function of time, with a random intercept for each pig. The model is specified using formulas. Since the random effects structure is not specified, the default random effects structure (a random intercept for each group) is automatically used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         Mixed Linear Model Regression Results\n",
      "========================================================\n",
      "Model:            MixedLM Dependent Variable: Weight    \n",
      "No. Observations: 861     Method:             REML      \n",
      "No. Groups:       72      Scale:              11.3669   \n",
      "Min. group size:  11      Log-Likelihood:     -2404.7753\n",
      "Max. group size:  12      Converged:          Yes       \n",
      "Mean group size:  12.0                                  \n",
      "--------------------------------------------------------\n",
      "             Coef.  Std.Err.    z    P>|z| [0.025 0.975]\n",
      "--------------------------------------------------------\n",
      "Intercept    15.724    0.788  19.952 0.000 14.179 17.268\n",
      "Time          6.943    0.033 207.939 0.000  6.877  7.008\n",
      "Group Var    40.394    2.149                            \n",
      "========================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "data = sm.datasets.get_rdataset('dietox', 'geepack', cache=True).data\n",
    "md = smf.mixedlm(\"Weight ~ Time\", data, groups=data[\"Pig\"])\n",
    "mdf = md.fit()\n",
    "print(mdf.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the same model fit in R using LMER:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%%R\n",
    "data(dietox, package='geepack')\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%R print(summary(lmer('Weight ~ Time + (1|Pig)', data=dietox)))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Linear mixed model fit by REML ['lmerMod']\n",
    "Formula: Weight ~ Time + (1 | Pig)\n",
    "   Data: dietox\n",
    "\n",
    "REML criterion at convergence: 4809.6\n",
    "\n",
    "Scaled residuals: \n",
    "    Min      1Q  Median      3Q     Max \n",
    "-4.7118 -0.5696 -0.0943  0.4877  4.7732 \n",
    "\n",
    "Random effects:\n",
    " Groups   Name        Variance Std.Dev.\n",
    " Pig      (Intercept) 40.39    6.356   \n",
    " Residual             11.37    3.371   \n",
    "Number of obs: 861, groups:  Pig, 72\n",
    "\n",
    "Fixed effects:\n",
    "            Estimate Std. Error t value\n",
    "(Intercept) 15.72352    0.78805   19.95\n",
    "Time         6.94251    0.03339  207.94\n",
    "\n",
    "Correlation of Fixed Effects:\n",
    "     (Intr)\n",
    "Time -0.275\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the statsmodels summary of results, the fixed effects and random effects parameter estimates are shown in a single table.  The random effect for animal is labeled \"Intercept RE\" in the statsmodels output above.  In the LME4 output, this effect is the pig intercept under the random effects section.\n",
    "\n",
    "There has been a lot of debate about whether the standard errors for random effect variance and covariance parameters are useful.  In LME4, these standard errors are not displayed, because the authors of the package believe they are not very informative.  While there is good reason to question their utility, we elected to include the standard errors in the summary table, but do not show the corresponding Wald confidence intervals.\n",
    "\n",
    "Next we fit a model with two random effects for each animal: a random intercept, and a random slope (with respect to time).  This means that each pig may have a different baseline weight, as well as growing at a different rate. The formula specifies that \"Time\" is a covariate with a random coefficient.  By default, formulas always include an intercept (which could be suppressed here using \"0 + Time\" as the formula)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:567: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warn(\"Maximum Likelihood optimization failed to converge. \"\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2112: ConvergenceWarning: Retrying MixedLM optimization with lbfgs\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:567: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warn(\"Maximum Likelihood optimization failed to converge. \"\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2112: ConvergenceWarning: Retrying MixedLM optimization with cg\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           Mixed Linear Model Regression Results\n",
      "===========================================================\n",
      "Model:             MixedLM  Dependent Variable:  Weight    \n",
      "No. Observations:  861      Method:              REML      \n",
      "No. Groups:        72       Scale:               5.7891    \n",
      "Min. group size:   11       Log-Likelihood:      -2220.3890\n",
      "Max. group size:   12       Converged:           No        \n",
      "Mean group size:   12.0                                    \n",
      "-----------------------------------------------------------\n",
      "                 Coef.  Std.Err.   z    P>|z| [0.025 0.975]\n",
      "-----------------------------------------------------------\n",
      "Intercept        15.739    0.672 23.438 0.000 14.423 17.055\n",
      "Time              6.939    0.085 81.326 0.000  6.772  7.106\n",
      "Group Var        30.266    4.271                           \n",
      "Group x Time Cov  0.746    0.304                           \n",
      "Time Var          0.483    0.046                           \n",
      "===========================================================\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:567: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warn(\"Maximum Likelihood optimization failed to converge. \"\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2118: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2130: ConvergenceWarning: Gradient optimization failed, |grad| = 31.645802\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "md = smf.mixedlm(\"Weight ~ Time\", data, groups=data[\"Pig\"], re_formula=\"~Time\")\n",
    "mdf = md.fit()\n",
    "print(mdf.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the same model fit using LMER in R:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%R print(summary(lmer(\"Weight ~ Time + (1 + Time | Pig)\", data=dietox)))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Linear mixed model fit by REML ['lmerMod']\n",
    "Formula: Weight ~ Time + (1 + Time | Pig)\n",
    "   Data: dietox\n",
    "\n",
    "REML criterion at convergence: 4434.1\n",
    "\n",
    "Scaled residuals: \n",
    "    Min      1Q  Median      3Q     Max \n",
    "-6.4286 -0.5529 -0.0416  0.4841  3.5624 \n",
    "\n",
    "Random effects:\n",
    " Groups   Name        Variance Std.Dev. Corr\n",
    " Pig      (Intercept) 19.493   4.415        \n",
    "          Time         0.416   0.645    0.10\n",
    " Residual              6.038   2.457        \n",
    "Number of obs: 861, groups:  Pig, 72\n",
    "\n",
    "Fixed effects:\n",
    "            Estimate Std. Error t value\n",
    "(Intercept) 15.73865    0.55012   28.61\n",
    "Time         6.93901    0.07982   86.93\n",
    "\n",
    "Correlation of Fixed Effects:\n",
    "     (Intr)\n",
    "Time 0.006 \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The random intercept and random slope are only weakly correlated $(0.294 / \\sqrt{19.493 * 0.416} \\approx 0.1)$.  So next we fit a model in which the two random effects are constrained to be uncorrelated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10324316832591753"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    ".294 / (19.493 * .416)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:567: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warn(\"Maximum Likelihood optimization failed to converge. \"\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2112: ConvergenceWarning: Retrying MixedLM optimization with lbfgs\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:567: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warn(\"Maximum Likelihood optimization failed to converge. \"\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2112: ConvergenceWarning: Retrying MixedLM optimization with cg\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           Mixed Linear Model Regression Results\n",
      "===========================================================\n",
      "Model:             MixedLM  Dependent Variable:  Weight    \n",
      "No. Observations:  861      Method:              REML      \n",
      "No. Groups:        72       Scale:               5.8015    \n",
      "Min. group size:   11       Log-Likelihood:      -2220.0996\n",
      "Max. group size:   12       Converged:           No        \n",
      "Mean group size:   12.0                                    \n",
      "-----------------------------------------------------------\n",
      "                 Coef.  Std.Err.   z    P>|z| [0.025 0.975]\n",
      "-----------------------------------------------------------\n",
      "Intercept        15.739    0.672 23.416 0.000 14.421 17.056\n",
      "Time              6.939    0.084 83.012 0.000  6.775  7.103\n",
      "Group Var        30.322    4.025                           \n",
      "Group x Time Cov  0.000    0.000                           \n",
      "Time Var          0.462    0.040                           \n",
      "===========================================================\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/base/model.py:567: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warn(\"Maximum Likelihood optimization failed to converge. \"\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2118: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2130: ConvergenceWarning: Gradient optimization failed, |grad| = 20.942980\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "md = smf.mixedlm(\"Weight ~ Time\", data, groups=data[\"Pig\"],\n",
    "                  re_formula=\"~Time\")\n",
    "free = sm.regression.mixed_linear_model.MixedLMParams.from_components(np.ones(2),\n",
    "                                                                      np.eye(2))\n",
    "\n",
    "mdf = md.fit(free=free)\n",
    "print(mdf.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood drops by 0.3 when we fix the correlation parameter to 0.  Comparing 2 x 0.3 = 0.6 to the chi^2 1 df reference distribution suggests that the data are very consistent with a model in which this parameter is equal to 0.\n",
    "\n",
    "Here is the same model fit using LMER in R (note that here R is reporting the REML criterion instead of the likelihood, where the REML criterion is twice the log likelihood):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%R print(summary(lmer(\"Weight ~ Time + (1 | Pig) + (0 + Time | Pig)\", data=dietox)))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Linear mixed model fit by REML ['lmerMod']\n",
    "Formula: Weight ~ Time + (1 | Pig) + (0 + Time | Pig)\n",
    "   Data: dietox\n",
    "\n",
    "REML criterion at convergence: 4434.7\n",
    "\n",
    "Scaled residuals: \n",
    "    Min      1Q  Median      3Q     Max \n",
    "-6.4281 -0.5527 -0.0405  0.4840  3.5661 \n",
    "\n",
    "Random effects:\n",
    " Groups   Name        Variance Std.Dev.\n",
    " Pig      (Intercept) 19.8404  4.4543  \n",
    " Pig.1    Time         0.4234  0.6507  \n",
    " Residual              6.0282  2.4552  \n",
    "Number of obs: 861, groups:  Pig, 72\n",
    "\n",
    "Fixed effects:\n",
    "            Estimate Std. Error t value\n",
    "(Intercept) 15.73875    0.55444   28.39\n",
    "Time         6.93899    0.08045   86.25\n",
    "\n",
    "Correlation of Fixed Effects:\n",
    "     (Intr)\n",
    "Time -0.086\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sitka growth data\n",
    "\n",
    "This is one of the example data sets provided in the LMER R library.  The outcome variable is the size of the tree, and the covariate used here is a time value.  The data are grouped by tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sm.datasets.get_rdataset(\"Sitka\", \"MASS\", cache=True).data\n",
    "endog = data[\"size\"]\n",
    "data[\"Intercept\"] = 1\n",
    "exog = data[[\"Intercept\", \"Time\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the statsmodels LME fit for a basic model with a random intercept.  We are passing the endog and exog data directly to the LME init function as arrays.  Also note that endog_re is specified explicitly in argument 4 as a random intercept (although this would also be the default if it were not specified)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         Mixed Linear Model Regression Results\n",
      "=======================================================\n",
      "Model:             MixedLM Dependent Variable: size    \n",
      "No. Observations:  395     Method:             REML    \n",
      "No. Groups:        79      Scale:              0.0392  \n",
      "Min. group size:   5       Log-Likelihood:     -82.3884\n",
      "Max. group size:   5       Converged:          Yes     \n",
      "Mean group size:   5.0                                 \n",
      "-------------------------------------------------------\n",
      "              Coef. Std.Err.   z    P>|z| [0.025 0.975]\n",
      "-------------------------------------------------------\n",
      "Intercept     2.273    0.088 25.864 0.000  2.101  2.446\n",
      "Time          0.013    0.000 47.796 0.000  0.012  0.013\n",
      "Intercept Var 0.374    0.345                           \n",
      "=======================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "md = sm.MixedLM(endog, exog, groups=data[\"tree\"], exog_re=exog[\"Intercept\"])\n",
    "mdf = md.fit()\n",
    "print(mdf.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the same model fit in R using LMER:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%R\n",
    "data(Sitka, package=\"MASS\")\n",
    "print(summary(lmer(\"size ~ Time + (1 | tree)\", data=Sitka)))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "```\n",
    "Linear mixed model fit by REML ['lmerMod']\n",
    "Formula: size ~ Time + (1 | tree)\n",
    "   Data: Sitka\n",
    "\n",
    "REML criterion at convergence: 164.8\n",
    "\n",
    "Scaled residuals: \n",
    "    Min      1Q  Median      3Q     Max \n",
    "-2.9979 -0.5169  0.1576  0.5392  4.4012 \n",
    "\n",
    "Random effects:\n",
    " Groups   Name        Variance Std.Dev.\n",
    " tree     (Intercept) 0.37451  0.612   \n",
    " Residual             0.03921  0.198   \n",
    "Number of obs: 395, groups:  tree, 79\n",
    "\n",
    "Fixed effects:\n",
    "             Estimate Std. Error t value\n",
    "(Intercept) 2.2732443  0.0878955   25.86\n",
    "Time        0.0126855  0.0002654   47.80\n",
    "\n",
    "Correlation of Fixed Effects:\n",
    "     (Intr)\n",
    "Time -0.611\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now try to add a random slope.  We start with R this time.  From the code and output below we see that the REML estimate of the variance of the random slope is nearly zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ipython\n",
    "%R print(summary(lmer(\"size ~ Time + (1 + Time | tree)\", data=Sitka)))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Linear mixed model fit by REML ['lmerMod']\n",
    "Formula: size ~ Time + (1 + Time | tree)\n",
    "   Data: Sitka\n",
    "\n",
    "REML criterion at convergence: 153.4\n",
    "\n",
    "Scaled residuals: \n",
    "    Min      1Q  Median      3Q     Max \n",
    "-2.7609 -0.5173  0.1188  0.5270  3.5466 \n",
    "\n",
    "Random effects:\n",
    " Groups   Name        Variance  Std.Dev. Corr \n",
    " tree     (Intercept) 2.217e-01 0.470842      \n",
    "          Time        3.288e-06 0.001813 -0.17\n",
    " Residual             3.634e-02 0.190642      \n",
    "Number of obs: 395, groups:  tree, 79\n",
    "\n",
    "Fixed effects:\n",
    "            Estimate Std. Error t value\n",
    "(Intercept) 2.273244   0.074655   30.45\n",
    "Time        0.012686   0.000327   38.80\n",
    "\n",
    "Correlation of Fixed Effects:\n",
    "     (Intr)\n",
    "Time -0.615\n",
    "convergence code: 0\n",
    "Model failed to converge with max|grad| = 0.793203 (tol = 0.002, component 1)\n",
    "Model is nearly unidentifiable: very large eigenvalue\n",
    " - Rescale variables?\n",
    " ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run this in statsmodels LME with defaults, we see that the variance estimate is indeed very small, which leads to a warning about the solution being on the boundary of the parameter space.  The regression slopes agree very well with R, but the likelihood value is much higher than that returned by R."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             Mixed Linear Model Regression Results\n",
      "===============================================================\n",
      "Model:               MixedLM    Dependent Variable:    size    \n",
      "No. Observations:    395        Method:                REML    \n",
      "No. Groups:          79         Scale:                 0.0264  \n",
      "Min. group size:     5          Log-Likelihood:        -62.4834\n",
      "Max. group size:     5          Converged:             Yes     \n",
      "Mean group size:     5.0                                       \n",
      "---------------------------------------------------------------\n",
      "                     Coef.  Std.Err.   z    P>|z| [0.025 0.975]\n",
      "---------------------------------------------------------------\n",
      "Intercept             2.273    0.101 22.513 0.000  2.075  2.471\n",
      "Time                  0.013    0.000 33.888 0.000  0.012  0.013\n",
      "Intercept Var         0.646    0.914                           \n",
      "Intercept x Time Cov -0.001    0.003                           \n",
      "Time Var              0.000    0.000                           \n",
      "===============================================================\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "exog_re = exog.copy()\n",
    "md = sm.MixedLM(endog, exog, data[\"tree\"], exog_re)\n",
    "mdf = md.fit()\n",
    "print(mdf.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can further explore the random effects structure by constructing plots of the profile likelihoods. We start with the random intercept, generating a plot of the profile likelihood from 0.1 units below to 0.1 units above the MLE. Since each optimization inside the profile likelihood generates a warning (due to the random slope variance being close to zero), we turn off the warnings here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    likev = mdf.profile_re(0, 're', dist_low=0.1, dist_high=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a plot of the profile likelihood function.  We multiply the log-likelihood difference by 2 to obtain the usual $\\chi^2$ reference distribution with 1 degree of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '-2 times profile log likelihood')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn4AAAHoCAYAAADXB+KjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3iV5f3H8fc3mxEIkLB3CHsTQKYDERURcI866q6z1mptq1atttra+nMrWm3Vat2g4kRZgoBhb0LCngmBEAjZ9++PHGyKIeTASZ7knM/rup6LnOc55+ETLi/zzX0/3/s25xwiIiIiEvzCvA4gIiIiItVDhZ+IiIhIiFDhJyIiIhIiVPiJiIiIhAgVfiIiIiIhQoWfiIiISIiI8DpAbRAfH+/at2/vdQwRERGRY1q4cGGmcy6hvGsq/Cqhffv2pKSkeB1DRERE5JjMbNPRrmmqV0RERCREqPATERERCREq/ERERERChAo/ERERkRChwk9EREQkRKjwExEREQkRKvxEREREQoQKPxEREZEQocJPREREJESo8BMREREJESr8REREREKECj8RERGREKHCT0RERCREqPATERERCREq/ERERERChAo/ERERkRChwk9ExEPOOQqLS7yOISIhIsLrACIioSY7t5A5aZnMWpfB7NRMdu7Po0/rhgzrFM/QxHj6t4sjOiLc65giEoRU+ImIVLGi4hKWbNnHrNRMZqdmsHTLPkocxEZHMLRTE8b2bsGCDVk8N309z3y7nuiIMAa2b8zQTk0YlhhPz1YNCQ8zr78NEQkCKvxERKrAlqxcZqVmMGtdBnPX7yEnv4gwg96t47j1tCRGJsXTt00cEeH/feJmf14h89OzmJuWydz1e/jLF2uBtcTGRHBSxyYMS2zC0E7xJDWtj5kKQRHxnwo/EZEAOJBfxLy0PcxKLZ2+3ZB5EICWDWMY27sFI5ISGNapCXF1o456jwYxkYzu3ozR3ZsBkJGTz9y0TL5P28OctEy+XrULgITYaIYmlo4GDklsQpvGdav+GxSRoGDOOa8z1HjJyckuJSXF6xgiUoOUlDhWbt//46jeos17KSx21IkM56SOjRnZOYERSQkkJtQL2Ojclqxc5qZlMmf9Huam7SHzQD4AbRvXZVinJgz1FYLx9aMD8veJSO1kZgudc8nlXlPhd2wq/EQEYNf+vB8bMr5bn0nWwQIAurdowMjOCYxMimdA+0bV0pjhnCN19wHmrC8tBOenl04nA3RtHsvQxHiGJjZhcMfGxMZEVnkeEak5VPidIBV+IqEpr7CYBRuyfiz21u7KASC+fjQjk+IZ0Tme4Z0SSIj1foStqLiEFdv3M2d9JnPTMknZuJf8ohLCw4zerRv+ODXcv10jYiLVMSwSzFT4nSAVfiKhwTnHul0HmLUug1mpGSzYkEV+UQlR4WEM7NCIEUkJjExKoGvzWMJqeJdtXmExizbvZe760ucDl23NprjEER0RRnL7Rj+OCPZq1fB/GkxEpPZT4XeCVPiJBK+sgwXM9jVkzE7NYNf+0ufmOjWtz8ikBEZ0juekDk2oE1W7R8ly8gpZsCHL93xgJmt2lo5exkZHMLhjk9IRwU7xdG6mjmGR2q6iwk9dvSISUgqKSli0eS+zUzOYtS6TFduzcQ4a1olkeFJ86RRuUgIt4+p4HTWgYmMiGdWtGaO6lXYMZx7I5/u00iaRuWmZTFtd2jEcXz+KIYnxDPMVguoYFgkuGvGrBI34idRezjk27sn1PaeXwfdpezhYUEx4mNG/bVzp9G3nBHqF+CLJW/fmMtc3GjgnbQ8ZOaUjn20a12Fox3iG+rqGa8LzjCJSMU31niAVfiK1y/68Quauz2RWaum2aFv3HgJKi5iRvkJvSGITGqjbtVzOOdbvPsDctD3MWZ/JvPQ97M8r7Rju3Kw+QxPjGdYpnsEdG+vfUKQGUuF3glT4idRsxSWOpVv3MXtdJrNSM1iyZR/FJY760REMSWzy4/Rt+/h6XketlYpLHCu2Zf84LfzDxizyCksIM+jVOu7HaeEB6hgWqRFU+J0gFX4iNc/2fYd+7L79LjWT/XlFmEHvVg1/nL7t1zaOSHWsBlx+UTGLN+9j7vrSaeGlW/ZRVOKIigjjouTW3De2uwpAEQ+puUNEgoZzjkmz0nnsizU4B80aRDOmR3NGdk5gWKd4Gtc7+pZoEhjREeGc1LEJJ3Vswq8o3a7uhw1ZfLVqJ2/O28zCTft47rJ+dEyo73VUETmCRvwqQSN+IjVDYXEJD0xZwdsLtjC2VwvuOD2JpKZafqQmmb5mN3e+u4SiYsdj5/finN4tvY4kEnIqGvHTHIiI1Ar78wq55p8/8PaCLdxyaiLPXNqPzs1iVfTVMKd2bcrU20fQuVl9bn1rMfdPXkF+UbHXsUTER4WfiNR4W7JyOf/5uXyftoe/XNCbu8d0rfE7Z4SyVnF1eOfGIVw/ogNvzNvE+S/MZfOeXK9jiQgq/ESkhlu8eS8Tn5/Drv15vH7NIC5KbuN1JKmEyPAwfj+2O5OuGMDmPbmMfWY2X6zY6XUskZCnwk9EaqzPlu/gkknzqBMVzoc3D2Nop3ivI4mfzujRnKm3j6BjfD1uenMhD32ykoKiEq9jiYQsFX4iUuM453hxZho3/3sRPVo2YPLNw+jUVB2itVWbxnV596YhXD20Pa/N2ciFL33P1r2a+hXxggo/EalRCotL+O2Hy3ns8zWc07sFb11/Ek3qa5uw2i46IpwHz+3B85f3J333AcY+/R3TVu3yOpZIyKmRhZ+ZXWhmK82sxMySy5wfbWYLzWy578/TylybYWZrzWyJ72h6lHv/1szW+947pjq+HxGpnOxDhfz8tR/4zw9buPXUTjx9ST8tBBxkzu7Vgk9uG06ruDpc93oKf/psNYXFmvoVqS41dQHnFcB5wEtHnM8ExjnntptZT+BLoFWZ65c754664J6ZdQcuAXoALYFpZtbZOae1BkQ8tiUrl2v++QMb9xzkrxf05kI1cQSt9vH1+PDmoTwydRWTZqWzcNNenrm0Hy3j6ngdTSTo1cgRP+fcaufc2nLOL3bObfe9XAnEmJk/c0Djgf845/KdcxuA9cCgE08sIififzt3B6voCwExkeE8MqEXT1/ajzU79jP26dlMX7vb61giQa9GFn6VdD6w2DmXX+bca75p3vut/FVdWwFbyrzeyv+OGIpINTvcuVs3KoIPbx7GkMQmXkeSanRun5Z8fNtwmjWI4eev/cBfvlhDkaZ+RaqMZ4WfmU0zsxXlHOMr8dkewOPAjWVOX+6c6wWM8B1XlPfRcs6Vu2edmd1gZilmlpKRkXHsb0hE/OKc44UZpZ27PVs15KObh6pzN0QlJtRn8i3DuGRgG56fkcZlr8xn1/48r2OJBCXPCj/n3OnOuZ7lHFMq+pyZtQY+Aq50zqWVud823585wFuUP4W7FSg7h9Qa2F7O+3DOTXLOJTvnkhMSEvz75kSkQoXFJdz7wXIe/2IN4/q05N/XDVbnboiLiQznsfN78/eL+rB8azZnPzWb71IzvY4lEnRq1VSvmcUBU4HfOufmlDkfYWbxvq8jgXMobRA50sfAJWYWbWYdgCRgQdUnF5HDsg8VcvVrC3gnZQu3ndaJpy7uq85d+dF5/Vvz8a3DaFwviitenc+TX6+juKTciRkROQ41svAzs4lmthUYAkw1sy99l24FOgH3H7FsSzTwpZktA5YA24CXffc618weBnDOrQTeBVYBXwC3qKNXpPpsycrl/BfmsmBDFk9c2Ie7zuiiPXflJ5KaxTLl1mFM7NeKp75J5cpX55ORk3/sD4rIMZlz+k3qWJKTk11KylFXiRGRSli0eS/X/yuFohLHiz8boCYOOSbnHO+lbOX+KStoUCeSpy/pp/9uRCrBzBY655LLu1YjR/xEJLhMXbaDSyfNo35MBB/ePFQ/vKVSzIyLBrZh8i3DiI2O4PJX5vHst6mUaOpX5Lip8BORKuOc47np67nlrUX0atWQj24eRmKCOnfFP91aNODj24ZzTu+WPPHVOq7+5w/sOaCpX5HjocJPRKpEYXEJv/lgGX/9ci3n9mnJm9cNpnG9KK9jSS1VPzqCpy7py6MTezIvfQ9jn/6OHzZmeR1LpNZR4SciAZd9qJCrXl3Auylbuf20Tjx1iTp35cSZGZcPbseHvxhKdGQYl0yax4sz0zT1K+IHFX4iElCHO3d/2JjF3y7sw6/O6EL5G+mIHJ+erRryyW3DGdOjGY99vobrXk9h78ECr2OJ1Aoq/EQkYBZu2suE5+aQkZPPG9cO5vwBrb2OJEGqQUwkz13Wn4fO7cHs1AzGPj2bRZv3eh1LpMZT4SciAfHpsu1c+vJ/O3dP6qjOXalaZsZVQ9vzwS+GEhZmXPTi97wyOx0tUyZydCr8ROSEHO7cvfWtxfRprc5dqX69W8cx9bYRnNq1KY9MXc2NbywkO7fQ61giNZIKPxE5bgVFJdzzfmnn7vi+6twV7zSsG8mkKwZw39hufLtmN2Ofmc2yrfu8jiVS46jwE5Hjkp1buufuewu3cvuoJP7v4r5ER6hzV7xjZlw3oiPv3jSEkhLHBS98z7/mbtTUr0gZKvxExG+b9+Ry3gtz/tu5O7qzOnelxujfthFTbx/B8KR4/vDxSm59azH78zT1KwIq/ETETws37WXi83PIPFCgzl2psRrVi+KVK5O596yufLFyJ+c+8x0rt2d7HUvEcyr8RKTSPlla2rkbGxPBR+rclRouLMy46eRE/nPDSRwqLGbi83P59/xNmvqVkKbCT0SO6XDn7m1vl3bufnjzMDqqc1dqiYHtG/PZ7SMY3KExv/9oBb98ZwkH84u8jiXiCRV+IlKhgqIS7vZ17k5Q567UUk3qR/Ovnw/irtGd+WTpdsY9+x1rdu73OpZItVPhJyJHlZ1buufu+wu3cseoJJ5U567UYmFhxm2jknjzusHk5BUx4bk5vJuyxetYItVKhZ+IlOtw527Kpiz+flEf7lTnrgSJoYnxTL19OP3aNOKe95dx17tLyS3Q1K+EBhV+IvITCzdlMeH5Oew5WMCb1w7mvP7q3JXg0jQ2hjevG8zto5L4cPFWxj87h9RdOV7HEqlyKvxE5H98vHQ7l748nwYxEXx08zAGq3NXglR4mPGr0Z15/ZpBZB0s4Nxn5/Dhoq1exxKpUir8RAQo7dx99ttUbn97MX1bx/HRzcPoEF/P61giVW5EUgKf3TGCXq0b8qt3l/Kb95eRV1jsdSyRKqHCT0QoKCrh1+8t44mv1jGxXyveuG4QjdS5KyGkWYMY3rpuMDefksg7KVuY+Pxc9uUWeB1LJOBU+ImEuOzcQq58dT4fLNrKL09P4u8X9VHnroSkiPAw7jmzK69enUza7gPc8tYiCotLvI4lElAq/ERC2KY9B5n4whwWbdrHkxf34Zenq3NX5LSuzXh0Yk/mrN/DI5+u8jqOSEBFeB1ARLyRsjGLG95YSIlzvHndYAZ1aOx1JJEa48LkNqzblcPLszfQuXkslw9u53UkkYBQ4ScSgqYs2cbd7y+jVVwdXr16oJo4RMpx71ndSN19gD9MWUnH+PoMSVSHu9R+muoVCSHOOZ75JpU7/rOEvm3i+PAXQ1X0iRxFeJjx9KX9aNekLjf/eyGb9+R6HUnkhKnwEwkRxSWOu99fxt++9nXuXqvOXZFjaRATyStXDaTEwXWv/0BOXqHXkUROiAo/kRDxyux03l+4ldtHqXNXxB8d4uvx/OX9Scs4yJ3vLKG4xHkdSeS4qfATCQGpu3L429frGNOjGXeenqTOXRE/DesUzwPndGfa6t088dVar+OIHDc1d4gEuaLiEn793lLqR0fw6MReKvpEjtOVQ9qxdlcOL8xIo3Oz+kzspz2spfbRiJ9IkHtpVjpLt2bzx/E9ia8f7XUckVrLzHjo3B4M7tCY33ywnMWb93odScRvKvxEgtianfv5v2nrGNu7BWN7t/A6jkitFxkexgs/G0CzBtHc8MZCdmQf8jqSiF9U+IkEqcLiEu56dykN60Tyx/E9vY4jEjQa14vilSsHkptfxA2vL+RQQbHXkUQqTYWfSJB6fnoaK7fv55EJvWisZVtEAqpL81ieuqQfK7Znc88Hy3BOnb5SO6jwEwlCK7dn88y3qYzv25Izezb3Oo5IUDq9ezPuHtOFT5Zu57np672OI1Ip6uoVCTIFRaVTvI3qRfHguB5exxEJar84OZF1O3N44qt1dGoaq1+0pMbTiJ9IkHn221TW7MzhTxN7aWcOkSpmZjx2fm/6tInjV+8uYfWO/V5HEqmQCj+RILJ8azbPzUjjvP6tGN29mddxREJCTGQ4L18xgNiYCK77VwqZB/K9jiRyVCr8RIJEflExd723hPj6UfzhHE3xilSnpg1iePnKZDIP5POLNxdSUFTidSSRcqnwEwkST01LZd2uAzx2fm8a1o30Oo5IyOndOo4nLuzDDxv3cv/kFer0lRpJzR0iQWDx5r28ODONi5PbcGqXpl7HEQlZ4/q0ZN2uHJ75dj1dmsdyzfAOXkcS+R8a8ROp5fIKi/n1e0tp3iCG35/Tzes4IiHvztM7M6ZHMx6ZuoqZ6zK8jiPyP1T4idRyf/96HWkZB3n8gt40iNEUr4jXwsKMv1/Ul87NYrn1rUWkZRzwOpLIj1T4idRiCzdl8fLsdC4b3JYRSQlexxERn3rREbxyVTJR4WFc/68UsnMLvY4kAqjwE6m1DhUU8+v3ltGyYR1+d7ameEVqmtaN6vLiFQPYsjeXW99eRFGxOn3Feyr8RGqpv365lg2ZB/nrBb2pH60+LZGaaGD7xjwyoSezUzN59LPVXscRUVevSG00P30Pr83dwJVD2jG0U7zXcUSkAhcPbMvanQd4dc4GujSL5ZJBbb2OJCHsqIWfmZUAfi9C5JwLP6FEIlKh3IIi7n5/GW0a1eU3Z3b1Oo6IVMLvzu7K+owD3D9lBR0T6jOoQ2OvI0mIqmjE72F+WvhNAHoCXwJrAQO6AGcAy4EpVZBRRMp4/PM1bM7K5Z0bTqKepnhFaoWI8DCeubQfE5+bw01vLmTKLcNo07iu17EkBB31p4Zz7sGyr83sGqA50NM5t/aIa92A6cDmKsgoIj5z12fyr+838fNh7RncsYnXcUTEDw3rRPLKVclMeG4O17+ewvu/GKrnc6Xa+dPccQ/w7JFFH4BzbjXwHPCbQAUTkf91IL90irdDfD3uGaMpXpHaqGNCfZ69rD/rduVw5ztLKCnRtm5Svfwp/NoBhyq4nut7j4hUgT99tprt2Yd44sLe1InSo7QitdXIzgncN7Y7X6/axd+/Xud1HAkx/hR+64AbzCzuyAtm1gi4gdLn/kQkwGaty+Ct+Zu5fkRHBrTTQ+Eitd3Ph7XnkoFteHb6eqYs2eZ1HAkh/jxc8DtgMpBqZm9QWgg6oCvwMyCO0uYPEQmg/XmF3PvBMhIT6vGr0Z29jiMiAWBmPDy+J+kZB7nn/WW0b1KPPm1+Mq4iEnCVHvFzzk0FxlDawPFL4HngBeAO37mzfO8RkQB69NPV7NyfxxMX9iEmUlO8IsEiKiKMF37Wn/j60dzwRgq79ud5HUlCgF87dzjnvnXODQBaAEOAoUAL59wA59y0qggoEsqmr93NOylbuPHkRPq1beR1HBEJsCb1o3nlqmRy8oq44fUU8gqLvY4kQe64tmxzzu1yzs13zs1zzu0KdCgRgezc0inezs3q88vTk7yOIyJVpFuLBjx5cV+Wbs3mNx8swzl1+krV8avwM7M4M3vMzFaY2UEzO+D7+k/lNX2IyPF7+NNVZB4o4IkL+xAdoSlekWA2pkdzfn1GZ6Ys2c4LM9O8jiNBrNKFn5m1AhZTup4fwFTgc0obPO4FFplZy4AnFAlBX6/axQeLtnLzKYn0bq3fqURCwS2ndmJcn5b89cu1fL1Kk2lSNfwZ8fsz0Aw4xznX0zl3kXPuQudcL2Cs79qfqyKkSCjZe7CA3320nK7NY7ntNE3xioQKM+OvF/SmV6uG/PI/i1mzc7/XkSQI+VP4nQk85Zz77MgLzrnPgWeAswIVTCRUPfjJSvYeLOBvF/UhKuK4HsMVkVoqJjKcSVckUy86guv+lcKeA/leR5Ig489PlVigolUmt/reIyLH6YsVO5iyZDu3nZZEj5YNvY4jIh5o3jCGSVcmszsnn1/8exEFRSVeR5Ig4k/htxa4wMx+8hkzCwcuQDt3iBy3PQfy+f1HK+jRsgE3n5rodRwR8VDfNnH89YLeLNiQxR8+XqlOXwkYf3bueBp4BfjWzJ7kv0VeV0oXcR4BXBvYeCKh44GPV7I/r5B/Xz+YyHBN8YqEuvF9W7F2Zw7Pz0ija/NYrhra3utIEgQqXfg55141s6bAH4APy1wyIB/4nXPun4GNJxIaPl22nanLdnD3mC50bd7A6zgiUkP8+owurNt1gIc/XUViQn2GJ8V7HUlqOfN3+NjMmgCjgXa+UxuBr51zWYGNVnMkJye7lJQUr2NIkMrIyeeMJ2fSpnFdPvzFUCI02iciZRzIL+L85+eyI/sQU24dTof4el5HkhrOzBY655LLu+b3Txjn3B7n3H+cc4/7jncCXfSZ2YVmttLMSswsucz50Wa20MyW+/48rcy1GWa21syW+I6m5dy3iZlN9y08/WwgM4scD+cc901ezsH8Yv52YR8VfSLyE/WjI3jlqmTCw4xr//UD2YcKvY4ktZjfP2XM7Ewze9bMpprZp76vzwhwrhXAecCsI85nAuN8awdeBbxxxPXLnXN9fcfucu6bB9wP/DrAeUWOy8dLt/Plyl386ozOJDVTU7yIlK9N47q88LMBbN6Ty+1vL6a4RM0ecnz82bkjysymULpjx83AQGCw7+vPzWyymUUFIpRzbrVz7icdws65xc657b6XK4EYM4v2474HnXPfUVoAinhq9/48Hpiykn5t47h+REev44hIDXdSxyY8PL4nM9dl8OfPVnsdR2opf0b8/gCMA/4GNHXONXXOJQAJwBPAuZSOplWX84HFzrmyq1u+5pvmvd/MrBqziPjFOcfvPlpOXmExT1zYh/Aw/ecqIsd22eC2XD20Pa98t4F3f9jidRyphfwp/C4D3nTO3eOcyzx80vfM32+AN4GfVfZmZjbNzFaUc4yvxGd7AI8DN5Y5fblvCniE77iislmO8nfcYGYpZpaSkZFxIrcS+YkPF21j2urd3D2mC4kJ9b2OIyK1yH1juzG8Uzy/n7yclI1B21cpVcSfwq8lMLeC698DLSp7M+fc6b49f488plT0OTNrDXwEXOmcSytzv22+P3OAt4BBlc1ylHyTnHPJzrnkhISEE7mVyP/YmZ3Hg5+sZGD7Rvx8WAev44hILRMRHsazl/WjVVwdbnpzIVv35nodSWoRfwq/7cBJFVwfDOw4sTgVM7M4Sp8x/K1zbk6Z8xFmFu/7OhI4h9IGEZEaxTnHvR8uo7C4hL9eoCleETk+cXWjeOWqgeQXlXD96ws5mF/kdSSpJfwp/N4GrjCzR8ys0eGTZtbIzP5I6dTqW4EIZWYTzWwrMASYamZf+i7dCnQC7j9i2ZZo4EszWwYsoXRP4Zd99zrXzB4uc++NwN+Bq81sq5l1D0Rmkcp4L2UrM9ZmcO+ZXWmvtbhE5AR0alqfZy7tx9qd+7nr3aWUqNNXKqHSCzj7umc/As4EHHD4wbcESnfv+AKYeESzRVDQAs4SCNv2HeLMJ2fRvWUD3r7+JMI02iciAfDK7HQembqa20cl8avRnb2OIzVARQs4+7NlWz5wtpmdA4wF2vsubQQ+cc59doI5RYKWc457P1hGsXP89YI+KvpEJGCuHd6BtTtzePqbVDo3q885vVt6HUlqsEoXfoc55z4FPq2CLCJB6+0FW5idmskfJ/SkbZO6XscRkSBiZjwysSfpmQf59XtLad+kHj1bNfQ6ltRQ2h9KpIptycrl0amrGNapCZcPaut1HBEJQtER4bz4swE0rhvF9a+nsHu/9imQ8vlV+Pn2yn3HzH4wszQzSz/iSDv2XURCR0mJ4573lwHw+Pm9NcUrIlUmITaal69KZl9uITe8sZC8wmKvI0kN5M+WbXdS2sBxMqVds7OAmUccR+6tKxLS3py/ie/T93DfOd1p3UhTvCJStXq0bMiTF/dhyZZ9/O7D5VS2gVNChz/P+N1JaXF3pnOuoIryiASNTXsO8ufP1jAiKZ5LBrbxOo6IhIgze7bgV6M78/ev19GleSw3npzodSSpQfyZ6o0H3lHRJ3JsJSWOu99bRkSY8fj5vdHW0SJSnW47rRNje7fgsS/WMDtV247Kf/lT+C0COlZVEJFg8s+5G1mwMYsHxnWnZVwdr+OISIgxM564oA8d4uvx+49W6Hk/+ZE/hd8vKd25Y3RVhREJBukZB/jLl2s4rWtTLhjQ2us4IhKi6kSF88j4nmzOyuX56eu9jiM1xFGf8TOzr8o5nQN84dv2bCNw5K8Qzjk3JmDpRGqZ4hLH3e8vIyo8jD+f10tTvCLiqaGd4pnQtyUvzkxnQr9WdEyo73Uk8VhFI36dgaQjjihgs+9zHcu5rr1iJKS9+t0GFm7ay0Pje9CsQYzXcURE+N3YbkRHhnH/lBXq8pWjj/g559pXYw6RWm/97hz++tVaRndvxoS+rbyOIyICQNPYGO4e04UHpqzk46XbGa//P4U07dwhEgBFxSXc9d4y6kaF8+jEnpriFZEa5fLB7ejduiGPTF3N/rxCr+OIh1T4iQTApNnpLN2yj4fH96RprKZ4RaRmCQ8zHpnQk8wD+fz9q3VexxEPHbXwM7MSMysys6gyr4uPcRRVX3SRmmHtzhz+7+tUzurZnHG9W3gdR0SkXL1bx3HFSe14/fuNrNiW7XUc8UhFO3c8DDig6IjXIuJTWFzCr99bSv2YCP44QVO8IlKz3XVGFz5bvpPff7ScD28eRrj2Dw85FTV3PFjRaxGBF2eksXxbNs9f3p/4+tFexxERqVDDOpHcN7Ybv3xnCW8t2MwVJ7XzOpJUMz3jJ3KcVm3fz9PfpjKuT0vO7qUpXhGpHcb3bcnQxCb85Ys1ZOTkex1HqllFCziPPJ4bOs8w01cAACAASURBVOdmHX8ckdqhoKh0irdhnSgePreH13FERCrNzHh4fE/OemoWf/5sNX+/uK/XkaQaVfSM3wz8e6bPfO8PP5FAIrXBc9PXs2rHfiZdMYBG9aK8jiMi4pdOTetz48hEnp2+nguT2zAksYnXkaSaVFT4nVptKURqkRXbsnlu+nom9mvFGT2aex1HROS43HJqJyYv2cb9U1bw2e0jiIrQ01+hoKLmjpnVGUSkNsgvKuaud5fSuF4UfxjX3es4IiLHrU5UOA+P78E1/0zh5dnp3HJqJ68jSTU4rvLezJLMbJiZNQx0IJGa7OlvUlm7K4c/n9eLuLqa4hWR2u20rs0Y06MZz3ybypasXK/jSDXwq/Azs4vNbBOwBpgFDPCdjzezVDO7sAoyitQIS7fs44UZaVwwoDWjujXzOo6ISED8YVwPwsx46JOVXkeRalDpws/MxgNvA5uB+ylt5gDAOZcJrAauCHRAkZogr7CYu95bStPYGO4/R1O8IhI8WsbV4ZenJzFt9W6+WrnT6zhSxfwZ8bsPmOWcGwG8VM71+UCfgKQSqWGenLaO9bsP8Nj5vWhYJ9LrOCIiAfXzYR3o0iyWhz5ZRW6Bdl8NZv4Ufj2Adyu4vhPQ/JcEnYWb9vLyrHQuGdiGU7o09TqOiEjARYaH8cjEnmzbd4invkn1Oo5UIX8KvzwgpoLr7YF9J5RGpIbJKyzm7veW0qJhHX4/tpvXcUREqszA9o25cEBr/jF7A2t35ngdR6qIP4Xfd8Cl5V3wdfdeA3wbiFAiNcXbCzaTnnmQP5/Xi9gYTfGKSHD77dndqB8Twf2TV+CcP3s4SG3hT+H3INDDzKYD5/nOJZvZrcASoAHwx8DGE/FOYXEJr8zeQHK7RozsnOB1HBGRKte4XhT3ntmVBRuz+GDRNq/jSBWodOHnnFsEjAGa89/mjseAp4ECYIxzbnXAE4p45LPlO9i27xA3npzodRQRkWpzUXIb+reN40+frWZfboHXcSTA/FrHzzk32znXDegHXEzp1O9AoKtzbm4V5BPxhHOOF2emk5hQj1Fd1dAhIqEjLMx4ZEIvsg8V8vgXa72OIwHmzzp+rQ9/7Zxb6px7zzn3jnNuofM9CGBmJ1dFSJHqNjs1k9U79nPjyETCwuzYHxARCSLdWzbg6qHteXvBZhZt3ut1HAkgf0b8pplZ/NEumtkZwNQTjyTivUmz0mkaG834fi29jiIi4ok7R3emeYMY7vtoBUXFJV7HkQDxp/ALB74ub39eMxsHfAzMC1QwEa+s2JbNd+szuWZ4B6Ijwr2OIyLiifrRETwwrjurduzn9e83eR1HAsSfwu90oAnwuZnVO3zSzC4APgCmA+cENp5I9XtpVjr1oyO4bHBbr6OIiHjqrJ7NOblzAn//eh279ud5HUcCwJ+u3k3AaKAjMMXMos3sCkr37/0MGO+c038VUqttycpl6rLtXDa4LQ20bp+IhDgz46Fze1BQXMLDn67yOo4EgL9dvWspXdKlP6XTuq9ROtp3vnNOPd9S670yO53wMOPnw9p7HUVEpEZoH1+PW07pxNRlO5i1LsPrOHKC/Cr8oLSjFzgLSATeBC51zhUHOphIdcs6WMA7KVsY37cVLRrW8TqOiEiNcdMpHekQX48Hpqwgr1A/8muzoxZ+ZlZoZgXlHcBsoB5wGZBf5lp+dQUXCbQ3vt9EXmEJN4zs6HUUEZEaJToinD+O78nGPbm8ODPN6zhyAiIquPZvQBv1SUg4VFDMv77fyKiuTencLNbrOCIiNc7wpHjG9WnJ8zPSmNC3Fe3j6x37Q1LjHLXwc85dXY05RDz1/sItZB0s0GifiEgF7hvbjelrdnP/lBW8fs0gzLTAfW3j9zN+IsGmqLiEl2dvoG+bOAZ1aOx1HBGRGqtZgxjuOqMzs1Mzmbp8h9dx5DgcdcTPzEYCOOdmlX19LIffL1JbfLFyJ5uzcvnd2V3126uIyDFccVI73l+4lYc/WcXJnROI1dJXtUpFz/jNAJyZ1fEt1TKDip/5M991bXUgtYZzjpdmptMhvh6juzf3Oo6ISI0XER7GoxN7MfH5OTz5dSoPjOvudSTxQ0WF36kAZdbnO7Xq44hUr+/T97B8WzZ/mtiL8DCN9omIVEbfNnFcNqgt/5y7gfMHtKJHy5/s5io1VEXNHTMrei0SDF6amU58/SjO69/K6ygiIrXKPWO68sWKndw3eQUf3DSUMP3yXCuouUNC1uod+5m5LoOrh7YnJlJPKIiI+KNh3Uh+d3Y3Fm/ex39+2OJ1HKmkipo7HjiO+znn3B9PII9ItZk0K526UeH87KR2XkcREamVzuvfindTtvD4F2sY06MZTepHex1JjqGiZ/wePI77OUCFn9R42/Yd4uOl27lqSHvi6kZ5HUdEpFYyMx6Z0JOznprNnz9fwxMX9vE6khzDUad6nXNhx3FovkxqhVe/2wDAtSM6eJxERKR2S2oWy/UjO/L+wq0s2JDldRw5Bj3jJyEnO7eQtxds5tw+LWkVV8frOCIitd5tp3WiVVwd7pu8nMLiEq/jSAVU+EnIeXP+JnILirl+hLZnExEJhLpRETx4bg/W7TrAP3wzKlIzqfCTkJJXWMxrczYwsnMC3Vs28DqOiEjQGN29Gad3a8ZT01LZtu+Q13HkKFT4SUj5cNE2Mg8UcNNIjfaJiATag+eW7uLx0McrPU4iR6PCT0JGcYnj5dnp9GrVkCGJTbyOIyISdFo3qsvto5L4atUuvlm9y+s4Ug4VfhIyvl61iw2ZB7nx5I6YaYV5EZGqcO3wDiQ1rc8fPl7JoYJir+PIEVT4SUhwzvHizDTaNq7LmT2aex1HRCRoRUWE8ccJPdm69xDPfJvqdRw5QkULOP+PSuzk4YA8YCsw0zm3/USCiQTSDxv3smTLPh4e34OIcP2+IyJSlU7q2ITz+rfi5dnpnNe/FZ2axnodSXwqXfhRupOH83195DzZkeeLzexF4HbnnEPEYy/NTKNR3UguHNDG6ygiIiHhd2d345vVu7lv8grevv4kPWJTQ/gz9NEaWAq8CSQDDX3HIODfwBIgCRgA/Ae4GfhNIMOKHI/UXTl8s2Y3Vw1tT50obS4jIlId4utHc8+ZXZiXnsXkJdu8jiM+/hR+zwBpzrmrnHOLnHM5viPFOXclsAH4k3NusXPuCuAb4OoqyCzil0mz0omJDOPKIe29jiIiElIuHdiWPm3ieHTqarJzC72OI/hX+J1OaTF3NN8AY8q8/hRofxyZRAJmZ3Yek5ds4+LkNjSuF+V1HBGRkBIWZjw6oSdZBwv461drvI4j+Ff4lQB9Krjeh/8+63dYrt+JRALotTkbKC5xXKft2UREPNGzVUOuHNKef8/fzNIt+7yOE/L8KfwmA9eb2b1m9mN7jpnFmtlvget87zlsCLAuMDFF/Lc/r5B/z9/M2b1a0KZxXa/jiIiErLvO6ExC/Wjum7yC4hL1fHrJn8LvTuB74E/AXjPbYWY7gL3Ao8B833swsxhgH/D3wMYVqby35m/mQH4RN45M9DqKiEhIi42J5P5zurN8WzZvztvkdZyQVunCzzm3DxgJnA+8SmmH71LgH8B5wHDfe3DO5TnnfuGce/d4QpnZhWa20sxKzCy5zPnRZrbQzJb7/jytzLUZZrbWzJb4jqbl3Peon5fgkl9UzKvfbWBYpyb0at3Q6zgiIiHvnN4tGJEUzxNfrmX3/jyv44Qsv1aydaU+cs7d4Jw703fc6JybHOD1+lZQWkzOOuJ8JjDOOdcLuAp444jrlzvn+vqO3eXc91iflyAxZcl2dufka7RPRKSGMDMeOrcH+UUlPDJ1tddxQpY/CzgDYGaNgFFAB0qbOTYA3xwe7QsE59xq39915PnFZV6uBGLMLNo5l1/J+57Q56V2KClxTJqVTrcWDRiRFO91HBER8emYUJ+bTknk6W9SuXhgG4Z10v+jq5tfI35m9itKt2R7B3gc+AvwHrDNzO4MfLwKnQ8sPqJoe803zXu/HXuJ8PI+L0Hg2zW7Wb/7ADed3FErxYuI1DA3n5JIuyZ1uX/yCvKLir2OE3IqXfiZ2VXAE5Tu0HEx0BPoBVwELAaeMLMr/bjfNDNbUc4xvhKf7UFp4XljmdOX+6ZwR/iOK/z8/JHvucHMUswsJSMjo7LfltQAL81Ko1VcHc7u1cLrKCIicoSYyHAeHt+T9MyDTJqZ7nWckOPPVO+dwHfAqc65siX6SjP7CJgO/Ap4vTI3c86d7sff/SMzaw18BFzpnEsrc79tvj9zzOwtSreS+0mWo32+nHyTgEkAycnJ6j2vJRZu2ssPG/fywDndiQz3a0BbRESqycmdExjbqwXPTl/P+L6taNtES25VF39+MnYB3j2i6APAd+5d33uqjJnFAVOB3zrn5pQ5H2Fm8b6vI4FzKG0QqdTnJXhMmpVGwzqRXDywjddRRESkAvef052IMOOBj1cQ2P5QqYg/hV8O0LKC66187zlhZjbRzLZSugj0VDP70nfpVqATcP8Ry7ZEA1+a2TJKp6K3AS/77nWumT18jM9LEEjLOMBXq3Zx5ZB21Iv2u29JRESqUfOGMdw5ujMz1mbw5cqdXscJGVbZKtvM3gQmABOcc9OOuDaK0l07JjvnjvpsXW2VnJzsUlJSvI4hx/DbD5fxwaJtzL33NOLrR3sdR0REjqGouIRxz85hX24B0351sn5pDxAzW+icSy7vmj8jfvdSukvHl2a21Mz+4zuWAF/5rv32xOOK+G93Th4fLNzGBQNaq+gTEaklIsLDeGRCT3Zk5/F/07TLa3XwZ+eOrUBf4EkgChjvO6KBvwH9fO8RqXb/nLORwpISrh/R0esoIiLihwHtGnHpoDa8Omcjq3fs9zpO0PN35449zrlfO+e6Oefq+I5uzrl7nHN7qiqkSEUO5BfxxrxNnNmjOR3i63kdR0RE/HTPmK40rBPJfZNXUFKiRo+qpPUupNb7z4LN5OQVccNIjfaJiNRGjepFce9ZXVm4aS/vL9TkYVU66lOUZvbAcdzPOef+eAJ5RPxSWFzCP77bwOAOjenXtpHXcURE5Dhd0L8176Vs4c+fr2Z092Y0qhfldaSgVFH7zIPHcT8HqPCTavPJ0u3syM7jTxN7eR1FREROQFiY8ciEXox9ejaPfb6Gxy/o7XWkoHTUqV7nXNhxHOHVGV5Cm3OOl2am07lZfU7pkuB1HBEROUFdmsdy7fAOvJOyhZSNWV7HCUp6xk9qrRnrMli7K4cbRiZiZl7HERGRALh9VBItG8Zw3+QVFBWXeB0n6Kjwk1rrpZlpNG8Qw7l9KtpQRkREapN60RE8MK4Ha3bm8M+5G72OE3RU+EmttHTLPualZ3Ht8A5EReg/YxGRYDKmRzNO69qUJ79ex47sQ17HCSr6iSm10qRZ6cTGRHDJoDZeRxERkQAzMx46twfFzvHwJ6u8jhNUVPhJrbMx8yCfr9jBz05qR2xMpNdxRESkCrRpXJfbTkvi8xU7mb52t9dxgoYKP6l1XvkunYiwMH4+tL3XUUREpApdN6IDHRPq8YcpK8krLPY6TlA4rsLPzGLMrJWZaXVFqVaZB/J5L2UrE/u1ommDGK/jiIhIFYqOCOeR8T3ZnJXLizPTvI4TFPwq/MxsuJnNBnKAzcBw3/l4M/vGzM6ogowiP3r9+03kF5VwvbZnExEJCUM7xXNWz+b8Y/YGsnMLvY5T61W68DOz4cA3QHPgFeDHhdOcc5m+19cGOqDIYbkFRbz+/UZGd29Gp6b1vY4jIiLV5PZRSeTkF/HqnA1eR6n1/BnxewRYBfQE7i/n+kxgYCBCiZTn3R+2sC+3kJtO1mifiEgo6daiAWN6NOPVORvIPqRRvxPhT+GXDPzTOZdP6Z68R9pK6WigSMAVFZfw8uwNJLdrxIB2jb2OIyIi1ez2UUnk5BXxLy3qfEL8KfxKKL/gO6wlkHticUTKN3X5DrbtO8QNerZPRCQk9WjZkNO7NeMf320gJ0+jfsfLn8LvB+Dc8i74unt/BswNRCiRspxzvDQznY4J9Ti9WzOv44iIiEfuGJVE9qFCjfqdAH8Kvz8Bp5jZ65RO+wK0MbNzgFlAB997RAJqzvo9rNqxnxtHdiQszI79ARERCUq9WjfktK5NeeW7DRzIL/I6Tq1U6cLPOfcNcBlwNvCZ7/SrwMdAZ+Ay59y8gCeUkPfSrDQSYqOZ0K+V11FERMRjd4xKYl9uIa9/v9HrKLVShD9vds69a2afAqMpLfbCgPXAl865A1WQT0Lcim3ZzE7N5DdndiU6ItzrOCIi4rE+beI4pUsCr8zewFVD2lMv2q9SJuT5vXOHcy7XOTfFOfdX59zjzrkPVPRJVZk0K516UeFcNrit11FERKSGuH1UElkHC3hz3iavo9Q62qtXaqwtWblMXb6Dywa3pWGdSK/jiIhIDdG/bSNGJMUzaVY6uQV61s8fRy38zKzEzIr9PPSvLwHzj+82YMA1wzt4HUVERGqYX56exJ6DBbw1f7PXUWqViibGH6bidftEqszegwW888MWxvdtRYuGdbyOIyIiNcyAdo0Z1qkJL85M5/LB7agTpefAK+OohZ9z7sFqzCHyP96Yt4lDhcVasFlERI7qjlGdueil73l7wWbNDlWSnvGTGievsJh/zt3IaV2b0qV5rNdxRESkhhrUoTEndWzMizPTyCss9jpOrXDUET8zGwngnJtV9vWxHH6/yPF6b+FWsg4WaLRPRESO6Y5Rnbn05Xn8Z8Fmrh6mUb9jqegZvxmAM7M6zrmCw68reL/5rmuSXY5bcYnjldnp9GkTx+AOjb2OIyIiNdyQxCYM6tCYF2amccmgtsREqgypSEWF36kAvqIP4DTU7CFV7MuVO9m0J5d7z+yKmbZnExGRY7tjVBKXvzKf91K2cMWQ9l7HqdEqKvz6AF8cfuGcm1HlaSSkOed4aWYa7ZvU5Ywezb2OIyIitcTQxCYkt2vE8zPSuGhgG+30VIGKmjueBJIPv/Ct03dZ1UeSUDUvPYulW7O5fmRHwsM02iciIpVjZtw+Kokd2Xm8v3Cr13FqtIoKv31AfJnX+kksVeqlWWnE14/i/P6tvY4iIiK1zIikePq1jeP56WkUFJV4HafGqmiqdx5wn5m1A7J9584zs04VfMY55/4YsHQSMtbs3M+MtRncNbqzHswVERG/mRl3jEri6td+4MNFW7lkkPZ4L09Fhd8twD+BOygdGXTAeb7jaBygwk/8NmlWOnUiw7liSDuvo4iISC11cucE+rRuyLPT13P+gNZEhmu54iMd9V/EObfROXcKEAO0pnSq9zagTQWHymvx2/Z9h/h4yXYuGdSGuLpRXscREZFaysy44/Qktu49xEeLtnkdp0aqaMQPAOdcEbDdzB4CZjrn9C8pAfXqdxtwwLXabkdERE7QqV2a0qtV6ajfef1bEaFRv/9R6X8N59xDzrkVAGaWYGYDzSzZzBKqLp4Eu+zcQt5esJlxvVvQulFdr+OIiEgtd7jDd3NWLpOXbPc6To3jVxlsZkPMbB6wk9Lmj/nATjOba2YnVUVACW5vzt/EwYJibhiZ6HUUEREJEqd3a0r3Fg14bvp6iorV4VtWpQs/X2H3LdAVeIHS5/1u933dHZhuZoOrIqQEp7zCYl6bs5ERSfF0b9nA6zgiIhIkDo/6bcg8yCfLNOpX1jGf8SvjEWAHMNQ5t7PsBTN7BJjre8/owMWTYDZ58TYyD+Rz08l9vY4iIiJB5ozuzejaPJZnvl3PuX1aaWMAH3+megcDLx1Z9AH4zk3yvUfkmEpKHJNmpdOzVQOGJjbxOo6IiASZsLDSUb/0jIN8qlG/H/lT+DnfcTSaRJdK+3r1LtIzD3LjyETM9FuYiIgE3pk9mtOlWemoX3FJRSVM6PCn8PsBuNHM4o+84Dt3I7AgUMEkeDnneHFmGm0a1+Gsns29jiMiIkEqLMy4bVQn1u8+wOcrdngdp0bw5xm/B4BvgLVm9jqw1ne+K3AFUNf3p0iFUjbtZfHmfTx0bg+tryQiIlXqrJ4t6NQ0lae/SeXsni0IC/Fn/fxZx28OcAawgdJt3J73HbcDacAZzrm5VRFSgstLM9NoVDeSC5Nbex1FRESCXHiYcdtpnVi36wBfrPxJm0LI8Wu4xTk3yzmXDLQAhviOFs65Qc652VURUILL+t05TFu9myuHtKdulD8DziIiIsfnnN4t6ZhQj6e/SaUkxJ/1q1ThZ2Z1zWyhmd0E4Jzb5Zyb7zt2VW1ECSaTZqUTExnGlUPaeR1FRERCxOFRvzU7c/hqVWiXLZUq/JxzuUBH1LkrJ2DX/jw+WryNi5Lb0KR+tNdxREQkhIzr3ZIO8aWjfs6F7qifP1O904FTqiiHhIBX52yguMRx3fCOXkcREZEQExEexi2ndmLVjv1MW73b6zie8afwux3oY2ZPmllnM9MDWlJpOXmFvDVvM2f1akHbJnW9jiMiIiFoQt+WtGtSl6e+WReyo37+FH4bgM6UFoCrgXwzKzjiyK+SlFLrvb1gMzn5Rdw4UqN9IiLijcOjfiu27Wf62tAc9fNn1O7fVLxzh0i5CopK+Md3Gxia2ITereO8jiMiIiFsYr9WPPNtKk9NS+XULk1DbveoShd+zrmrqzCHBLEpS7axa38+f7mgj9dRREQkxEWGh3HLKZ2498PlzFiXwaldmnodqVpp2wSpUiUljkmz0unaPJaRST/Z7U9ERKTande/Na3i6vDUtNDr8PWr8DOzODN71MyWmlm271jqO9eoqkJK7TV97W5Sdx/gppMTQ244XUREaqaoiDBuPjWRJVv2MTs10+s41arShZ+ZJQLLgN8C4cA0SvfuDfedW+Z7j8iPXpqZTsuGMYzt3cLrKCIiIj+6YEBrWjaM4akQW9fPnxG/p4FGlO7J29M5d75z7jznXE9gDBDne48IAIs272XBxiyuHdGRyHA9VSAiIjVHdEQ4vzglkYWb9jI3bY/XcaqNPz+NTwb+zzk37cgLzrmvKS36Tg5UMKn9Js1Mp2GdSC4Z2MbrKCIiIj9x0cA2NG8QE1LP+vlT+B0AMiq4vsv3HhHSMw7w5aqdXHFSO+pFa61vERGpeaIjwrnp5I4s2JjFvPQsr+NUC38Kv7eBn5lZ1JEXzCwGuAJ4K1DBpHZ7efYGIsPDuGpoe6+jiIiIHNUlg9rSNDaap75Z53WUauFP4fcJEAksMbPbzexMMxtjZncACylt8vjEzIaWPaoitNRsGTn5fLBoKxcMaE1CbLTXcURERI4qJjKcm05OZF56FvPTg/9ZP3/m4Mo+2/d//HcXDzvKe8z3nvDjiya11ZvzNlFYXML1I7Q9m4iI1HyXDW7L8zPSePrbVP7dsYnXcaqUP4Xfz6sshQQN5xxTlmxjaGITOsTX8zqOiIjIMZWO+nXkkamrSdmYRXL7xl5HqjL+bNn2r6oMUpaZXQg8CHQDBjnnUnznRwOPAVFAAXC3c+5b37UZQAvgkO82Zzjndh9x30HApMMvgQedcx9V6TcTYlZs28/GPbn84hQt6SgiIrXHZYPb8sKMNJ76JpU3rh3sdZwqU1MXV1sBnAfMOuJ8JjDOOdcLuAp444jrlzvn+vqO3fzUCiDZOdcXOBN4yczUchpAnyzbTmS4MaZHc6+jiIiIVFrdqAhuGNmR2amZLNq81+s4VaZGFn7OudXOubXlnF/snNvue7kSiDGzSncPOOdynXNFvpcx/Pc5RQmAkhLHp0u3MyIpgbi6P2n+FhERqdF+dlI7GteL4qlpqV5HqTI1svCrpPOBxc65/DLnXjOzJWZ2vx1lY1gzG2xmK4HlwE1lCkE5QYu37GV7dh7j+mh7NhERqX3qRUdw3YgOzFyXwZIt+7yOUyU8K/zMbJqZrSjnGF+Jz/YAHgduLHP6ct8U8AjfcUV5n3XOzXfO9QAGAr/1rUFY3t9xg5mlmFlKRkZF61bLYZ8s3UF0RBind2vmdRQREZHjcuWQ9sTVjeTpb4Jz1M+zws85d7pvz98jjykVfc7MWgMfAVc659LK3G+b788cSheSHnSMv381cBDoeZTrk5xzyc655ISEBP++uRBUXOL4dNkOTuvalNiYSK/jiIiIHJf60RFcN7wD367ZzfKt2V7HCbhaNdVrZnHAVOC3zrk5Zc5HmFm87+tI4BxKGzmO/HyHw80cZtYO6AJsrIboQW9++h4yD+RzTu+WXkcRERE5IVcNbU/DOpE8FYSjfpUu/Mysh5mdd8S5U83sGzNbaGZ3BSqUmU00s63AEGCqmX3pu3Qr0Am43/cs3xIzawpEA1+a2TJgCbANeNl3r3PN7GHf54cDS81sCaWjhjc75zIDlTuUfbJsB3Wjwjmta1Ovo4iIiJyQ2JhIrh3egWmrd7FiW3CN+plzlWtsNbOpvvef7XvdClgD5AEZlI6e/dw593oVZfVMcnKyS0lJ8TpGjVVYXMLAR6dxcucEnrqkn9dxRERETlj2oUKGP/4tQxOb8NIVyV7H8YuZLXTOlRvan6ne/sDMMq8vp3Q7tr7Oue7A58Atx51Saq3v1meyL7eQcZrmFRGRINGwTiTXDOvAlyt3sXrHfq/jBIw/hV8jYFeZ12cBMw83VQCfAJ0DFUxqj0+Wbic2JoIRneO9jiIiIhIw1wzrQGx0BM98GzzP+vlT+GUCrQHMrB6lz99NK3M9Av/2/pUgkFdYzNcrd3Fmj+ZER4R7HUdERCRgGtaN5Oph7fls+U7W7szxOk5A+FP4zQZ+4Wvw+D8gEv6/vTuPk6sq8z/++WbtJCSE7AuEEEICJCSCYRcIi6wJ6qijjsA4yqCMoOg48xsGcdxm0XHGBUVlGFFA1BEZmU5YFNkXEZB0JyEESEgg6c4G2UhCtn5+f9zbUBbVa3X3reX7fr3qVV3nnnvrObcqt56cc8+9/F/O8qkkkyqsijzw3Hq27tzD3JkeWmTSigAAIABJREFU5jUzs8rzsXccxKB+vflOhfT6dSTx+0dgB3Ar8DHg6xHxPICk3sD7+NNzAK0K1NY1MGxQP044eHjWoZiZmXW5oQP78ZcnTOSOhY08v7b8e/3anfhFxIvAocDbgIMi4sqcxQOBS4F/7drwrJRt37WH3y1ZxznTx9Cnd1ldEtLMzKzdLj5pEgP69uaae1/IOpSidejXOiL2RER9RKzMK98aEbdHxIoujc5K2u+WrGPH7r0e5jUzs4o2bFA/Ljp+IrX1Dbyw7rWswylKhxI/ScMkfUXSI5Kel3R8Wj5c0hckHdo9YVopqq1rYPSQ/hw9cVjWoZiZmXWrvz7pIGr69OZ795V3r19H7txxAMldMf4eGAxMAgYARMQrwIfwdfyqxpbXd3P/0vWce8RYevdS1uGYmZl1q+H79OfC4w/k9gWrWb6+fHv9OtLj93WghuQcv9OA/F/729NyqwK/XbyWXXubPMxrZmZV469PmkS/Pr343n3Lsg6l0zqS+L0T+E5ELAEK3eftReCALonKSl5tfQPjhw7gyAOGZh2KmZlZjxg5uD8fPvZAfr1gNStf2ZZ1OJ3SkcRvELCuleX7FBmLlYlXt+3i4ec3MHfmOCQP85qZWfX4+MmT6NNLfLdMZ/h2JPFbChzXyvJzgUXFhWPl4K5Fa9jTFMyZMTbrUMzMzHrUqCE1fOiYCdz29GpefnV71uF0WEcSvx8CF0j6KNB8b66QNFjSN4HZwLVdHJ+VoHn1DUwaMYhp44ZkHYqZmVmPu3T2wfTupbKc4duRCzh/H7gOuB5YnhbfCmwEPg1cExE3d3mEVlLWbXmdx5a/whwP85qZWZUaPaSGDx59ALc+tYpVG8ur16+jF3C+DDiRJPm7E/gD8H3gpIi4ouvDs1Jzx8JGImCuh3nNzKyKXTr7YHpJXHt/ec3w7dPRFSLiMeCxbojFykBtfSOHjhnMIaMHZx2KmZlZZsbuO4A/P3p/fvHEy3zy1MmMHzog65DaxTdYtXZbvWkHT63c6Gv3mZmZAZfOngzAD8qo16+jt2y7WNKjktZI2ilpV95jZ3cFatmbX98A4Nm8ZmZmwPihA3jf2w/gF0+8TOPmHVmH0y7tHuqV9A3gM0ADyVDvpu4KykpTbV0jM/fflwOHD8o6FDMzs5LwN7MP5pdPvswPH1jOF8+flnU4berIOX4fJZnQ8a6I2NtN8ViJWrFhGwtXb+aqcw/LOhQzM7OSccCwgbz3qP255Q8vcensgxk9pCbrkFrVkaFeAbVO+qrTvHSY9zwP85qZmf2JT546mb1NwQ8eKP1z/TqS+N0NHN1dgVhpq61r5OiJ+zGuTGYtmZmZ9ZQJwwfyniPHc8vjL7Fu6+tZh9OqjiR+nwKOlvRPksZ3V0BWepau2crStVs9m9fMzKwFl506md17m7jugeVtV85QR+7csQ64CfgC8JKk3Z7VWx3m1TfQS3DOdA/zmpmZFTJxxCDe/bbx3Pz4StZvLd10qCOzer8MXEUyq/dJPKu3KkQE8+obOf7g4Ywc3D/rcMzMzErWZadN5tcLVnP9Q8u5skQnQ3ZkVu/H8azeqrO4YQsvbtjGx0+elHUoZmZmJW3SyH04f+Y4bnxsJZecPInh+5Reh0lHzvGrwbN6q05tXQN9eomzp4/JOhQzM7OSd9lpk3l9z16uf/jFrEMpqCOJ32/wrN6q0jzMe9IhIxg6sF/W4ZiZmZW8yaMGM2fGOG58dAUbt+3KOpy36Ejidxnwdklf8qze6vDHlzaxetMOz+Y1MzPrgMtPm8z23Xu5/uHSm+HbkcRvFTAd+Dye1VsVausa6NenF+88fHTWoZiZmZWNKaMHc+70sfzk0ZVs2l5avX4dmdzxUyC6KxArLXubgvkLGzl16kgG1/TNOhwzM7Oycvnpk5m/sJEfPfwinz1zatbhvKHdiV9EfKQb47AS8/iLr7B+604P85qZmXXCoWOGcPa0MdzwyAo+dtIk9h1QGp0oHRnqtSoyr76Rgf16c9qho7IOxczMrCx96vRD2LpzDzc8UjozfFvs8ZN0MkBEPJj7ui3N9a187d7bxJ0LGznjsNEM7NeRswHMzMys2eHjhnDm4aP50cMv8tF3HMSQEjh1qrVf9fuBkDQgInY1v26lvtLlvbssOsvEIy9sYOP23cyZ4Vu0mZmZFeNTpx/Cb55Zy08eWcHlpx+SdTitJn6nAqRJ3xuvrfLNq29kcE0fTpk6MutQzMzMytr08ftyxmGjuP7hF/nIiRMznzDZYuIXEQ+09toq0849e7l70RrOmj6G/n3ceWtmZlasT51+COd/9xFufGwlnzx1cqaxtHtyh6R7JZ3eyvJTJd3bNWFZVh5Yup6tO/d4Nq+ZmVkXmbH/UE6dOpLrH1rOtp17Mo2lI7N6ZwOtXcl3FHBKUdFY5mrrG9lvYF9OOHh41qGYmZlVjE+dfggbt+/mxsdWZhpHRy/n0trkjoOB14qIxTK2fdce7nlmLeccMZa+vX2lHzMzs65y5IT9OHnKSBY1bM40jlav1SHpQuDCnKIrJf1VgapDgSOB33ZhbNbD7n12HTt272XuDA/zmpmZdbUfXHBU5pdJa+vdhwHNc48DGAMMzqsTwDaSW7p9vkujsx5VW9fAqMH9OeagYVmHYmZmVnGyTvqgjcQvIr4NfBtAUhNwRUTc0hOBWc/a+vpu7lu6nr84ZgK9eynrcMzMzKwbdORevT7pq4L99pm17NrT5Nm8ZmZmFczJnAHJMO/4oQM4asLQrEMxMzOzbuLEz9i4bRcPPb+BOTPGInmY18zMrFI58TPuWryGPU3hYV4zM7MK58TPmFffwEEjBjFt3JCsQzEzM7Nu5MSvyq3b+jqPLXuFuR7mNTMzq3hO/KrcnQvX0BQwx8O8ZmZmFa9diZ+koZKOkXRwK3UOknRR14VmPaG2roGpowczZXT+dbnNzMys0rSZ+En6Z2At8BjwnKQ/SDqqQNUTgBu6OD7rRg2bdvDkyo3MnTk261DMzMysB7Sa+En6c+BK4AHgk8C/ABOARyVd0P3hWXeaX98IwBzfm9fMzKwqtHXnjiuA+yLizOYCSf8J/AL4saQREfGt7gzQuk9tfQNHjN+XiSMGZR2KmZmZ9YC2hnoPBX6VWxARG4GzgRuB/5D0lW6KzbrRig3bqF+12cO8ZmZmVaStHr+mQoUR0QR8VNJG4CpJ+wGPd3Vw1n3mL0yGec/zMK+ZmVnVaCvxew54B3BtoYUR8beStgD/BMzp4tisG9XWNTDrwP0YP3RA1qGYmZlZD2lrqPdO4F2ShrdUISK+BHwGOKArA7Pu89zarTy7ZitzZniY18zMrJq01eP3I+BVYBTwSkuVIuLbklYCM7swNusm8+oa6CU414mfmZlZVWk18YuI1cD32rOhiPg18OuuCMq6T0Qwr76R4yYNZ9TgmqzDMTMzsx7U6Vu2SaqRdJGk0V0ZkHWvxQ1bWL5hG3N9izYzM7OqU8y9evcluVPHtC6KxXpAbX0DfXqJs6eNyToUMzMz62HFJH4A6pIorEdEBPPqGnnHISPYb1C/rMMxMzOzHlZs4hddEoX1iKdf3sTqTTuY62v3mZmZVSX3+FWR2roG+vXpxTun+bRMMzOzalRM4rceOAh4pItieYOk90taLKlJ0qyc8ndKekrSwvT5tJxl90taKmlB+hjVyvYnSHpN0ue6OvZStbcpmF/fyOwpIxlS0zfrcMzMzCwDbV3Hr0XpbdtWdmEsuRYBfwb8MK98AzA3IhokTQfuBsbnLP9wRDzZju1/k+Ti1FXjiRWvsm7rTs/mNTMzq2KdTvy6U0QsAZCUX/50zsvFQI2k/hGxs73blvRuYDmwrQtCLRu1dQ0M6Nub0w9rsSPUzMzMKlyx5/hl6b3A03lJ3w3pMO/Vys8aAUmDgP8HfKmngiwFu/c2ceeiNZx+2CgG9ivJXN/MzMx6QGZZgKR7gEIXk7sqIm5vY91pwNeAM3OKPxwRqyUNBn4FXAjcmLfql4BvRsRrBfLC/Pe4BLgEYMKECa3WLXWPLnuFV7ft8jCvmZlZlcss8YuIMzqznqT9gf8FLoqIZTnbW50+b5V0C3AMb038jgXeJ+nrwFCgSdLrEfHdAvFdB1wHMGvWrLK+bM28ugYG9+/DKVNGZh2KmZmZZaisxv0kDQXmA1dGxCM55X2AoRGxQVJfYA5wT/76EXFSzjpfBF4rlPRVkp179nLX4jWcOW0MNX17Zx2OmZmZZagkz/GT9B5Jq4DjgfmS7k4XXQZMBq7Ou2xLf+BuSfXAAmA18F/pts6X9OWeb0VpePC5DWx9fQ9zZo7NOhQzMzPLmCLKehSzR8yaNSuefLI9V4kpPZ/62dM8+Px6nrjqDPr2Lsk838zMzLqQpKciYlahZc4EKtiOXXu5Z8lazpk+1kmfmZmZOfGrZPc+u47tu/Yy18O8ZmZmhhO/ilZb18DIwf059qDhWYdiZmZmJcCJX4Xa+vpu7l26jvOOGEvvXq1fs9DMzMyqgxO/CnXPkrXs2tPkYV4zMzN7gxO/ClVb18j4oQM48oD9sg7FzMzMSoQTvwq0afsuHnxuPefNGEsvD/OamZlZyolfBbpr0Rr2NAVzZ/jevGZmZvYmJ34VaF59IxOHD2T6+CFZh2JmZmYlxIlfhVm/dSePLtvA3JnjkDzMa2ZmZm9y4ldh7lzUSFPAHA/zmpmZWR4nfhWmtq6BKaP3YeqYwVmHYmZmZiXGiV8Fady8gydWbPSkDjMzMyvIiV8FmV/fCMCcmU78zMzM7K2c+FWQ2roGpo8fwkEjBmUdipmZmZUgJ34VYuUr26hbtdnDvGZmZtYiJ34VYl46zHveDN+b18zMzApz4lchausaOGrCUPbfb2DWoZiZmVmJcuJXAZ5fu5Vn12xlrid1mJmZWSuc+FWA2vpGJDjvCA/zmpmZWcuc+JW5iGBefQPHHTScUUNqsg7HzMzMSpgTvzL3TOMWlq/fxpyZ7u0zMzOz1jnxK3O1dY307iXOme7Ez8zMzFrnxK+MNQ/zvmPyCIYN6pd1OGZmZlbinPiVsQUvb2LVxh2ezWtmZmbt4sSvjNXWNdKvdy/OnDY661DMzMysDDjxK1NNTcH8hQ2cMnUkQ2r6Zh2OmZmZlQEnfmXqiRWvsnbLTg/zmpmZWbs58StTtfUNDOjbmzMOG5V1KGZmZlYmnPiVoT17m7hj4RpOO2wUA/v1yTocMzMzKxNO/MrQo8te4dVtu5g7w8O8ZmZm1n5O/MrQvPoG9unfh9lTR2YdipmZmZURJ35lZueevdy1aA1nThtNTd/eWYdjZmZmZcSJX5l56LkNbHl9j4d5zczMrMOc+JWZ2voGhg7sy4mTR2QdipmZmZUZJ35lZMeuvdzzzFrOmT6Gfn380ZmZmVnHOHsoI/ctXce2XXuZ42FeMzMz6wQnfmWktq6BEfv057hJw7MOxczMzMqQE78y8drOPdz77DrOO2IMvXsp63DMzMysDDnxKxP3PLOWnXuafG9eMzMz6zQnfmWitq6BsfvWcNSE/bIOxczMzMqUE78ysGn7Lh58fj1zZoyll4d5zczMrJOc+JWBuxevYffe8DCvmZmZFcWJXxmYV9/IgcMHcsT4fbMOxczMzMqYE78St+G1nTzywgbmzBiL5GFeMzMz6zwnfiXuzoWNNAUe5jUzM7OiOfErcbV1jRwyah+mjh6cdShmZmZW5pz4lbDGzTt4YuWrzJ05zsO8ZmZmVjQnfiVsfn0jETBnxtisQzEzM7MK4MSvhNXWNzJt3BAmjdwn61DMzMysAjjxK1EvvbKdupc3eVKHmZmZdRknfiVq3sIGAM47wsO8ZmZm1jWc+JWo2rpGjpwwlAOGDcw6FDMzM6sQTvxK0AvrXmNJ4xbmzvAwr5mZmXUdJ34laF59AxKc59m8ZmZm1oWc+JWYiKC2roFjJg5j9JCarMMxMzOzCuLEr8QsadzKsvXbPJvXzMzMupwTvxJTW99A717inOljsg7FzMzMKowTvxISEcyrb+DEySMYvk//rMMxMzOzCuPEr4TUrdrMy6/u8C3azMzMrFs48SshtXUN9O0tzprmYV4zMzPrek78SkRTUzC/vpFTpoxi3wF9sw7HzMzMKpATvxLx5MqNrNnyOnNnepjXzMzMuocTvxJRW9dATd9enHHY6KxDMTMzswrlxK8E7NnbxB0LGzn90NEM6t8n63DMzMysQjnxKwEbXtvFpJGDfNFmMzMz61YlmfhJer+kxZKaJM3KKX+npKckLUyfT8tZdr+kpZIWpI9RBbY7UdKOnDo/6Kk2tWbMvjX88hMncLYv2mxmZmbdqFTHFRcBfwb8MK98AzA3IhokTQfuBsbnLP9wRDzZxraXRcTbui5UMzMzs/JQkolfRCwBkJRf/nTOy8VAjaT+EbGzB8MzMzMzK0slOdTbTu8Fns5L+m5Ih3CvVn7W+KaDJD0t6QFJJ/VAnGZmZmYlIbMeP0n3AIVOarsqIm5vY91pwNeAM3OKPxwRqyUNBn4FXAjcmLdqIzAhIl6R9Hbg15KmRcSWAu9xCXAJwIQJE9rbLDMzM7OSlVniFxFndGY9SfsD/wtcFBHLcra3On3eKukW4BjyEr+0d3Bn+vdTkpYBU4C3nBcYEdcB1wHMmjUrOhOrmZmZWSkpq6FeSUOB+cCVEfFITnkfSSPSv/sCc0gmiOSvP1JS7/TvScAhwPKeiN3MzMwsayWZ+El6j6RVwPHAfEl3p4suAyYDV+ddtqU/cLekemABsBr4r3Rb50v6crr+yUC9pDrgVuATEfFqz7XMzMzMLDuK8ChmW2bNmhVPPtnWVWLMzMzMsifpqYiYVWhZSfb4mZmZmVnXc+JnZmZmViWc+JmZmZlVCSd+ZmZmZlXCiZ+ZmZlZlXDiZ2ZmZlYlnPiZmZmZVQknfmZmZmZVwomfmZmZWZVw4mdmZmZWJXzLtnaQtB5Y2UqVEcCGHgqn1Ljt1aua2++2V6dqbjtUd/vLre0HRsTIQguc+HUBSU+2dE+8Sue2V2fbobrb77a77dWomttfSW33UK+ZmZlZlXDiZ2ZmZlYlnPh1jeuyDiBDbnv1qub2u+3VqZrbDtXd/oppu8/xMzMzM6sS7vEzMzMzqxJO/Foh6WxJSyW9IOkfCiyfLWmzpAXp4ws5y4ZKulXSs5KWSDq+Z6MvXpHt/4ykxZIWSfqZpJqejb44bbU9rTM7bfdiSQ90ZN1S1tm2SzpA0n3p932xpE/3bOTFK+ZzT5f1lvS0pHk9E3HXKvJ7X9bHvCLbXtHHO0l/l3OcXyRpr6Rh7Vm31HW27WV9vIsIPwo8gN7AMmAS0A+oAw7PqzMbmNfC+j8BLk7/7gcMzbpNPdV+YDzwIjAgff0/wEeyblMXt30o8AwwIX09qr3rlvKjyLaPBY5K/x4MPFctbc9Z/lnglpaOC6X8KLb95XzMK/J7X/HHu7z6c4F7O7NuqT2KbHvZHu/c49eyY4AXImJ5ROwCfg68qz0rShoCnAz8N0BE7IqITd0WaffodPtTfYABkvoAA4GGboixu7Sn7X8B3BYRLwFExLoOrFvKOt32iGiMiD+mf28FlpD8KJaLYj53JO0PnAdc30PxdrVOt78CjnlFffZU/vEu14eAn3Vy3VLT6baX8/HOiV/LxgMv57xeReEP9XhJdZLulDQtLZsErAduSId9rpc0qJvj7Wqdbn9ErAa+AbwENAKbI+I33R1wF2pP26cA+0m6X9JTki7qwLqlrJi2v0HSROBI4PFuirM7FNv2bwF/DzR1b5jdppj2l/sxr9Ntr5LjHQCSBgJnA7/q6Lolqpi25y6bSBkd75z4tUwFyvKnQP+R5LYoM4FrgF+n5X2Ao4DvR8SRwDag3M596HT7Je1H8r+mg4BxwCBJF3RjrF2tPW3vA7ydpIfnLOBqSVPauW4pK6btyQakfUgOjldExJbuCrQbdLrtkuYA6yLiqW6OsTsV89mX+zGvmM++Go53zeYCj0TEq51YtxQV0/ZkA2V4vHPi17JVwAE5r/cnr/s+IrZExGvp33cAfSWNSNddFRHN2f+tJAfFclJM+88AXoyI9RGxG7gNOKFnwu4SbbY9rXNXRGyLiA3Ag8DMdq5byoppO5L6khwEfxoRt/VAvF2pmLafCJwvaQXJcNFpkm7u/pC7VLHf+3I+5hXT9mo43jX7IG8O83Z03VJUTNvL93iX9UmGpfog+d/dcpL/xTWf9Dktr84Y3rwW4jEkXf3Nrx8CpqZ/fxH496zb1FPtB44FFpOc6yKSk74vz7pNXdz2w4DfpXUHAouA6e1Zt5QfRbZdwI3At7JuR0+3Pa/ObMpzckdR7S/nY16R3/uKP96l9fYFXgUGdXTdUn0U2fayPd71wQqKiD2SLgPuJpn586OIWCzpE+nyHwDvAy6VtAfYAXww0m8EcDnwU0n9SL5Yf9XjjShCke1/XNKtJEPBe4CnKaOrnren7RGxRNJdQD3JOV3XR8QigELrZtKQTiim7ZLeAVwILJS0IN3kP0bSG1zyiv3cy10XtL9sj3ld8G++oo93adX3AL+JiG1trduzLei8YtpO0stflsc737nDzMzMrEr4HD8zMzOzKuHEz8zMzKxKOPEzMzMzqxJO/MzMzMyqhBM/MzMzsyrhxM/MiibpYkmR3q+2akg6WNIdkjam7b8s65jaIulhSfdkHUdHSOqT7t/PZx2LWblz4mdWQSTdJmm3pJGt1Pl0+iM6tydjq1A/AWaRXLD4QqCc7tFqZlXIiZ9ZZbmJ5Gr0H2ylzgXABuCuLnzfG4ABEbGqC7dZ0iTVkNya68aI+HZE3BwRz2Udl5lZa5z4mVWW+SS3Fip4k3hJU0h6qH4eyX1FiyJpEEBE7I2I14vdXpkZRXLbpk2dWVmJgV0bkplZ65z4mVWQiNgF/A9wjKRDClS5MH2+qblA0rsl1UpaJWlX+vw9SUNyV5T01XSI+HBJP5a0AViRLnvLOX6STpH0C0krJe2UtFbSjZLG5W23ed2TJH1N0hpJOyTdLWlifgMkTZZ0c1pvp6QXJf2wOQlN6wyR9A1JK9I6L6XbrmnPfkxjuU/Sa5K2SvqtpGNz9wWwMn35lTT+Pa1sr/kctR9Ier+kOuB14KPp8o9KuienTcsl/XN6+7Pc7dws6XVJYyTdmsa2MW1/TV5dSbpS0svp/nw0tw15dYensTWm779E0mclqYU2nC+pLt3uAkknp3XOk/R0GuMSSWe2c3+/X9IfJG2WtE3S85K+2471DpR0i6QN6XsukHRRXp3Jadz/IOlvJC1L6/5R0hkFtlnUd8es1PlevWaV5ybgE8CHSc49y/UXwHMR8Yecso8Bu4BrgI3AkcDFJDegP6XA9v8HeAn4AjCowPJmHwCGkty3dB0wFbgEOFbSzAI9hN8kuefzvwAjgc+lbTmpuYKk6cBDJP9pvQ54HtgfeC+wH7BN0gDgfpIbr18HLAPeBnwWmAbMaSVmJJ1Kcu/Ol4Gvpu/1CeABSbMj4vfAL0mGy78J3ArcTnL/1racBLwf+B5wLbAkLb8cWATcCWwjuQ/oP6Rt+8u8bfRK46sH/g44jmS/rgOuzqn3ZeDzwG+B/yPZ/3eS9FAuz2lvDcn+Ogz4PvAscB7wH8AE4Iq89z8emJvGvwP4e2C+pL9O17kW2J6W3yrpgIjY3NIOSZPDXwD3Av9Icr/bScBZLa2TrjcKeJTkO3YN0EByisNPJA2LiG/lrfJBkl7aa0m+7x9P454dEY+l2yzqu2NWFiLCDz/8qLAH8ALwfF7ZiUAAn88rH1hg/Y+kdY/NKftqWvarAvUvTpft38Z2Z6f1PlBg3QeAXjnln0vLp+aUPUiSVEwpsO3me49fTZKQHJ63/G/S7Z3axr5bALwCjMwpGw9sBR7NKZtYaH+2sM0+ad0m4G0FlhfaV18E9gJjc8puTrfz1by684DGnNejSJKb3+bt0+Z9cE9O2RVp2cdy9yVwWxrv1Lw27Mrd/yRJYJD0YE4uUH5xG/vmGpLTE3q3Y/99PqfsW2nZ6Tll/YDHSZLnoWnZ5Jy4D86pOxrYAjyUU1bUd8cPP8rh4aFes8p0MzBZ0nE5ZReQ/Hj9NLdiRGyHN4YG95U0Ang4Xfz2Atv+fnsCaN5uuu3B6XYXkSRQhbb7w4jI7TV7IH2elG5jNEmP2Y+jwCSKiIj0zw8AjwDrJI1ofgDNlzA5raWYlQxVzwR+EhHrc7a9GrgFOF7S8Faa3ZZHI2JBgdibP4Nekoam8T5I0rt3ZIHtXJv3+gFgjN48Z/AsoC/wnbx9ej3J/s81h6T38sc58QTw7yQJ4Hl59e/L2/+PNccQES8UKJ9UIP5cm4DBQLuGhXPMAZ6OiN/lxL2LJCEcyFs/59qIWJZTdy3wM+BESUPT4k5/d8zKhRM/s8p0c/p8AUB6rtifAw9HxIu5FZWcs1cLvEbyI7yeZAgVkmG0fMsKlL2FpP0l/VTSZpKelfXpY3AL212Z93pj+jwsfZ6cPi9s462nAKfnvF/zY2m6fFQr605Mn58tsOyZvDqdUXDfSTpB0n0kvZkbSeJtTmjy99XuiGjIK2veV/ulzwemz0tzK6WJ0Yq8dSeS9A7vzStvbu9BeeUv5b3e1Eb5frTuu2mcd6TnON4i6QOS2joV6UBa/5zy416aXzEtE2/ur2K+O2Zlwef4mVWgiHhB0mPAByRdAZxLkkDdlFsv7el4gGRo7GqSIeLtJENm8yn8n8Mdbb2/pN4kvSQjgK+R/BhvI+lx/GUL281PPN7YXN5ztFAvt/69wL+2sHx1G+u3tt32vH9r3rLvJB1MkuQ9T3Iu2Uskw6YTgP/mrfuqtXMJ27OvVKCsNfnbaOlzauvzK7zxiLWSjiRJuM4i6fn7EPCkpJMjos3vWwvvlx93e/ZFd313zEqGEz+zynUTyZDgWSQ9fztJkq5cp5MkZ++OiEcmsD/wAAADeklEQVSaCyUdXuR7v41kMsEFEfHG0LKSmbf7dnKbzcOIM9qotwwYHBGduTvFivT50ALLmsvyeyaL9W6gBjgnHVIGQNK5RWxzRfp8KG/23jb3/B4IrM2rO0NS77xev8PyttVtIrm00F3pA0mXA98B/oy8UxNyrKT1z2lFC+W5ppAkhM2faTHfHbOy4KFes8r1C5IT2j9Jcj5UbUTkX3Ouufco/1jwd0W+d2vb7WiPEwARsYZkRu9HlFyP8E/kXHrk58DRks4vUKdG0uBW3mMVyeSOi9Jzu5rXG0cyS/rRiHilM/G34i37Ku0x/WwR2/wNyezYyyXlfgYXkwy156olmUX9xmVQ0n3ZPLlmfhFxtKmFcyafTp8LnRLQrBY4Mp2F3bytvsCnSXqt782rPzftXW2uO5qkZ/HRnH8Xnf7umJUL9/iZVaiIeFXSHSQ9SpA3zJt6iGRG5U8lXUMyxDiXpBewGItJeui+peRafBtIZvQew5vno3XGZSQxPyHpOuA5YBzJ5VzOBVYBXyeZkHCbpJuBP5AMXU8lOc/xXbw5eaWQz5JcLuX3kv6LJFG9lGSyxN8WEXtL7gT+DbgzbVMfkkkGfTu7wXT49N+BK4G7JN1O2gMLvJhX/TqShPC6dMh1Kcm+PBf4dkQUOjeuK/04PeXgXpJh7lEkl895jSS5a8m/knyetel3t4Fkvx0HfKbAf3KeAR6S9D1gN8nlXGqA/5dTp9jvjlnJc+JnVtluIkn8XiFJMP5ERGxIhxS/AfwTSeJ3B0nvT2Nn3zQidkmaQ3Kdu+Zk6X6SWZEPFbHdeiUXIf4S8FfAPiQ/+HeTJpQRsV3SbJIf9A+Q9Oq8RnLtuu+QJKWtvcd96YV9v0RyrcIAfk9yCZrfdzb2Vt7vWUnvJrlczr8Bm0mulfgjkt7HzrqK5LzKS0lmQz8NnENyrb3c99+R7q9/JkluhpEkh58D/rOI92+vG0kuZH1J+t4bSGbWfiUi8ieMvCEi1kk6gWSfXUJyTcmlwF9GxI0FVvk5ySSjvyW5PuIzwJzcUxyK/e6YlYPm616ZmZlVHEmTSc5zvDIi/i3reMyy5nP8zMzMzKqEEz8zMzOzKuHEz8zMzKxK+Bw/MzMzsyrhHj8zMzOzKuHEz8zMzKxKOPEzMzMzqxJO/MzMzMyqhBM/MzMzsyrhxM/MzMysSvx/0DW6enxwqBEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(likev[:,0], 2*likev[:,1])\n",
    "plt.xlabel(\"Variance of random slope\", size=17)\n",
    "plt.ylabel(\"-2 times profile log likelihood\", size=17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a plot of the profile likelihood function. The profile likelihood plot shows that the MLE of the random slope variance parameter is a very small positive number, and that there is low uncertainty in this estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '-2 times profile log likelihood')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAHoCAYAAADAJXJqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3yV5f3G8c/3ZLIDIayEvfeeIm5QEVTQirtubW0dte6Bs1r156h1tVqtVkUBq4iIxQEqS6bsEXZYSdiEkHX//sjBIiUhh5yT55yc6/16nVc9z3MSLirIxfM89/c25xwiIiIiErl8XgcQERERkfJRoRMRERGJcCp0IiIiIhFOhU5EREQkwqnQiYiIiEQ4FToRERGRCBfrdQCv1a1b1zVr1szrGCIiIiLHNHfu3CznXMqRx6O+0DVr1ow5c+Z4HUNERETkmMxs/dGO65ariIiISIRToRMRERGJcCp0IiIiIhFOhU5EREQkwqnQiYiIiEQ4FToRERGRCKdCJyIiIhLhVOhEREREIpwKnYiIiEiEU6ETERERiXAqdCIiIiIRToVOREREJMKp0ImIiIhEOBU6ERERkQinQiciIiIS4VToRERERCKcCp2IhBXnHAWFRV7HEBGJKLFeBzgaM7sQGA20B/o45+b4j58BPAnEA3nAH51zX/vPfQs0BA74v81g59z2ik0uIsdjT24+01dnM21VJtNWZrJtTy6ntavPyJ5pnNw2hbgY/d1TRKQ0YVnogMXACOC1I45nAcOcc5vNrBMwGUg97Pylh8qfiISvwiLH4ozdTFuZybRVmczbsIvCIke1+Bj6t6zLyW1TmLRoK18s2UqdavEM79qIkT3S6JRaEzPzOr6ISNgJy0LnnFsG/M9/uJ1z8w97uwRINLME59zBCownIsdh255cpq7M5LtVWXy/KpOdOfkAdEqtyQ2DWjCoTQo9mtQmPrb4atxDwzoybWUm4+dl8N6sDbw1fR1t6ldnZI80zuueSv2aiV7+dEREwoo557zOUCL/bdQ7jnbVzcwuAG50zp1+2GeTgUJgHPCYK+EnZ2bXA9cDNGnSpOf69etDkl8kmuXmF/Ljuh3FV+FWZrFi214A6lZPYFCbugxqncLA1nWpWz3hmN9rd04+E37azLh5m5i/YRc+g4GtUxjZI5XBHRpQJT4m1D8dEZGwYGZznXO9/ue4V4XOzKYADY5y6j7n3Cf+z3zLUQqdmXUEPqX4Obl0/7FU51yGmdWguNC965z757Fy9OrVy82Zo7u0IuXlnCM9cx9TV2YxbWUms9Zmk5tfRHyMj17NajOoTQqDWqfQvmGNct02XZO5j/HzMvh4fgYZuw5QPSGWoZ0bMqJHKn2a19EtWRGp1MKu0JXF0QqdmaUBXwNXOed+KOHrfg30cs7dfKwfQ4VO5Pjtzsnn+9XFBe67VZls3p0LQIu61YoLXJu69GuRTNX44D/dUVTkmLk2m/HzMvh80RZy8gppXKcKI7qnMaJHKk2TqwX9xxQR8VqlKHRmlgRMBR5xzo077HOxQJJzLsvM4oD3gSnOuVeP9WOo0ImUXUFhEQs3/Xcxw8KNuyhyUCMhlhNa1WVQmxRObF2XxnWqVmiunLwCvli8lfHzMvghPQvnoHez2ozokcbQLg2pmRhXoXlEREIlogqdmZ0P/AVIAXYBC5xzQ8zsfuAeYNVhHx8M7AemAXFADDAFuN05V3isH0uFTqR0m3cd+LnAfb8qiz25BZhBl7QkTmpdXOK6NU4iNkxGi2zedYCP52cwbt4m1mTuJyHWx+CODRjRI5UTW9UNm5wiIscjogpdRVKhE/mlA3mFzFyb7V/MkEl65n4A6tdMYFDrFAa1SWFgq7rUrhbvcdLSOedYuGk34+dt4tOFm9mVk09KjQTO757KiB6ptGtQ0+uIIiIBU6ErgQqdRDvnHCu27f15NersdTvIKygiIdZHn+Z1OKlNcYlrXa96xC44OFhQyDfLtzNuXgbfLN9OQZGjY6OajOyRxvBujcq00lZEJByo0JVAhU6i0e4D+Uz1X4H7blUm2/YUj3JsXa+6fzFDCn2b1yExrvKNA8ned5AJCzczbl4GizJ2E+szTm6bwogeaZzWvh4JsZXv5ywilYcKXQlU6CTaLNuyhyvfnM32vQepVSWOga3qMqhNXU5snUKjpCpex6tQK7ftZdy8Tfx7fgbb9hT//zGsa0NG9Eije+OkiL0iKSKVlwpdCVToJJrMWpPNtf+cQ7X4WJ4f1Y3ezeoQ41NpKSxyfL86i/HzNjF5yVZy84toUbcaI3qkcn6PNFKjrOiKSPhSoSuBCp1Eiy+XbOXm9+eTVrsK71zTVyWlBHtz85m0aCtj521i9todmEH/FsmM6JHGWZ0aUC0hLHdMFJEooUJXAhU6iQYf/riRu8f/ROe0JP7x697UCfMVquFiQ3YOH8/PYPz8TazPzqFqfAxndmrAyB5p9G+RjE9XN0WkgqnQlUCFTioz5xyvTl3DU18s58TWdXn1sp66wnQcnHPMXb+TcfM28dnCLew9WECjWomc3yOVET3SaJlS3euIIhIlVOhKoEInlVVRkePxz5fxxvdrGd61Ec9c2JX4WA3VLa/c/EL+s3Qb4+ZtYtrKTIocdGucxMgeqQzr2oikqrr6KSKho0JXAhU6qYzyC4u4c+xPfDw/g18PaMaD53TQ7cEQ2L4nl08WbGbcvE0s37qX+Bgfp7arx8ieaZzcNoU47UohIkGmQlcCFTqpbHLyCvjNv+bx7YpM7hjcht+e0krjN0LMOcfSLXsYNzeDTxZkkL0/j+Rq8Qzv1oiRPdLo2Kim/h2ISFCo0JVAhU4qk105eVz11o8s3LiLx8/vzMV9mngdKerkFxYxbWUm4+ZtYsrS7eQVFtG2fo3iESjdU6lXM9HriCISwVToSqBCJ5XFlt0HuOKN2azfkcOLo7pxZqeGXkeKertz8pnwU/Et2fkbduEzOLF1CiN6pDKkY4NKuROHiISWCl0JVOikMli9fR9XvDGLPbkF/O2KXvRvmex1JDnCmsx9jJ+XwcfzM8jYdYAaCbGc3bkhI3um0btZbd2SFZEyUaErgQqdRLr5G3Zy9Vs/EuPz8dZVvemUWsvrSFKKoiLHzLXZjJubwaTFW8jJK6RxnSqM6J7GBT3TaFynqtcRRSSMqdCVQIVOItnUlZnc+M5cUmok8M41fWiaXM3rSBKAnLwCvli8lfHzMvghPYs4n4+XL+3B6R3qex1NRMJUSYVOa+pFItQnCzK45q0faVa3GmNv6q8yF4Gqxscyokca717bl+/vOpV2DWtw47tzmbRoi9fRRCTCqNCJRKB//LCWWz5YQI+mtRlzQz/q1dDKyUiXmlSFd6/tS5e0Wtz8/nw+WZDhdSQRiSAqdCIRxDnHM5NX8PCEpQzuUJ9/Xt2HmolxXseSIKmZGMc/r+lLr6a1uW3MAsbO3eR1JBGJECp0IhGisMhx78eLeOmb1Yzq3ZiXL+2hsReVUPWEWN66qg8DWtblj2MX8t6sDV5HEpEIoEInEgFy8wv5zb/m8v7sjdx8Siv+NKIzsdpWqtKqEh/D36/sxcltUrj340W8PX2d15FEJMzpTwSRMLcnN58r35zN5CXbeGhYB+4Y0lYzy6JAYlwMr17ek8Ed6vPQp0v427Q1XkcSkTCmQicSxrbvzWXUazOZu34nL4zqxlUnNPc6klSghNgY/nppD4Z2acjjny/jpa9XeR1JRMJUrNcBROTo1mfv5/I3ZpO592Dx7be29byOJB6Ii/HxwkXdSIjx8cyXKzlYUMTtZ7TRVVoR+QUVOpEwtGTzbq5880cKiop477q+dG9S2+tI4qHYGB9PX9iVuBgff/l6NXkFRdx9VjuVOhH5mQqdSJiZuSab696eQ43EWD64vj+t6tXwOpKEgRif8acRnYmP9fHatDUcLCjioWEdVOpEBFChEwkrXyzeyu8/mE+TOlX559V9aJRUxetIEkZ8PuORczsSH+vjje/XkldYxGPndsLnU6kTiXYqdCJh4oPZG7j340V0bZzEm1f2pna1eK8jSRgyM+4f2p74WB+vfJtOXkERT43sQoxKnUhUU6ET8Zhzjpe/TefpySs4qU0Kr1zWg6rx+q0pJTMz7hzSloRYH89PWUV+YRHPXthVswlFopj+1BDxUFGR45HPlvLW9HWc163Rzw++ixyLmXHr6W2Ij/Xx5y9WkFdQxAujuhMfq18/ItFIhU7EI3kFRdzx0UI+XbiZawY2576z2+tZKAnYb05uRXyMj8cmLiP/X3P566U9SIjVlnAi0UZ/lRPxwP6DBVz7zzl8unAzd53ZjvuHqszJ8bv2xBY8em5HpizbznX/nEtufqHXkUSkgqnQiVSwHfvzuOTvs/h+VSZPjezMTSe31OgJKbfL+zfjqZGd+W5VJle/9SM5eQVeRxKRCqRCJ1KBMnYd4MJXp7Nsyx5evawnF/Vu4nUkqUQu6t2E//tVV2auyebKN2ezNzff60giUkFU6EQqyKpte7nglels33uQd67uw+CODbyOJJXQ+d3TePHi7szbsIvL35jN7gMqdSLRQIVOpALMXb+TC1+bQUGRY8z1/enbItnrSFKJndOlES9f2oMlm3dz6d9nsnN/nteRRCTEVOhEQuybFdu57O+zSKoSx7gbB9ChUU2vI0kUGNKxAa9f3ouV2/Zx8d9mkrXvoNeRRCSEVOhEQujf8zO47u05tEipxkc3DqBJclWvI0kUOaVdPd64shfrsvcz6vWZbN+T63UkEQkRFTqREHnj+7XcOmYBvZvV4YPr+5FSI8HrSBKFTmydwltX9WHzrgNc9PpMtuw+4HUkEQkBFTqREHh+ykoe/WwpZ3VqwD+u6k2NxDivI0kU69cimXeu6UPW3oP86rUZbNyR43UkEQkyFTqRIFu4cRcvfLWKEd1TeemSHiTGaWq/eK9n0zq8e21fdufkc9FrM1iXtd/rSCISRCp0IkFUVOR46NMl1K2ewMPndiRGuz9IGOnaOIn3r+/HgfxCfvXaDFZv3+d1JBEJEhU6kSAaN28TCzbu4u4z2+k2q4Sljo1q8cH1/SlyMOr1GazYutfrSCISBCp0IkGyJzefp75YTvcmSZzfPdXrOCIlatugBmNu6EeMzxj1+gwWZ+z2OpKIlJMKnUiQvDBlFdn783hkeCd8utUqYa5lSnU+vKE/VeNjueRvM1mwcZfXkUSkHFToRIJg1ba9vD19HaN6N6ZzWi2v44iUSdPkaoy5oR9JVeO57O+zmLNuh9eRROQ4qdCJlJNzjtETllA1PoY7Brf1Oo5IQNJqV2XMDf2oVyOBK96czYz0bK8jichxUKETKafJS7byw+ps/jC4LcnVNTxYIk/DWlX44Pp+pCZV4df/mM20lZleRxKRAKnQiZTDgbxCHv1sGe0a1ODSvk28jiNy3OrVTOSD6/vRIqU61749h6+Xb/M6kogEQIVOpBxenZpOxq4DjB7ekdgY/XaSyJZcPYH3r+tL2wY1uOGduXyxeKvXkUSkjPQnkMhx2rgjh1enpnNOl4b0a5HsdRyRoEiqGs+/rutLp9Ra/Pa9eUxYuNnrSCJSBip0IsfpsYlL8Zlx39D2XkcRCaqaiXG8c01fejapzS0fzGfc3E1eRxKRY1ChEzkO363KZPKSbdx8aisa1qridRyRoKueEMtbV/emf8tk7hi7kA9mb/A6koiUQoVOJED5hUU8PGEpTZOrcu2Jzb2OIxIyVeNjeePK3gxqncLd4xfxzxnrvI4kIiVQoRMJ0NvT17F6+z4ePKcDCbExXscRCanEuBhev6Inp7evz4OfLOHv363xOpKIHIUKnUgAtu/N5fkpqzi5bQqntqvndRyRCpEQG8Mrl/VgaOeGPDZxGX/9ZrXXkUTkCLFeBxCJJE9NWsHBgkIePKcDZtqvVaJHXIyPF0Z1Iy7GeHryCg4WFHHb6a31+0AkTKjQiZTR3PU7GTdvEzee1JIWKdW9jiNS4WJjfDz7q27Exfh48atV5BUUcdeZbVXqRMKACp1IGRQVOUZ/uoT6NRP43amtvI4j4pkYn/HUyC7Ex/p4dWq6rliLhAkVOpEy+HDORhZl7OaFUd2olqDfNhLdfD7jsfM6ER/r4x8/rCOvoIhHz+2Ez6dSJ+IV/ckkcgy7c/L58+QV9G5Wm+FdG3kdRyQsmNnPK71fnZpOfmERfxrRhRiVOhFPqNCJHMNzU1ayKyeP0cP76LaSyGHMjLvObEt87H+fqXvmwq7a11jEAyp0IqVYvnUP78xczyV9m9CxUS2v44iEHTPj9jPakBDr4+nJK8grLOKFUd2JU6kTqVAqdCIlcM7x0CdLqJEYyx/OaOt1HJGw9ttTWpEQ6+OxicvIK5jHXy/trsHbIhVIf4USKcHERVuYtXYHdwxuS+1q8V7HEQl7157YgkfO7ciUZdu44Z255OYXeh1JJGqo0IkcRU5eAY9PXEbHRjW5uE8Tr+OIRIwr+jfjTyM6M3VlJte8/SM5eQVeRxKJCip0Ikfx8jfpbNmdy8PDO2rVnkiALu7ThGcu6MqM9Gx+/eaP7DuoUicSaip0IkdYn72f16et4bxujejVrI7XcUQi0sieabwwqjtzN+zk8jdmsftAvteRRCo1FTqRIzz62VLiYox7zm7vdRSRiDasayP+ekkPFmfs5rK/z2JXTp7XkUQqLRU6kcN8s2I7U5Zt53entaZ+zUSv44hEvDM7NeC1y3uyYtteRr0+k+x9B72OJFIplVjozKzIzAoDfQUjlJldaGZL/Bl6HXb8DDOba2aL/P976mHn4s3sdTNbaWbLzWxkMLJI9MgrKOLRCUtpUbcaV5/Q3Os4IpXGqe3q8/crerE2az+jXp/J9j25XkcSqXRKm0P3COCOOHYe0AmYDKwADGgLDAYWAZ8EKddiYATw2hHHs4BhzrnNZnYoR6r/3H3AdudcGzPzAXr4SQLy5g9rWZO1n7eu6k18rC5eiwTToDYpvHVVH655+0dGvT6Tf13Xl4a1qngdS6TSKLHQOedGH/7ezK4GGgCdnHMrjjjXHvgG2BCMUM65Zf7ve+Tx+Ye9XQIkmlmCc+4gcDXQzv+5IorLn0iZbNuTy1++WsXp7etzctt6XscRqZT6t0zmn1f34df/+JFfvTaD967tR+M6Vb2OJVIpBHIZ4k7gpSPLHPxcwP4K3BWsYGUwEpjvnDtoZkn+Y4+a2Twz+8jM6ldgFolwf/p8GflFjgfO0UIIkVDq1awO717bl905+Yx6fSbrs/d7HUmkUgik0DUFDpRyPsf/mTIxsylmtvgor3PL8LUdgaeAG/yHYoE04AfnXA9gBvBMKV9/vZnNMbM5mZmZZY0sldSP63bw7wWbuf7EFjRNruZ1HJFKr1vjJN67rh85eQVc9JoWSogEQyCFbiVw/WFXw35mZrWB6yl+rq5MnHOnO+c6HeVV6nN4ZpYGfAxc4ZxL9x/OprhQfux//xHQo5Qf+3XnXC/nXK+UlJSyRpZKqLCoeL/WRrUS+c0pLb2OIxI1OqXW4p1r+rJjfx53jVuEc0c+si0igQik0N0LtARWmdn/mdmNZnaDmT1HcdlrQfHChJDxl8mJwD3OuR8OHXfF/yWYAJzsP3QasDSUWaRyeH/2BpZu2cO9Q9tTNb60NUIiEmydUmtx11ntmLJsG+/NDsoj2CJRq8yFzjk3ERhC8cKHW4GXgVeAW/zHzvJ/ptzM7Hwz2wT0Byaa2WT/qZuBVsADZrbA/zr0BPtdwGgz+wm4HPhDMLJI5bVzfx7PfLmC/i2SGdq5oddxRKLSVQOacWLrujz62VJWb9/ndRyRiGXHc5nbv+CgGcVjS9Y657YFOVeF6dWrl5szZ47XMcQD9/97Ee/P3sjE3w+kXYOaXscRiVrb9+Qy5PlpNEqqwse/OUFjg0RKYWZznXO9jjx+XL9rnHPbnHOznHMzI7nMSfRasnk3783awOX9mqrMiXisXs1EnhrZhSWb9/Dsf8r8KLaIHCagQmdmSWb2pH816n4z2+f/5yeOtlhCJBw55xj96RKSqsZz2+ltvI4jIsDgjg24pG8TXp+2humrNUZUJFBlLnRmlgrMp3geHRQvTphE8W4SdwPzzKxR0BOKBNmnCzfz47qd3DmkLbWqxnkdR0T87h/anuZ1q3H7hwvZlZPndRyRiBLIFbo/AfWBc/zjRX7lnLvQOdcZGOo/96dQhBQJlv0HC3ji82V0SavFr3o19jqOiBymanwsL47qTvb+g9wzXqNMRAIRSKE7E3jBOff5kSecc5OAvwBnBSuYSCj85evVbNtzkNHDO+Lz2bG/QEQqVKfUWvxhcFsmLd7KR3M3eR1HJGIEUuhqABmlnN/k/4xIWFqTuY83vl/DyB5p9GhS2+s4IlKC609sQf8WyYz+dAnrsrQ1mEhZBFLoVgAXmNn/fI2ZxQAXEMBOESIVyTnHI58tJSE2hrvOaut1HBEphc9nPPurrsTF+LhlzALyC4u8jiQS9gIpdC8Cg4CvzexcM2vnf50HTAFOBJ4PRUiR8vpq2Xa+XZHJrae3pl6NRK/jiMgxNEqqwhPnd2bhxl28+NUqr+OIhL0y73XknHvTvyvDQ8D4w04ZcBC41zn3VnDjiZRfbn4hj3y2lFb1qnPlgGZexxGRMhrapSHfrkjjr9+s5sTWKfRpXsfrSCJhK6A5dM65J4E04FKK93a9F7gYSHXOPRX8eCLl98b3a9mwI4fRwzoSF6MJ9CKR5KHhHWlcpyq3jVnA7gP5XscRCVsB/+nmnMt2zn3gnHvK/xrjnNsRinAi5bV51wFe+no1Z3ZswMDWdb2OIyIBqp4Qy/MXdWPrnlwe/GSx13FEwlbAhc7MzjSzl8xsopl95v/nwaEIJ1JeT3y+jCLnuG9oe6+jiMhx6t6kNree1ppPFmzm3/NLG7YgEr0C2Ski3sw+oXiHiN8AvYG+/n+eZGb/NrP40MQUCdyM9Gw++2kLN57UksZ1qnodR0TK4TentKJ3s9o88O/FbNyR43UckbATyBW6h4BhwLNAPedcPedcCpACPAMMBx4IfkSRwBUUFvHwhCWkJlXhppNbeh1HRMopxmf836+6AXDbmAUUaJSJyC8EUuguAd51zt3pnPt552T/M3V3Ae8ClwU7oMjx+NesDSzfupcHzmlPYlyM13FEJAga16nKY+d3Ys76nbzybbrXcUTCSiCFrhEwvZTzM4CG5YsjUn7Z+w7y7JcrGNiqLkM6NvA6jogE0bndUjmvWyOe/2oV8zfs9DqOSNgIpNBtBvqVcr4vsKV8cUTK75kvV5CTV8jo4R0w036tIpXNI+d1okHNRG4ds4B9Bwu8jiMSFgIpdO8Dl5vZY2b280aYZlbbzB4FLgfeC3ZAkUD8tGkXH/y4kSsHNKNVPW0tLFIZ1UyM4/lR3di4I4eHP13idRyRsBBIoXsYmEzxMOEsM9tqZluBLOA+/7lHgh9RpGyKihwPfbqE5GoJ3HJ6a6/jiEgI9W5Wh5tPacVHczcx8SfdHBIJZOuvg8DZZnYOMBRo5j+1DpjgnPs86OlEAjB+fgbzN+zi6Qu6UDMxzus4IhJivzutNdNWZXHP+J/o3iSJRklVvI4k4pnj2SniM+fcTc65s/yvm1TmxGt7c/N5ctJyujdJYmSPNK/jiEgFiIvx8cKobhQWOW7/cAGFRc7rSCKe0caWUim8+NUqsvcf5OHhHfH5tBBCJFo0Ta7G6OEdmblmB3/7bo3XcUQ8U+ZbrgBmdgZwLdACqAMc+Senc85piqtUqNXb9/KPH9ZxUa/GdElL8jqOiFSwC3qm8e2KTJ79cgUntKxL57RaXkcSqXCBbP11G/AFcBKQAUwDph7xmhaCjCIlcs7x8ISlVImP4Y4hbb2OIyIeMDMeP78TdasncMuY+eTkaZSJRJ9ArtDdRnFpO9M5lxeiPCIBmbxkG9+tyuKhYR2oWz3B6zgi4pGkqvE8+6uuXPr3WTw2cRlPnN/Z60giFSqQZ+jqAmNU5iRc5OYX8tjEpbStX4PL+zX1Oo6IeGxAy7rcMKgl783awJdLtnodR6RCBVLo5lH87JxIWHht6ho27TzA6OEdiY3R+h4RgdvPaEOn1JrcNe4ntu/J9TqOSIUJ5E/BWyneKeKMUIURKatNO3N4+dvVDO3SkP4tk72OIyJhIj7WxwujunMgv5A/fLSQIo0ykShR4jN0ZvblUQ7vBb4ws3UUDxQuPOK8c84NCVo6kRI8PnEZZnDf2e29jiIiYaZlSnUeOKcD9328mH9MX8c1A5t7HUkk5EpbFNEGONpfbTZQfGVPt1/FEz+szmLS4q384Yw2mgwvIkd1SZ8mfLM8k6cmLWdAy2TaN6zpdSSRkCqx0DnnmlVgDpEye+4/K0mrXYXrBunvFCJydGbGUyM7c+YL33HLB/P59OaBJMbFeB1LJGT0JLlElLVZ+5mzfieX9Wuq/ziLSKmSqyfwzIVdWbltH09OWu51HJGQUqGTiDJu7iZ8Bud3T/U6iohEgJPapHD1Cc15a/o6vlm+3es4IiFTYqEzsyIzKzCz+MPeFx7jpfHcEjKFRY5x8zYxqE0K9Wsmeh1HRCLEnWe2pV2DGvxx7EKy9h30Oo5ISJS2KOIRihdFFBzxXsQTM9Kz2bI7l/uGamWriJRdYlwML4zqzrCXvufOsT/xxpW9MDtyK3KRyFbaoojRpb0XqWhj526kZmIsp7ev73UUEYkwbRvU4N6z2jF6wlLenbmey/s38zqSSFDpGTqJCHty8/liyVaGd2ukxRAiclyuHNCMk9qk8NjEZazattfrOCJBVdpg4UHH8w2dc9OOP47I0X3+0xZy84u4oGdjr6OISIQyM56+sAtnPf8dv/9gAf/+7QASYvUXRKkcSnuG7lsCe2bO/J/X7w4JurFzN9GqXnW6ptXyOoqIRLB6NRL58wVduObtOTwzeQX3De3gdSSRoCit0J1SYSlESnFo9tzdZ7XTg8wiUm6nta/P5f2a8rfv1nJSm3oMbF3X60gi5VbaooipFRlEpCSaPSciwXbv2e2ZsSab2z9cwORbB1G7WrzXkUTK5bgWRZhZazM7wcx0/0tCSrPnRCQUqsTH8MKobuzMyePu8bo+4MIAACAASURBVD/hnKZySWQLqNCZ2UVmth5YDkwDevqP1zWzVWZ2YQgyShQ7NHtuZI80r6OISCXTsVEt7hzSjslLtjHmx41exxEplzIXOjM7F3gf2AA8QPEiCACcc1nAMuDyYAeU6DZ27kZqJMZyRgfNnhOR4LtmYHNOaJXMwxOWsiZzn9dxRI5bIFfo7gemOedOBF47yvlZQNegpBLhsNlzXTV7TkRCw+cznr2wGwlxPm75YAF5BUVeRxI5LoEUuo7Ah6Wc3wroMooEzX9nz+l2q4iEToNaiTw5oguLMnbz/JSVXscROS6BFLpcoLSn0psBu8qVRuQwY+duomVKNbo1TvI6iohUcmd2asCo3o15ZWo6M9dkex1HJGCBFLrvgYuPdsK/2vVq4OtghBI5NHvugp6NNXtORCrEA+d0oFlyNW4fs4DdOflexxEJSCCFbjTQ0cy+AUb4j/Uys5uBBUBN4NHgxpNopdlzIlLRqiXE8vxF3di+9yD3/nuRRplIRClzoXPOzQOGAA3476KIJ4EXgTxgiHNuWdATStQ5NHvuxNYpNKil2XMiUnG6Nk7itjPaMPGnLYyfl+F1HJEyK23rr//hnPsOaG9mXYE2FBfC1cA8p7/KSJAcmj1379ntvY4iIlHoxpNaMnVlJg9+sphezWrTNLma15FEjimQOXQ/LzV0zi10zn3knBvjnJt7qMyZ2UmhCCnRRbPnRMRLMT7juYu64fMZt41ZQEGhRplI+AvkGbopZlbiDsZmNhiYWP5IEs00e05EwkFqUhWeOL8z8zbs4qVvVnsdR+SYAil0McB/jrZ/q5kNAz4FZgYrmEQnzZ4TkXAxrGsjRvRI5cWvVjF3/Q6v44iUKpBCdzqQDEwys58fKDCzC4BxwDfAOcGNJ9FGs+dEJJw8PLwjqbWrcOuYBezN1SgTCV+BrHJdD5wBtAA+MbMEM7uc4v1dPwfOdc7lhiamRAPNnhORcFMjMY7nL+rO5l25PPTpEq/jiJQokCt0OOdWUDy6pAfFt1f/QfHVuZHOubzgx5NootlzIhKOejatze9ObcX4eRl8unCz13FEjiqgQgfFK1yBs4CWwLvAxc65wmAHk+ii2XMiEs5uPqUVPZokcd/Hi8jYdcDrOCL/o8RCZ2b5ZpZ3tBfwHVANuAQ4eNi5gxUVXCqXQ7PntBhCRMJRbIyP5y/qjnNw25gFFBZp9KqEl9IGC/8L0K9YqRCaPSci4a5JclUeObcjt3+4kFenpvPbU1p5HUnkZyUWOufcryswh0SxQ7PnRvZI0+w5EQlr53dP5ZsVmTz3n5UMbFWXrlqRL2Ei4GfoRIJNs+dEJFKYGY+d14l6NRK4dcwC9h8s8DqSCFDKFTozGwTgnJt2+PtjOfR5kbLS7DkRiSS1qsTxfxd14+K/zeTRz5by5MguXkcSKfUZum8BZ2ZV/CNJvqX0Z+rMf173zKTMDs2eu+vMdpo9JyIRo1+LZG46qSUvf5vOyW1TOLNTQ68jSZQrrdCdAnDYfLlTQh9Hoo1mz4lIpLr19DZ8vzqLu8cvolvj2hq5JJ4qbVHE1NLei5SXZs+JSCSLj/Xx/EXdGPri9/zhowW8c3VffD7daRBvaFGEeEaz50Qk0rVIqc5Dwzrww+ps3vh+rddxJIqVtijiweP4fs4592g58kgU0ew5EakMLurdmG9WbOfPk5czoFUyHRvV8jqSRKHSnqEbfRzfzwEqdHJMmj0nIpWFmfHkiC6c+cI0bvlgARNuHkiVeP13TSpWibdcnXO+43jpV7CUiWbPiUhlUrtaPM9e2I3V2/fxxOfLvI4jUSgsn6EzswvNbImZFZlZr8OOn2Fmc81skf9/T/Ufr2FmCw57ZZnZ8979DORYNHtORCqbga3rct2JzXln5nq+WrbN6zgSZcKy0AGLgRHAkUOKs4BhzrnOwJXAOwDOub3OuW6HXsB6YHxFBpayOzR77oKejTV7TkQqlTuGtKV9w5rcOfYntu/N9TqORJGwLHTOuWXOuRVHOT7fObfZ/3YJkGhmCYd/xsxaA/WA70KfVI6HZs+JSGWVEBvDi6O6se9gAX/86CecK20ev0jwhGWhK6ORwHzn3MEjjl8MjHGl/C4ys+vNbI6ZzcnMzAxpSPklzZ4Tkcqudf0a3D+0PVNXZvL29HVex5Eo4VmhM7MpZrb4KK9zy/C1HYGngBuOcnoU8H5pX++ce90518s51yslJeX4fgJyXDR7TkSiwWX9mnJqu3o8MWk5K7bu9TqORAHPCp1z7nTnXKejvD4p7evMLA34GLjCOZd+xLmuQKxzbm4Io0s5aPaciEQDM+PPF3ShZmIst3wwn9z8Qq8jSSUXUbdczSwJmAjc45z74SgfuZhjXJ0T7+z1z54b3rWRZs+JSKVXt3oCT1/YleVb9/LnL/7nsXCRoCptsPAvlGHnCAfkApuAqYctXgiYmZ0P/AVIASaa2QLn3BDgZqAV8ICZPeD/+GDn3Hb/P/8KOPt4f1wJrc8XafaciESXU9rW49cDmvHmD2s5qW0KJ7XRYz4SGlbWFThmVkRxaQM4ctbEkccLgVeB35e2OCEc9OrVy82ZM8frGFHhwlens2N/HlNuP0njSkQkauTmFzL8pe/ZmZPPF7ecSHL1hGN/kUgJzGyuc67XkccDueWaBiwE3gV6AbX8rz7Av4AFQGugJ/AB8BvgrvLFlspiXdZ+flyn2XMiEn0S42J4YVR3dufkc9e4RRplIiERSKH7C5DunLvSOTfPP8x3r3NujnPuCmAt8IR/VtzlwFfAr0OQWSLQuHmaPSci0at9w5rcdVY7pizbxicLjvuJJJESBVLoTqe4pJXkK2DIYe8/A5odRyapZIqKHOPmavaciES3qwY0o0taLf40aRn7DxZ4HUcqmUAKXRHQtZTzXfnvs3SH5AScSCqdGWuy2azZcyIS5Xw+46FhHdm25yCvfJt+7C8QCUAghe7fwHVmdreZ1Th00MxqmNk9wLX+zxzSH1gZnJgSycbO3aTZcyIiQM+mtTm/eyqvf7eGjTt0zUOCJ5BCdxswA3gC2GlmW8xsC7ATeByY5f8MZpYI7AL+L7hxJdLszc1n0uItmj0nIuJ315ntiPUZj09c5nUUqUTKPIfOObfLzAYB5wFnAU38p9YDk4BPDo0occ7lAjcFOatEIM2eExH5pQa1EvntKa14evIKpq/OYkCrul5HkkqgzIUOwF/YPva/RI5p7NxNtEypRrfGSV5HEREJG9cMbM77szfw8ISlTPz9QGJjImrjJglDAf8KMrPaZnaBmf3RzO4ws5H+LblEfkGz50REji4xLob7h7Znxba9vD97g9dxpBII6Aqdmd0OPAok8svdIg6Y2f3OueeCGU4im2bPiYiUbEjHBgxomcyz/1nJsK6NSKoa73UkiWBlvkJnZlcCz1C8I8RFQCegM8X7p84HnjGzK0IRUiLPodlzAzV7TkTkqMyMB4d1YM+BfJ77j4ZCSPkEusr1e2CQc26sc26pc26Jc24scBLwA3B7KEJK5NHsORGRY2vXoCaX9WvKu7M2sGLrXq/jSAQLpNC1BT50zhUeecJ/7EP/Z0R+nj03WLPnRERKddvpbaieEMsjny3RPq9y3AIpdHuBRqWcT/V/RqLcodlzwzR7TkTkmGpXi+cPg9vww+psvly6zes4EqECKXRfAr83s9OPPGFmpwE3A5ODFUwil2bPiYgE5pI+TWhbvwaPTVxKbv7/3AgTOaZACt3dFO8KMdnMFprZB/7XAorL3k7gnlCElMgydu4mWqRUo7tmz4mIlElsjI8Hh3Vg444DvPH9Wq/jSAQqc6Fzzm0CugHPAfHAuf5XAvAs0N3/GYli/509l6bZcyIiATihVV2GdKzPX79ZzdbduV7HkQgT0GBh51y2c+4O51x751wV/6u9c+5O51x2qEJK5Dg0e25Ed91uFREJ1H1nd6CgyPHnL5Z7HUUijPYakaDR7DkRkfJpklyV605szvj5GczbsNPrOBJBStwpwswePI7v55xzj5Yjj0SwQ7Pn7j67vddRREQi1m9ObsXYuZt4+NMlfPybE/D59PiKHFtpW3+NPo7v5yjeGkyikGbPiYiUX7WEWO4+qx23jVnIuHmbuLBXY68jSQQo8Zarc853HC8NHYtSmj0nIhI853ZNpXuTJJ76YgV7c/O9jiMRQM/QSVBo9pyISPD4fMboYR3J2neQl75Z7XUciQAqdBIUmj0nIhJcXRsncUHPNN78fi1rs/Z7HUfCnAqdlJtmz4mIhMadZ7YlPsbH4xOXeR1FwpwKnZSbZs+JiIRGvRqJ/O601kxZto1pKzO9jiNhTIVOykWz50REQuuqE5rRLLkqj3y2lPzCIq/jSJhSoZNyOTR7ToshRERCIyE2hvuHdmD19n28M2O913EkTB1XoTOzRDNLNbP4YAeSyKLZcyIioXda+3oMapPCc1NWkr3voNdxJAwFVOjMbKCZfQfsBTYAA/3H65rZV2Y2OAQZJUxp9pyISMUwMx48pz0H8gp59j8rvY4jYajMhc7MBgJfAQ2AvwM/L2d0zmX5318T7IASvjR7TkSk4rSqV4Mr+jfj/dkbWLJ5t9dxJMwEcoXuMWAp0Al44CjnpwK9gxFKIoNmz4mIVKxbTm9N7arxPDxhKc45r+NIGAmk0PUC3nLOHaR4z9YjbaL46p1EAc2eExGpeLWqxHHH4LbMXruDzxdt9TqOhJFACl0RRy9yhzQCcsoXRyKFZs+JiHjjot6Nad+wJk98vowDeYVex5EwEUih+xEYfrQT/tWulwHTgxFKwptmz4mIeCfGZ4we1oGMXQd4fdoar+NImAik0D0BnGxm/6T49itAYzM7B5gGNPd/Rio5zZ4TEfFW3xbJDO3SkFemriZj1wGv40gYKHOhc859BVwCnA187j/8JvAp0Aa4xDk3M+gJJexo9pyIiPfuOasdzsGTk5Z7HUXCQEBz6JxzHwJNgBHAXcC9wIVAE+fc2ODHk3Cj2XMiIuEhrXZVbjypJRMWbmb22h1exxGPBbxThHMuxzn3iXPuaefcU865cc65faEIJ+FHs+dERMLHjSe1pFGtRB6esITCIo0xiWbay1UCotlzIiLho0p8DPec3Z4lm/fw0ZyNXscRD5VY6MysyMwKA3wVVGR4qViaPSciEn7O6dKQPs3q8PTkFew+kO91HPFIbCnnHqH0uXMSZTR7TkQk/JgZDw7rwLCXvucvX63i/nM6eB1JPFBioXPOja7AHBLmNHtORCR8dUqtxajejXlr+jpG9WlCq3rVvY4kFUzP0EmZaPaciEh4u2NwW6rEx/DYxKVeRxEPlHiFzswGATjnph3+/lgOfV4qF82eExEJb8nVE7jltNY8NnEZXy/fxqnt9N/raFLaM3TfAs7Mqjjn8g69L+Xz5j+v4WSVzKHZcyN6pGn2nIhIGLuifzPem72BRz9bxsBWKcTH6kZctCit0J0C4C9zAKeiRRJRSbPnREQiQ3ysjwfP6cCv//Ejb01fy/WDWnodSSpIaYWuK/DFoTfOuW9DnkbCkmbPiYhEjpPb1uPUdvV48avVnN89jZQaCV5HkgpQ2rXY54Beh97458xdEvpIEk40e05EJPLcP7Q9BwsKeXqy9nmNFqUVul1A3cPe60/zKKTZcyIikadFSnWuOqE5H83dxE+bdnkdRypAabdcZwL3m1lTYLf/2Agza1XK1zjn3KNBSyeeKipyjJ+XodlzIiIR6HentmL8vE08PGEpY2/sr7sslVxphe63wFvALRRfyXPACP+rJA5QoaskZq7JJmPXAe46q53XUUREJEA1EuO4c0g77hz3E58u3My53VK9jiQhVOItV+fcOufcyUAikEbxLdffAY1LeTUJcV6pQJo9JyIS2S7omUbn1Fr86fPl5ORpu/XK7JgDapxzBc65zcDDwFTnXEZpr9BHloqwNzefzxdvYVjXRpo9JyISoXw+Y/TwDmzdk8sr36Z7HUdCqMwTB51zDzvnFgOYWYqZ9TazXmaWErp44pVJi7Zq9pyISCXQs2kdzuvWiNemrWHjjhyv40iIBDRC2sz6m9lMYCvFiyZmAVvNbLqZ9QtFQPGGZs+JiFQed53Vjhgznvh8mddRJETKXOj8he1roB3wCsXP0/3e/88dgG/MrG8oQkrFWpe1n9nrdmj2nIhIJdGwVhV+e0pLJi3eyvT0LK/jSAgEcoXuMWAL0M45d7Nz7mXn3F+dczdTXPK2+D8jEW78vE2YwfndtSJKRKSyuPbEFqTVrsIjE5ZSUFjkdRwJskAKXV/gNefc1iNP+I+97v+MRLCiIse4eRkMbFWXhrWqeB1HRESCJDEuhvvObs/yrXt5/8eNXseRIAuk0Dn/qySq+5XAodlzWgwhIlL5nNmpAf1bJPPslyvYlZPndRwJokAK3Y/ADWZW98gT/mM3ALODFUy8MXbuJmokxDKkYwOvo4iISJCZGQ8O68CeA/k8P2WV13EkiErbKeJIDwJfASvM7J/ACv/xdsDlQFX//0qEOjR77vzuaZo9JyJSSbVvWJNL+zblnZnrubhPE9o2qOF1JAmCQObQ/QAMBtZSvB3Yy/7X74F0YLBzbnooQkrF0Ow5EZHocPsZbaieEMsjny3BudKeppJIEdAcOufcNOdcL6Ah0N//auic6+Oc+y4UAaXifL18O6lJVejRRLPnREQqs9rV4rn9jDb8sDqbL5du8zqOBEGZCp2ZVTWzuWZ2I4Bzbptzbpb/pV8JlUBhkWPGmmwGtEzW7DkRkShwad8mtKlfnccnLiM3v9DrOFJOZSp0zrkcoAVayVppLduyh90H8jmh1f+seRERkUooNsbHQ8M6smFHDm98v9brOFJOgdxy/QY4OUQ5xGOHJof3b5nscRIREakoJ7Sqy+AO9fnrN6vZtifX6zhSDoEUut8DXc3sOTNrY2aBrJCVMDc9PZuWKdWoXzPR6ygiIlKB7h/agYJCx1NfLPc6ipRDIIVuLdCG4mK3DDhoZnlHvA6GJKWEVH5hEbPX7mBAS91uFRGJNk2Sq3Ltic0ZPy+D+Rt2eh1HjlMgV9n+Rek7RQSNmV0IjAbaA32cc3P8x88AngTigTzgj865r/3nLgbu9WfcDFzmnNMOxGXw06Zd5OQVMkC3W0VEotJvTmnF2LmbGD1hKR/fNACfT4vjIk2ZC51z7tchzHGkxcAI4LUjjmcBw5xzm82sEzAZSPXf/n0B6OCcyzKzPwM3U1wK5Rh+WJ0NQL8WKnQiItGoekIsd5/Vjts/XMj4+RmaRxqBAppDV1Gcc8uccyuOcny+c26z/+0SINHMEgDzv6pZ8cyNmhRfpZMymJ6eRYeGNaldLd7rKCIi4pHzuqXSrXEST32xnH0HC7yOIwEKqNCZWZKZPW5mC81st/+10H+sdqhClmAkMN85d9A5lw/cBCyiuMh1AN6o4DwRKTe/kHnrd3FCK12dExGJZj6fMXp4RzL3HuSlr1d7HUcCVOZCZ2YtgZ+Ae4AYYArFe7vG+I/95P9MWb/fFDNbfJTXuWX42o7AU8AN/vdxFBe67kCjw3KW9PXXm9kcM5uTmZlZ1siV0tz1O8krLNKCCBERoVvjJEb2SOPN79eyLmu/13EkAIFcoXsRqE3xnq2dnHMjnXMjnHOdgCFAkv8zZeKcO93/fY58fVLa15lZGvAxcIVzLt1/uJv/e6a74k3pPgQGlPJjv+6c6+Wc65WSklLWyJXS9PQsYnxG7+Z1vI4iIiJh4K4z2xIXYzw2cZnXUSQAgRS6k4DnnXNTjjzhnPsPxWXupGAFOxozSwImAvc453447FQG0MHMDrWzMygerSLHMD09m65ptaieoLGCIiIC9WomcvOprZmybBvTVkb3XaxIEkih2weU9m92m/8z5WZm55vZJqA/MNHMJvtP3Qy0Ah4wswX+Vz3/QomHgWlm9hPFV+yeCEaWymxvbj4/bdqt260iIvILVw9sRtPkqjzy2VLyC7XrZyQIpNC9D1xmZv+zFNLMEoHLgfeCEco597FzLs05l+Ccq++cG+I//phzrppzrtthr+3+c68659o757o454Y557KDkaUym712B4VFTvPnRETkFxJiY7h/aAdWb9/HuzPXex1HyiCQ+2wTKN7LdYGZvQqspHiIbzvgeuAgMMHMfvHsmnNuenCiSrBNT88mPtZHj6YVvUBZRETC3ent63Fi67o895+VDO/aiOTqCV5HklIEUugOf3buef67a4SV8Bnzfybm+KJJqE1Pz6ZX09okxulfkYiI/JKZ8eA5HTjzhe/4v/+s5PHzO3sdSUoRSKG7KmQppMLt2J/Hsi17uGNwG6+jiIhImGpdvwZX9G/K29PXcWnfpnRoVNPrSFKCQLb+ejuUQaRizVxT/Ihhfy2IEBGRUtx6Whv+PT+Dhycs4YPr+1G8IZOEm7Dc+ktCb3p6FtXiY+iSVsvrKCIiEsZqVY3jD4PbMmvtDiYt3up1HCmBCl2Umr46mz7N6xAXo18CIiJSuov7NKFdgxo8PnEZufmFXseRo9Cf5lFoy+4DrMnazwmtdLtVRESOLca/z2vGrgO8Pm2N13HkKFTootCM9EPPz2n+nIiIlE2/FskM7dyQl79dzeZdB7yOI0dQoYtC09OzSaoaR/sGWq0kIiJld8/Z7XAOnpy03OsocgQVuijjnGNGejb9WyTj82mlkoiIlF1a7arccFJLPl24mR/X7fA6jhymzIXOzDqa2Ygjjp1iZl+Z2Vwz+0Pw40mwbdiRQ8auA9ruS0REjsuNJ7WgYa1ERn+6hMIid+wvkAoRyBW6PwPXHnpjZqnAp0AXoArwZzO7IrjxJNh+WK35cyIicvyqxsdyz9ntWbJ5Dx/N2eh1HPELpND1AKYe9v5Sirf16uac6wBMAn4bxGwSAtPTs6hfM4GWKdW8jiIiIhFqWJeG9G5Wm6cnr2BPbr7XcYTACl1tYNth788CpjrnMvzvJwDaRyqMHXp+bkDLupr0LSIix83MeGhYR3bk5PHilFVexxECK3RZQBqAmVUD+gNTDjsfS2B7w0oFW7ltH9n78zSuREREyq1Tai0u6tWYt6avIz1zn9dxol4ghe474Cb/wojngTiKn6E7pC2QcbQvlPAwPT0LQAsiREQkKO4Y0pYqcTE8+tlSr6NEvUAK3b3AAWAscA3wZ+fcKgAziwEu4JfP2EmYmZ6eTZM6VUmrXdXrKCIiUgnUrZ7ALae35tsVmXyzfLvXcaJamQudc24t0A7oBjR3zt1z2OmqwE3An4IbT4KloLCImWuydXVORESC6or+zWiRUo1HP1tKXkGR13GiVkCDhZ1zBc65n5xz6484vtc594lzbl1Q00nQLNm8h725BXp+TkREgio+1scD53RgTdZ+3p6+zus4USugQmdmdczsUTP7wcxWmVl///FkM3vQzNqFJqaU13T//q0DNH9ORESC7JS29TilbQovfrWKzL0HvY4TlQLZKaIxsAC4E6gBtKB4oDDOuWzgYjSHLmxNT8+iTf3qpNRI8DqKiIhUQg+c04ED+YU8M3mF11GiUqA7RSRS/AzdqcCRg8w+8R+XMJNXUMSP63bo6pyIiIRMi5TqXHVCMz6cu5FFm3Z7HSfqBFLozgBedM4tA462edtaoHFQUklQLdi4i9z8Ij0/JyIiIfW701qTXC2ehycswTnt81qRAil01YDS1iRXL2cWCZHp6VmYQb/mKnQiIhI6NRPj+OOQtsxZv5NPF272Ok5UCaTQrQD6lXL+bGBx+eJIKExfnU2nRrWoVTXO6ygiIlLJXdizMZ1Ta/Gnz5eTk1fgdZyoEUihew24zMyuBmL8x5yZ1TCz54CTgZeDnE/KKSevgPkbdzKgla7OiYhI6Pl8xkPDOrB1Ty6vfpvudZyoUea9V51zr5hZR+DvQI7/8FigFsXF8EXn3LvBjyjlMWfdTvILnRZEiIhIhenVrA7DuzbitWlruKxfU+rVTPQ6UqUX6GDhm4ETKC51k4DZwCvAic65W4MfT8preno2sT6jd7PaXkcREZEo8ofBbSgocrysq3QVosxX6A5xzs0AZoQgi4TAjPQsujdJomp8wP+qRUREjlvT5Gpc2DON92Zt4PpBLWiUVMXrSJVaQFfoJLLsPpDPoozd9NftVhER8cDNp7bC4Xjpm9VeR6n0At3661ozm25mW83soJnlHfHSfh9hZNaabIocDND8ORER8UBa7aqM6t2ED3/cyMYdOcf+AjluZb4PZ2bPALcBmym+5borVKEkOKanZ5MY56N7kySvo4iISJT67SmtGDNnIy9+tYqnL+zqdZxKK5AHq66meCHEuc65whDlkSCakZ5N72Z1SIiNOfaHRUREQqBBrUQu69uUt2es4zentKJ53WpeR6qUArnlasAElbnIkLn3ICu27dV2XyIi4rmbTm5JfIyPF6as9DpKpRVIoZsM9A5VEAmumWuyATR/TkREPJdSI4ErBjTlk4WbWbVtr9dxKqVACt3vgd5m9pCZpYYqkATH9PRsaiTE0qlRTa+jiIiIcMOgllSNi+H5Kau8jlIplbnQOee2A+8ADwIbzCxfq1zD1/T0LPq2qENsjCbTiIiI9+pUi+fqgc2ZuGgLSzfv8TpOpRPIKtdHgPsoXuU6B61yDVubduawPjuHK/s38zqKiIjIz64d2IK3pq/juSkr+dsVvbyOU6kEssr1BrTKNSLMSPc/P9dKCyJE5P/bu/PwOMor3+Pfo83yLu/GyLss1gAOBmyzY7aQECAJgcwQIBuBEAYCmbnhJrk3mSQzWQazBEhCMhm2ZCDrJYQwYGOzWTJgYhvigJf2ghcwlrzvlnTuH1UNTdOSWlJLVSX9Ps/TT0tvvVV13uq2ffwuVSLxMbBPKV84eQIzZy3jlXVbOapSt9UqlLaMx5WjVa6JUJuqZ0jfMqqH9486FBERkff4zInjqOhTysxZWvFaUuP7tgAAIABJREFUSG1J6J5Eq1xjz92pSdUzdeIQioos6nBERETeo395KV88ZSJPL93Ey2u2RB1Ot9GWhO7LwLFm9m2tco2vVXW7eGv7Xj3uS0REYuuK6WMZ2q+MmbOWRh1Kt9GWhG4dcCTwDbTKNbbmpXT/ORERibc+ZSVcfepE5q2of+e+qdIxbVkU8SvAOysQKYzaVB0HDSxn3JA+UYciIiLSrMumjuXnz61k5pPLePiLUzHTNKGOyDuhc/crOzEOKYCmJqc2Vc8Zh47QHwwREYm18tJirj29iv/zyBKeX1HHyZOGRR1Soumus93I62/tYMvuA5o/JyIiiXDJcaMZNbCcW55chrsGATui2R46MzsFwN2fzfy9Nen60vVqUnUATFNCJyIiCdCrpJjrZkzi5j+8ytylb3PGoSOiDimxWhpyfRpwM+vt7vvTv7dQ38LtxQWLTtqkNlXP+KF9GVXRO+pQRERE8vKJYyv5ydMpZs5axumHDNeUoXZqKaE7HSBM5t75XeKpobGJF1Zt5qPHjIo6FBERkbyVFhfxTzMm8dXfLuaJJRs598iRUYeUSM0mdO7+TEu/S7y8sn4bO/c1aP6ciIgkzoXHjOLuuSu4ddYyzj58hG6M3w55L4owszlmNqOF7aeb2ZzChCVtlX5+67QJSuhERCRZSoqLuP7MSSzduIPHXn0z6nASqS2rXE8DWpqtOBw4tUPRSLvVpOo4dGR/hvTrFXUoIiIibXb+UaOoHtGP22Yvo7FJK17bqq23LWnpCk8EdnYgFmmnvQcaWbB6i54OISIiiVVUZHzlzGpSm3bxyKL1UYeTOC3eWNjMPg18OqPoZjP7TI6qFcBkYFYBY5M8LXxjK/samjR/TkREEu2cI0Zy+EEDuP2p5Zx/9ChKi3W73Hy1dqUGA5PClwMjM35Pv6qA3gSPBruq0yKVZtWm6igyOH7C4KhDERERabeiIuPGs6pZU7+bP/x1XdThJEqLPXTufjtwO4CZNQE3uPuvuyIwyd+8VD0fqKxgQHlp1KGIiIh0yIzDhnP06ArueGoFF02upKxEvXT5yPsquXuRkrn42bmvgcVrt3KihltFRKQbMAt66dZv3cPDC9ZGHU5iKO1NuJdWb6ahybUgQkREuo1TJg1lythB3DVnBXsPNEYdTiIooUu42lQ9ZcVFHDt2UNShiIiIFISZcePZ1by1fS+/fuGNqMNJBCV0CVeTqmPymAp6l+kRuiIi0n1MnziUaROGcPfTKfbsVy9da5TQJdjW3ftZsmG7hltFRKRbuunsaup27uP+2tVRhxJ7SugSbP7KetxhepUWRIiISPczZdxgTqkexk+fSbFzX0PU4cSaEroEq0nV06esmKMrK6IORUREpFPceFY1W3Yf4N55q6IOJdbySujMrMLMjjeziS3UGW9mlxcuNGlNTaqe48YN1j16RESk2zpmdAVnHjace55dybY9B6IOJ7ZazQTM7HvARqAWWGZmL5rZB3NUnQ78V4Hjk2a8vX0vK97eqcd9iYhIt/eVs6rZvreB/3xevXTNaTGhM7NPAjcDzwDXAv8GjAFqzOyyzg9PmlO7sh5ACyJERKTbO2LUQD505Eh++fwqtuzaH3U4sdRaD90NwFx3P9vdf+ru3wQOA54F7jWzGzo9QsmpZkU9A8pLOHzUgKhDERER6XRfOauaXfsbuOe5lVGHEkutJXSHAr/PLHD3LcC5wP3ALWb2nU6KTVpQs7KOqROGUFxkUYciIiLS6apH9Of8o0Zx77zV1O3cF3U4sdNaQteUq9Ddm9z9s8BtwNfN7M48jiUFsnbzbtZu3qP5cyIi0qNcf+Yk9jU08tOnU1GHEjutJWHLgJOa2+juNwHfBr4EqKeui9Sk6gA4sUrz50REpOeYOKwfF02u5IH5a9i4fW/U4cRKawnd48AFZtZsV5C7fxv4CjC6UEGZ2cVmtsTMmsxsSkb5WWb2spm9Gr6fkbHtEjN7Jdzvh4WKJY5qUvUM7deLquH9og5FRESkS10/YxKNTc7dc1dEHUqstJbQ/RL4F2B4S5Xc/Xbg48C/FiiuvwEfI1h8kakOON/dPwBcATwAECacPwJmuPsRwAgzm1GgWGLF3alJ1TN94hDMNH9ORER6ljFD+nDxlEr++8W1rN+6J+pwYqPFhM7d17v7Xe7+WmsHcvf/F/bWdZi7v+buS3OUL3T3DeGvS4ByM+sFTACWufumcNtsggSz20lt2smmHfs0f05ERHqsL58xCYA756iXLq3dCxnMrNzMLjezEYUMqA0+Dix0933ACuBQMxtnZiXAhRRwCDhOalK6/5yIiPRsB1f05tLjR/PbBWt5o3531OHEQkdWpg4keDLEEe3Z2cxmm9nfcrwuyGPfI4AfAF+Ed26lcg3wMPAcsBpo9im+ZnaVmS0wswWbNm1qrloszVtRx8EVvRk9uHfUoYiIiETm2tOrKC4y7pizPOpQYqGjtxpp9yQudz/T3Y/M8XqkxROaVQJ/BC5393fWLbv7o+5+grtPA5YCzX7C7n6Pu09x9ynDhg1rbxO6XGOTM3/lZk6s0vw5ERHp2UYMKOeyqWP5w1/XsXLTzqjDiVxHEzovSBR5MrMK4DHgZnefl7VtePg+iOA2Kr/oyti6wmtvbmfbngMabhUREQGuOW0ivUqKuf0p9dJF1kPX4kHNLjKzdcA04DEzeyLc9GWgCvimmS0KX+kVuLeb2d+BecD33X1ZZ8QWpfT956ZpQYSIiAhD+/Xiiunj+NPiDSzbuCPqcCLVkYRuEzCeIIEqKHf/o7tXunsvdx/h7ueE5d91977ufkzG6+1w26fc/fDw9VChY4qDmlQ9E4f1ZcSA8qhDERERiYUvnjKBvmUl3Da72/XjtEm7E7rw8V9rwlWm0skONDbx4qrNGm4VERHJMKhvGZ89cRx/efUtlmzYFnU4kdHzVxNi8dqt7N7fqPvPiYiIZPncyRMYUF7CrbN67lw6JXQJUZOqxwymTlBCJyIikmlg71K+cPIEZr+2kcVrt0YdTiSU0CVETaqOww8awKC+ZVGHIiIiEjufOWk8g/qUMnNWz5xLp4QuAfYeaOSva7ZquFVERKQZ/XqV8MVTJ/LMsk28vGZz1OF0OSV0CfDymi3sb2zSgggREZEWXD5tLEP7lXHLkz2vl04JXQLUpOooLjKOGz846lBERERiq09ZCdecVkVNqp7a8NnnPYUSugSoSdVzdOVA+vUqiToUERGRWPvHE8YwYkAvZs5ainuXPtAqUkroYm7H3gO8sm4bJ1ZpuFVERKQ15aXFfPn0Kl5avYXnltdFHU6XUUIXcy+u2kxjk+txXyIiInn65HGjObiiN7fMWtZjeumU0MVcTaqespIiPjhmUNShiIiIJEKvkmKuO6OKxWu3Muf1t6MOp0sooYu5mlQ9U8YOory0OOpQREREEuPjx1YyZnAfZvaQXjoldDG2edd+Xntzu+4/JyIi0kalxUVcP2MSSzZs54klb0UdTqdTQhdj81cGS66n6f5zIiIibXbh5IOZMKwvt85aTlNT9+6lU0IXY/NW1NG3rJijKgdGHYqIiEjiFBcZN5xZzdKNO/jzq29GHU6nUkIXY7Wpek6YMITSYn1MIiIi7fGRDxzEISP6c9vsZTQ0NkUdTqdRphBTb27bw8q6XZo/JyIi0gFFRcZXzprEyk27eGTRhqjD6TRK6GIq/cgS3X9ORESkY845YiRHjBrA7U8t50A37aVTQhdTNal6KvqUctjIAVGHIiIikmhmxo1nVfPG5t38/uV1UYfTKZTQxZC7U5uqZ9qEIRQVWdThiIiIJN4Zhw7nmNEV/HjOCvY1NEYdTsEpoYuhNfW7Wb91j+bPiYiIFEi6l2791j385qW1UYdTcEroYqgmnD83vUr3nxMRESmUkycN5bhxg7hz7gr2HuhevXRK6GKoJlXHiAG9mDC0b9ShiIiIdBtBL90hbNy+j1+98EbU4RSUErqYSc+fmz5xKGaaPyciIlJI0yYOYfrEIfzk6RXs3t8QdTgFo4QuZpZt3En9rv26XYmIiEgnuensaup27uf+2jVRh1IwSuhipiZVB6AFESIiIp3k2LGDObV6GD97JsXOfd2jl04JXczUpOoZM7gPlYP6RB2KiIhIt3XjWdVs2X2A/3p+VdShFIQSuhhpaGxi/sp6TqxS75yIiEhnOnp0BWceNoKfP7eSbXsORB1Ohymhi5ElG7azY28D0ybqdiUiIiKd7cazqtm+t4H/fG5l1KF0mBK6GEnff27aBPXQiYiIdLbDRw3gvA+M5JfzVrNl1/6ow+kQJXQxUpOqo3pEP4b17xV1KCIiIj3CDWdWs2t/Az97Ntm9dEroYmJ/QxMvrd7MdA23ioiIdJnqEf356NGjuK9mNZt27Is6nHZTQhcTi9ZuZe+BJt1/TkREpItdP2MS+xoa+ekzqahDaTcldDExb0UdRQZTNX9ORESkS00Y1o+PfbCSB+evYeP2vVGH0y5K6GKiNlXPkQcPZGDv0qhDERER6XGunzGJxibnrrkrog6lXZTQxcDu/Q0sXLtFw60iIiIRGT24DxdPGc1DL65l/dY9UYfTZkroYmDB6i0caHQtiBAREYnQdWdUAXDnnOURR9J2SuhioCZVT0mRcdy4QVGHIiIi0mONqujNp44fzW8XrOON+t1Rh9MmSuhioDZVx+QxFfQpK4k6FBERkR7t2tOrKC4ybn8qWb10Sugitm3PAV5dv03DrSIiIjEwfEA5n546lj8uXEdq086ow8mbErqIvbCyniaH6VoQISIiEgtXnzaR8tJibp+dnF46JXQRq0nVU15axDFjKqIORURERICh/XpxxfRxPPrKBpa+tSPqcPKihC5ital6jhs3mF4lxVGHIiIiIqGrTp5A37ISbpu9LOpQ8qKELkKbduxj6cYduv+ciIhIzAzqW8ZnTxrP4397iyUbtkUdTquU0EVo/sp6AC2IEBERiaHPnTSeAeUl3Dor/nPplNBFqCZVT/9eJRw5akDUoYiIiEiWgb1LueqUCcx+bSOL126NOpwWKaGLUE2qjhMmDKGkWB+DiIhIHF154ngG9Sll5qx4z6VTJhGRdVt2s6Z+t25XIiIiEmP9epVw9akTeWbZJl5esznqcJqlhC4italw/lyVEjoREZE4u3zaOIb268UtT8a3l04JXURqU/UM6VtG9fD+UYciIiIiLehdVsyXTptITar+nQ6ZuFFCFwF3pyZVz9SJQygqsqjDERERkVb8wwljGDmgnJmzluLuUYfzPkroIrCqbhdvbd+r+XMiIiIJUV5azLVnVPHS6i08t7wu6nDeRwldBOaF3bUn6v5zIiIiiXHJlNEcXNGbW2Yti10vnRK6CNSm6hg1sJyxQ/pEHYqIiIjkqaykiH+aUcXitVuZ8/rbUYfzHkroulhTk1ObqmfaxKGYaf6ciIhIknzsg5WMHdKHmTHrpVNC18Vef2sHW3Yf0Pw5ERGRBCotLuL6GZNYsmE7Tyx5K+pw3qGErovVpIKJlNOU0ImIiCTSBccczMRhfbl11nKamuLRS6eErovVpuoZP7Qvoyp6Rx2KiIiItENxkXHDmdUs3biDP7/6ZtThAEroulRDYxMvrNqs4VYREZGE+/AHDuLQkf25bfYyGhqbog5HCV1XemX9Nnbua2C6blciIiKSaEVhL93KTbt4ZNGGqMNRQteV0o8LmTphcMSRiIiISEedc8QIjjx4ALc/tZwDEffSKaHrQjWpOg4d2Z8h/XpFHYqIiIh0kJlx41nVvLF5N79/eV2ksSih6yJ7DzSyYPUWDbeKiIh0I6cfMpxjRlfw4zkr2NfQGFkcSui6yMI3trKvoUkLIkRERLoRM+Oms6sZ0q+MTTv2RRZHSWRn7mFqUnUUFxknaP6ciIhIt3JS1VBOqor2CVBK6LpITaqeDxw8kP7lpVGHIiIiIgUUh0d5asi1C+zc18DitVs13CoiIiKdQgldF3hp9WYamlwLIkRERKRTKKHrArWpesqKizh27KCoQxEREZFuKJYJnZldbGZLzKzJzKZklB9vZovC12Izuyhj27Fm9qqZrTCzOywOA9qhmlQdk8dU0LusOOpQREREpBuKZUIH/A34GPBsjvIp7n4McC7wMzNLL+z4CXAVMCl8ndtFsbZo6+79LNmwXcOtIiIi0mlimdC5+2vuvjRH+W53bwh/LQccwMwOAga4e627O3A/cGGXBdyC+SvrcYcTq7QgQkRERDpHLBO6lpjZCWa2BHgVuDpM8A4GMp+5sS4si1xNqp4+ZcUcVVkRdSgiIiLSTUV2Hzozmw2MzLHp6+7+SHP7ufsLwBFmdhhwn5k9DuSaL+ctnPsqguFZxowZ06a426p3aTHnHDGSspLE5c4iIiKSEJEldO5+Zgf3f83MdgFHEvTIVWZsrgQ2tLDvPcA9AFOmTGk28SuEm887rDMPLyIiIpKsIVczG59eBGFmY4FDgNXu/iaww8ymhqtbLwea7eUTERER6U5imdCZ2UVmtg6YBjxmZk+Em04CFpvZIuCPwJfcvS7cdg3wC2AFkAIe7+KwRURERCJhwaLQnmvKlCm+YMGCqMMQERERaZWZvezuU7LLY9lDJyIiIiL5U0InIiIiknBK6EREREQSTgmdiIiISMIpoRMRERFJOCV0IiIiIgmnhE5EREQk4ZTQiYiIiCScEjoRERGRhFNCJyIiIpJwSuhEREREEk4JnYiIiEjCKaETERERSTgldCIiIiIJp4ROREREJOHM3aOOIVJmtglYU4BDDQXqCnCcpFL71X61v+dS+9V+tb/rjHX3YdmFPT6hKxQzW+DuU6KOIypqv9qv9qv9UccRFbVf7Y9D+zXkKiIiIpJwSuhEREREEk4JXeHcE3UAEVP7eza1v2dT+3s2tT8GNIdOREREJOHUQyciIiKScD0uoTOzc81sqZmtMLOv5dhuZnZHuP0VM/tga/ua2WAzm2Vmy8P3QRnbbg7rLzWzczLKjzWzV8Ntd5iZZcXxCTNzMyvoypm4t9/MrjSzTWa2KHx9vie1P9z2STP7u5ktMbNf96T2m9mtGZ/9MjPb2sPaP8bM5prZwvD85/Ww9o81s6fCcz9tZpXdtP3fM7O1ZrYz6/y9zOzhcJ8XzGxcD2r7KWb2VzNrMLNPFKrdCWr/jRb8vf9K+GdgbJsb6e495gUUAylgAlAGLAYOz6pzHvA4YMBU4IXW9gV+CHwt/PlrwA/Cnw8P6/UCxof7F4fbXgSmhed5HPhQRgz9gWeB+cCUntR+4Ergzp76+QOTgIXAoPD34T2p/VmxXAf8sie1n2AuzjUZ+6/uYe3/LXBF+PMZwAPdtP1TgYOAnVnn/xLw0/DnS4GHe1DbxwFHAfcDnyjU556g9p8O9Al/vqY9n31P66E7Hljh7ivdfT/wEHBBVp0LgPs9MB+oMLODWtn3AuC+8Of7gAszyh9y933uvgpYARwfHm+Au9d68Ondn7EPwHcIvih7C9d0aKUNaXFof2dJQvu/ANzl7lsA3P3tHtb+TJ8C/rsA7U5LQvsdGBD+PBDYULDWJ6P9hwNPhT/PzRFfR8Si/QDuPt/d38wRY+axfgfMSPdedlDs2+7uq939FaCpAO3NloT2z3X33eGv84E29073tITuYGBtxu/rwrJ86rS074j0BxS+D8/jWOtyHcvMJgOj3f3PbWlYnmLf/tDHw27n35nZ6PyalpcktL8aqDazeWY238zOzbt1rUtC+4Fg6I3gf7Zz8mhXvpLQ/m8Bl5nZOuAvBL2UhZKE9i8GPh7+fBHQ38yG5NG2fMSl/XnF6O4NwDagEO1PQts7U9La/zmC3sI26WkJXa7/6WQv822uTj775nu+nOVmVgTcCtzUynHbK9btD98fBca5+1HAbN79308hJKH9JQTDrqcR9FD9wswqWjlPvpLQ/rRLgd+5e2Mr52iLJLT/U8C97l5JMAT0QPj3QiEkof1fBU41s4XAqcB6oKGV8+QrLu0v9D75SELbO1Ni2m9mlwFTgB/lUz9TT0vo1gGZPT6VvH9Io7k6Le27MeyaJXxPD5O1dKzKHOX9gSOBp81sNcFY+5+scAsj4t5+3L3e3feF5T8Hjs2zbfmIffvDbY+4+4Gwq34pQYJXCElof9qlFHa4taV48qnTVe3/HPAbAHevBcoJnhNZCLFvv7tvcPePuftk4Oth2bb8m9iiuLQ/rxjNrIRg2H1zK/vkIwlt70yJaL+ZnUnwvf9oxr+D+fMCTjyM+4ug92MlwVBOenLjEVl1Psx7J0a+2Nq+BJl05sTIH4Y/H8F7J0au5N2JkS+Fx09PCj4vR7xPU9hFEbFvP3BQRiwXAfN7WPvPBe4Lfx5K0G0/pKe0P9x2CLAagvtk9rDP/3HgyvDnwwj+ESjIdUhI+4cCReHP3wP+tTt+/hnny54Yfy3vXRTxm57S9ozyeyn8oojYtx+YTLB4YlK721nIi5aEF8EwxrLwwn09LLsauDr82YC7wu2vkpFQ5do3LB9CMJF3efg+OGPb18P6S3nvStYpwN/CbXeS4y9tCpzQJaH9wL8DS8I/DHOBQ3tY+w2YCfw9PP+lPan94bZvAd/viX/+CRYFzAu//4uAs3tY+z8RHmcZ8AugVzdt/w8JenGawvdvheXlBCt9VxCsBJ7Qg9p+XPj7LqAeWNLDPvvZwEaCP/eLgD+1tY16UoSIiIhIwvW0OXQiIiIi3Y4SOhEREZGEU0InIiIiknBK6EREREQSTgmdiIiISMIpoRORZpnZ583MzazNzxVMMjObaGZ/MbMtYfu/HHVMrTGz581sdtRxtIWZlYTX9xtRxyKSdEroRBLAzP5gZgfMbFgLda4P/3E8vytj66buI7hX2reATwNPRhqNiEgrlNCJJMMDBHcsv7SFOpcBdcD/FPC8/wX0dvd1rdbsJsysHJgO3O/ut7v7g+6+LOq4RERaooROJBkeI3im42W5NppZNUGP0kPufqCjJzOzvgDu3ujuezt6vIQZTnDX+K3t2dkCfQobkohIy5TQiSSAu+8neGj78WY2KUeVT4fvD6QLzOxCM3vUzNaZ2f7w/S4zG5C5o5l9NxyqPdzM7jWzOoJnqeacQ2dmp5rZw2a2xsz2mdlGM7vfzEZlHTe978lm9gMze8vM9pjZE2Y2LrsBZlZlZg+G9faZ2Soz+1k6uQzrDDCz/zCz1WGdN8Jjl+dzHcNY5prZTjPbYWazzOyEzGsBrAl//U4Yf0MLx0vPAfupmV1sZouBvcBnw+2fNbPZGW1aaWbfM7OyrOM8aGZ7zWykmf0ujG1L2P7yrLpmZjeb2drwetZktiGr7pAwtjfD879mZjeamTXTho+a2eLwuIvM7JSwzofNbGEY42tmdnae1/tiM3vRzLaZ2S4zW25md+ax31gz+7WZ1YXnXGRml2fVqQrj/pqZfcnMUmHdv1rwkPPsY3bouyMSdyVRByAieXuA4NmD/0gwtyvTPwDL3P3FjLLPAfuBHwNbCB7+/HngSODUHMf/DfAG8H+Avjm2p10CVAD3AG8DhwBXASeY2dE5evRuBfYA/wYMA74atuXkdAUzOxJ4juA/mfcQPBuxEvg4MAjYZWa9CZ5vPD6skwKOAW4keBj2R1qIGTM7HXgCWAt8NzzX1cAzZnaau88neI5mXRjz74BHCJ652JqTgYsJngV5N/BaWH4dwTNLHyd4RuWJBA/xrgSuyDpGURjfK8A/Ezwg/CqCa/zNjHr/CnwDmAX8ieD6P07Qo7gyo73lBNfrMOAnwOsEDyC/BRgD3JB1/mnA+WH8e4B/AR4zsy+E+9wN7A7Lf2dmo919W3MXJEz6HgbmAP8baAAmAOc0t0+433CghuA79mNgA8FUg/vMbLC735a1y6UEvap3E3zfvxjGfZq714bH7NB3RyQRCvnwW7300qtzXwQP7V6eVXYi4MA3ssr75Nj/yrDuCRll3w3Lfp+j/ufDbZWtHPe0sN4lOfZ9BijKKP9qWH5IRtmzBMlCdY5jp585/U2CROPwrO1fCo93eivXbhHBQ7+HZZQdDOwAajLKxuW6ns0csySs2wQck2N7rmv1LaAROCij7MHwON/Nqvtn4M2M34cTJC2zsq5p+hrMzii7ISz7XOa1BP4QxntIVhv2Z15/guTOCXocq3KUf76Va/NjgmkCxXlcv29klN0Wls3IKCsDXiBIiivCsqqMuCdm1B0BbAeeyyjr0HdHL72S8NKQq0iyPAhUmdnUjLLLCP5R+lVmRXffDe8M0Q00s6HA8+HmY3Mc+yf5BJA+bnjs/uFx/0aQGOU67s/cPbOX65nwfUJ4jBEEPVz3eo7FB+7u4Y+XAPOAt81saPoFpG/VcUZzMVswZHw0cJ+7b8o49nrg18A0MxvSQrNbU+Pui3LEnv4MisysIoz3WYLeuMk5jnN31u/PACPt3Tl55wClwB1Z1/QXBNc/00cIehvvzYjHgR8RJHYfzqo/N+v616ZjcPcVOcon5Ig/01agP5DX8GyGjwAL3f2pjLj3EyR6fXj/5/you6cy6m4E/hs40cwqwuJ2f3dEkkIJnUiyPBi+XwYQzsX6JPC8u6/KrGjBnLhHgZ0E/7huIhjKhGA4K1sqR9n7mFmlmf3KzLYR9IRsCl/9mznumqzft4Tvg8P3qvD91VZOXQ3MyDhf+rU03D68hX3Hhe+v59j296w67ZHz2pnZdDObS9D7uIUg3nSikn2tDrj7hqyy9LUaFL6PDd+XZlYKE57VWfuOI+jNbcwqT7d3fFb5G1m/b22lfBAtuzOM8y/hHMJfm9klZtbaVJ+xtPw5Zce9NLtiWGa8e7068t0RSQTNoRNJEHdfYWa1wCVmdgNwHkFi9EBmvbBn4hmCIapvEgzV7iYYunqM3P+Z29Pa+c2smKBXYyjwA4J/ZHcR9BD+tpnjZicU7xwu692bqZdZfw7w781sX9/K/i0dN5/zt+Sba80ZAAAD4klEQVR9187MJhIkb8sJ5mq9QTB8OQb4T95/rVqaq5fPtbIcZS3JPkZzn1Nrn1/ug7tvNLPJBInUOQQ9dZ8CFpjZKe7e6vetmfNlx53Pteis745IbCihE0meBwiG5s4h6KnbR5BMZZpBkHRd6O7z0oVmdngHz30MwST8y9z9nSFeC1aiDmznMdPDeUe1Ui8F9Hf39jwNYXX4fmiObemy7J7EjroQKAc+FA7tAmBm53XgmKvD90N5t7c13VM7FtiYVfcoMyvO6qU7LOtYncaDW+j8T/jCzK4D7gA+RtYUgQxraPlzWt1MeaZqgkQv/Zl25LsjkggachVJnocJJoJfSzDf6FF3z75nWrq3J/vP+D938NwtHbetPUQAuPtbBCtcr7TgfnrvkXGLjYeA48zsoznqlJtZ/xbOsY5gUcTl4dyp9H6jCFYN17h7fXvib8H7rlXYw3ljB475JMFq0evMLPMz+DzBkHemRwlWFb9zu4/wWqYXpTzWgTha1cycxIXhe66h+bRHgcnhquT0sUqB6wl6medk1T8/7A1N1x1B0BNYk/Hnot3fHZGkUA+dSMK4+2Yz+wtBDxBkDbeGniNYYfgrM/sxwVDf+QS9dh2xhKBH7TYL7iVXR7DC9Xjene/VHl8miPklM7sHWAaMIrhtyXnAOuCHBBP5/2BmDwIvEgwhH0Iwj/AC3l30kcuNBLcFmW9mPydIQK8hWGRwUwdib87jwPeBx8M2lRBMzi9t7wHDYcwfATcD/2NmjxD2mAKrsqrfQ5Do3RMOfS4luJbnAbe7e665Z4V0bzj0P4dguHk4wW1idhIkbc35d4LP89Hwu7uB4LpNBb6S4z8vfweeM7O7gAMEty0pB/5XRp2OfndEYk8JnUgyPUCQ0NUTJA7v4e514dDefwD/lyCh+wtBb82b7T2pu+83s48Q3KctnQQ9TbBK8LkOHPcVC26O+23gM0A/gn/InyBMFN19t5mdRvAP9SUEvTA7Ce69dgdBstnSOeaGN5z9NsG99hyYT3Crlfntjb2F871uZhcS3Bbm+8A2gnv9/ZKgt7C9vk4wb/EagtXBC4EPEdwrLvP8e8Lr9T2CpGUwQdL3VWBmB86fr/sJbrB8VXjuOoKVpt9x9+yFFu9w97fNbDrBNbuK4J6IS4Er3P3+HLs8RLA45yaC+/v9HfhI5lSDjn53RJIgfX8nERGRxDCzKoJ5hDe7+/ejjkckappDJyIiIpJwSuhEREREEk4JnYiIiEjCaQ6diIiISMKph05EREQk4ZTQiYiIiCScEjoRERGRhFNCJyIiIpJwSuhEREREEk4JnYiIiEjC/X9OA23FBjiUaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "re = mdf.cov_re.iloc[1, 1]\n",
    "likev = mdf.profile_re(1, 're', dist_low=.5*re, dist_high=0.8*re)\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "plt.plot(likev[:,0], 2*likev[:,1])\n",
    "plt.xlabel(\"Variance of random slope\", size=17)\n",
    "plt.ylabel(\"-2 times profile log likelihood\", size=17)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
