{
 "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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXRV1fnG8e+bGUggQEKYp5AwTxJEZgUVFRFQsda5znVurdZa7WC11dZfrbOitVat1llUrAMOgKBgmGdCwjwmDCEQEpLc/fsjF5tiCLlwk5Pc+3zWOiu559xzeeJimZe9z7u3OecQERERkdAX4XUAEREREakdKvxEREREwoQKPxEREZEwocJPREREJEyo8BMREREJEyr8RERERMJElNcB6oOkpCTXsWNHr2OIiIiIHNW8efPynHPJlV1T4VcNHTt2JDMz0+sYIiIiIkdlZuuPdE1TvSIiIiJhQoWfiIiISJhQ4SciIiISJlT4iYiIiIQJFX4iIiIiYUKFn4iIiEiYUOEnIiIiEiZU+ImIiIiECRV+IiIiImFChZ+IiIhImFDhJyIiIhImVPiJiIiIhAkVfiIiIiJhQoWfiIiISJhQ4SciIiISJlT4iYiIiIQJFX4iIh5yzlFS5vM6hoiEiSivA4iIhJv8whJmZecxY3UuM7Py2La3iL5tmzC0SxJDUpM4oUMisVGRXscUkRCkwk9EpIaVlvlYuHEPM7LymJmVy6KNe/A5SIiNYkiX5ozt04q5a3fx5JdrePyLNcRGRTCwYzOGdGnO0NQkerVpQmSEef1jiEgIUOEnIlIDNu4qZEZWLjNW5zJ7zU4KikuJMOjTNpGbRqUxIi2Jfu0SiYr87xM3e4tKmJOzi9nZecxes5M/f7wKWEVCXBQndW7O0NTmDOmSRFqLeMxUCIpI4FT4iYgEwb7iUr7N3smMrPLp27V5+wFo3SSOsX1aMTwtmaFdmpPYMOaIn9E4LprTeqRwWo8UAHILipmdncc32TuZlZ3HZ8u3A5CcEMuQ1PLRwMGpzWnXrGHN/4AiEhLMOed1hjovIyPDZWZmeh1DROoQn8+xbMve70f15m/YTUmZo0F0JCd1bsaI9GSGpyWTmtwoaKNzG3cVMjs7j1lrdjI7eyd5+4oBaN+sIUO7NGeIvxBMio8Nyp8nIvWTmc1zzmVUek2F39Gp8BMRgO17i75vyPh6TR679h8EoEerxoxIT2ZEWhIDOjatlcYM5xxZO/Yxa015ITgnp3w6GaBbywSGpCYxJLU5gzo3IyEuusbziEjdocLvOKnwEwlPRSVlzF276/tib9X2AgCS4mMZkZbE8PQkhnVJJjnB+xG20jIfS7fsZdaaPGZn55G5bjfFpT4iI4w+bZt8PzV8QoemxEWrY1gklKnwO04q/ETCg3OO1dv3MWN1LjOycpm7dhfFpT5iIiMY2Kkpw9OSGZGWTLeWCUTU8S7bopIy5m/Yzew15c8HLt6UT5nPERsVQUbHpt+PCPZu0+R/GkxEpP5T4XecVPiJhK5d+w8y09+QMTMrl+17y5+b69IinhFpyQxPT+KkTs1pEFO/R8kKikqYu3aX//nAPFZuKx+9TIiNYlDn5uUjgl2SSE9Rx7BIfVdV4aeuXhEJKwdLfczfsJuZWbnMWJ3H0i35OAdNGkQzLC2pfAo3LZnWiQ28jhpUCXHRjO6ewuju5R3DefuK+Sa7vElkdnYe01aUdwwnxccwODWJof5CUB3DIqFFI37VoBE/kfrLOce6nYX+5/Ry+SZ7J/sPlhEZYZzQPrF8+jY9md5hvkjypt2FzPaPBs7K3kluQfnIZ7tmDRjSOYkh/q7huvA8o4hUTVO9x0mFn0j9sreohNlr8piRVb4t2qbdB4DyImaEv9AbnNqcxup2rZRzjjU79jE7eyez1uTxbc5O9haVdwynp8QzJDWJoV2SGNS5mf4bitRBKvyOkwo/kbqtzOdYtGkPM1fnMSMrl4Ub91Dmc8THRjE4tfn307cdkxp5HbVeKvM5lm7O/35a+Lt1uygq8RFh0Ltt4vfTwgPUMSxSJ6jwO04q/ETqni17Dnzffft1Vh57i0oxgz5tmnw/fdu/fSLR6lgNuuLSMhZs2MPsNeXTwos27qHU54iJiuCCjLbcM7aHCkARD6m5Q0RChnOOyTNyePDjlTgHKY1jGdOzJSPSkxnaJYlmjY68JZoER2xUJCd1bs5JnZvzc8q3q/tu7S4+Xb6NV77dwLz1e3jyov50To73OqqIHEYjftWgET+RuqGkzMdvpizltbkbGdu7FbeemkZaCy0/Upd8uXIHP3tjIaVljgfP683ZfVp7HUkk7FQ14qc5EBGpF/YWlXDli9/x2tyN3HhKKo//uD/pKQkq+uqYU7q1YOotw0lPieemVxdw73tLKS4t8zqWiPip8BOROm/jrkLOe2o232Tv5M/n9+GOMd3q/M4Z4axNYgNev24w1wzvxMvfrue8p2ezYWeh17FEBBV+IlLHLdiwm4lPzWL73iJeuvJELsho53UkqYboyAh+PbYHky8dwIadhYx9fCYfL93mdSyRsKfCT0TqrI+WbOXCyd/SICaSd24YypAuSV5HkgCd3rMlU28ZTuekRlz/yjx+/8EyDpb6vI4lErZU+IlIneOc45np2dzwr/n0bN2Y924YSpcW6hCtr9o1a8gb1w/miiEd+cesdUx69hs27dbUr4gXVPiJSJ1SUubjV+8s4cH/rOTsPq149ZqTaB6vbcLqu9ioSH53Tk+euvgEcnbsY+xjXzNt+XavY4mEnTpZ+JnZJDNbZmY+M8uocP40M5tnZkv8X0dVuPaVma0ys4X+o8URPvtXZrbG/94xtfHziEj15B8o4Sf/+I5/f7eRm07pwmMX9tdCwCHmrN6t+ODmYbRJbMDVL2Xyx49WUFKmqV+R2lJXF3BeCpwLPHvY+TxgnHNui5n1Aj4B2lS4frFz7ogL7plZD+BCoCfQGphmZunOOa01IOKxjbsKufLF71i3cz9/Ob8Pk9TEEbI6JjXinRuGcP/U5UyekcO89bt5/Mf9aZ3YwOtoIiGvTo74OedWOOdWVXJ+gXNui//lMiDOzAKZAxoP/Ns5V+ycWwusAU48/sQicjz+t3N3kIq+MBAXHcn9E3rz2I/7s3LrXsY+NpMvV+3wOpZIyKuThV81nQcscM4VVzj3D/80771W+aqubYCNFV5v4n9HDEWklh3q3G0YE8U7NwxlcGpzryNJLTqnb2vev3kYKY3j+Mk/vuPPH6+kVFO/IjXGs8LPzKaZ2dJKjvHVuLcn8BBwXYXTFzvnegPD/celld1ayblK96wzs2vNLNPMMnNzc4/+A4lIQJxzPP1VeedurzZNePeGIercDVOpyfG8d+NQLhzYjqe+yuai5+ewfW+R17FEQpJnhZ9z7lTnXK9KjilV3WdmbYF3gcucc9kVPm+z/2sB8CqVT+FuAirOIbUFtlTyPpxzk51zGc65jOTk5MB+OBGpUkmZj7veXsJDH69kXN/W/OvqQercDXNx0ZE8eF4f/npBX5ZsyuesR2fydVae17FEQk69muo1s0RgKvAr59ysCuejzCzJ/300cDblDSKHex+40MxizawTkAbMrfnkInJI/oESrvjHXF7P3MjNo7rw6I/6qXNXvnfuCW15/6ahNGsUw6UvzOGRz1ZT5qt0YkZEjkGdLPzMbKKZbQIGA1PN7BP/pZuALsC9hy3bEgt8YmaLgYXAZuA5/2edY2b3ATjnlgFvAMuBj4Eb1dErUns27irkvKdnM3ftLh6e1JfbT++qPXflB9JSEphy01Am9m/Do59ncdkLc8gtKD76jSJyVOac/iV1NBkZGS4z84irxIhINczfsJtr/plJqc/xzCUD1MQhR+Wc483MTdw7ZSmNG0Tz2IX99fdGpBrMbJ5zLqOya3VyxE9EQsvUxVv58eRviY+L4p0bhuiXt1SLmXHBwHa8d+NQEmKjuPj5b3niiyx8mvoVOWYq/ESkxjjnePLLNdz46nx6t2nCuzcMJTVZnbsSmO6tGvP+zcM4u09rHv50NVe8+B0792nqV+RYqPATkRpRUubjl28v5i+frOKcvq155epBNGsU43UsqafiY6N49MJ+PDCxF9/m7GTsY1/z3bpdXscSqXdU+IlI0OUfKOHyF+byRuYmbhnVhUcvVOeuHD8z4+JBHXjnp0OIjY7gwsnf8sz0bE39igRAhZ+IBNWhzt3v1u3i/yb15eend6XyjXREjk2vNk344OZhjOmZwoP/WcnVL2Wye/9Br2OJ1Asq/EQkaOat382EJ2eRW1DMy1cN4rwBbb2OJCGqcVw0T150Ar8/pyczs3IZ+9hM5m/Y7XUskTpPhZ+IBMWHi7fw4+f+27l7Umd17krNMjMuH9KRt386hIgI44JnvuH5mTlomTKRI1PhJyLH5VDn7k2vLqBvW3XuSu3r0zaRqTcP55RuLbh/6gque3ke+YUlXscSqZNU+InIMTtY6uPOt8o7d8f3U+eueKdJw2gmXzqAe8Z254uVOxj7+EwWb9rjdSyROkeFn4gck/zC8j1335y3iVtGp/G3H/UjNkqdu+IdM+Pq4Z154/rB+HyO85/+hn/OXqepX5EKVPiJSMA27Czk3Kdn/bdz97R0de5KnXFC+6ZMvWU4w9KS+O37y7jp1QXsLdLUrwio8BORAM1bv5uJT80ib99Bde5KndW0UQzPX5bBXWd24+Nl2zjn8a9ZtiXf61ginlPhJyLV9sGi8s7dhLgo3lXnrtRxERHG9SNT+fe1J3GgpIyJT83mX3PWa+pXwpoKPxE5qkOduze/Vt65+84NQ+mszl2pJwZ2bMZHtwxnUKdm/Prdpdz2+kL2F5d6HUvEEyr8RKRKB0t93OHv3J2gzl2pp5rHx/LPn5zI7ael88GiLYx74mtWbtvrdSyRWqfCT0SOKL+wfM/dt+Zt4tbRaTyizl2pxyIijJtHp/HK1YMoKCplwpOzeCNzo9exRGqVCj8RqdShzt3M9bv46wV9+Zk6dyVEDElNYuotw+jfril3vrWY299YROFBTf1KeFDhJyI/MG/9LiY8NYud+w/yylWDOPcEde5KaGmREMcrVw/iltFpvLNgE+OfmEXW9gKvY4nUOBV+IvI/3l+0hR8/N4fGcVG8e8NQBqlzV0JUZITx89PSeenKE9m1/yDnPDGLd+Zv8jqWSI1S4SciQHnn7hNfZHHLawvo1zaRd28YSqekRl7HEqlxw9OS+ejW4fRu24Sfv7GIX761mKKSMq9jidQIFX4iwsFSH794czEPf7qaif3b8PLVJ9JUnbsSRlIax/Hq1YO44eRUXs/cyMSnZrOn8KDXsUSCToWfSJjLLyzhshfm8Pb8Tdx2ahp/vaCvOnclLEVFRnDnGd144YoMsnfs48ZX51NS5vM6lkhQqfATCWPrd+5n4tOzmL9+D4/8qC+3narOXZFR3VJ4YGIvZq3Zyf0fLvc6jkhQRXkdQES8kbluF9e+PA+fc7xy9SBO7NTM60gidcakjHas3l7AczPXkt4ygYsHdfA6kkhQqPATCUNTFm7mjrcW0yaxAS9cMVBNHCKVuOvM7mTt2Mdvpyyjc1I8g1PV4S71n6Z6RcKIc47HP8/i1n8vpF+7RN756RAVfSJHEBlhPPbj/nRo3pAb/jWPDTsLvY4kctxU+ImEiTKf4463FvN/n/k7d69S567I0TSOi+b5ywfic3D1S99RUFTidSSR46LCTyRMPD8zh7fmbeKW0ercFQlEp6RGPHXxCWTn7udnry+kzOe8jiRyzFT4iYSBrO0F/N9nqxnTM4WfnZqmzl2RAA3tksRvzu7BtBU7ePjTVV7HETlmau4QCXGlZT5+8eYi4mOjeGBibxV9IsfossEdWLW9gKe/yiY9JZ6J/bWHtdQ/GvETCXHPzshh0aZ8/jC+F0nxsV7HEam3zIzfn9OTQZ2a8cu3l7Bgw26vI4kETIWfSAhbuW0vf5u2mrF9WjG2Tyuv44jUe9GRETx9yQBSGsdy7cvz2Jp/wOtIIgFR4ScSokrKfNz+xiKaNIjmD+N7eR1HJGQ0axTD85cNpLC4lGtfmseBg2VeRxKpNhV+IiHqqS+zWbZlL/dP6E0zLdsiElRdWybw6IX9WbolnzvfXoxz6vSV+kGFn0gIWrYln8e/yGJ8v9ac0aul13FEQtKpPVK4Y0xXPli0hSe/XON1HJFqUVevSIg5WFo+xdu0UQy/G9fT6zgiIe2nI1NZva2Ahz9dTZcWCfqHltR5GvETCTFPfJHFym0F/HFib+3MIVLDzIwHz+tD33aJ/PyNhazYutfrSCJVUuEnEkKWbMrnya+yOfeENpzWI8XrOCJhIS46kucuHUBCXBRX/zOTvH3FXkcSOSIVfiIhori0jNvfXEhSfAy/PVtTvCK1qUXjOJ67LIO8fcX89JV5HCz1eR1JpFIq/ERCxKPTsli9fR8PnteHJg2jvY4jEnb6tE3k4Ul9+W7dbu59b6k6faVOUnOHSAhYsGE3z0zP5kcZ7Tilawuv44iErXF9W7N6ewGPf7GGri0TuHJYJ68jifwPjfiJ1HNFJWX84s1FtGwcx6/P7u51HJGw97NT0xnTM4X7py5n+upcr+OI/A8VfiL13F8/W0127n4eOr8PjeM0xSvitYgI468X9CM9JYGbXp1Pdu4+ryOJfE+Fn0g9Nm/9Lp6bmcNFg9ozPC3Z6zgi4tcoNornL88gJjKCa/6ZSX5hideRRAAVfiL11oGDZfzizcW0btKAu8/SFK9IXdO2aUOeuXQAG3cXctNr8yktU6eveE+Fn0g99ZdPVrE2bz9/Ob8P8bHq0xKpiwZ2bMb9E3oxMyuPBz5a4XUcEXX1itRHc3J28o/Za7lscAeGdEnyOo6IVOFHA9uzats+Xpi1lq4pCVx4YnuvI0kYO2LhZ2Y+IOBFiJxzkceVSESqVHiwlDveWky7pg355RndvI4jItVw91ndWJO7j3unLKVzcjwndmrmdSQJU1WN+N3HDwu/CUAv4BNgFWBAV+B0YAkwpQYyikgFD/1nJRt2FfL6tSfRSFO8IvVCVGQEj/+4PxOfnMX1r8xjyo1DadesodexJAwd8beGc+53FV+b2ZVAS6CXc27VYde6A18CG2ogo4j4zV6Txz+/Wc9PhnZkUOfmXscRkQA0aRDN85dnMOHJWVzzUiZv/XSIns+VWhdIc8edwBOHF30AzrkVwJPAL4MVTET+177i8ineTkmNuHOMpnhF6qPOyfE8cdEJrN5ewM9eX4jPp23dpHYFUvh1AA5Ucb3Q/x4RqQF//GgFW/IP8PCkPjSI0aO0IvXViPRk7hnbg8+Wb+evn632Oo6EmUAKv9XAtWaWePgFM2sKXEv5c38iEmQzVufy6pwNXDO8MwM66KFwkfruJ0M7cuHAdjzx5RqmLNzsdRwJI4E8XHA38B6QZWYvU14IOqAbcAmQSHnzh4gE0d6iEu56ezGpyY34+WnpXscRkSAwM+4b34uc3P3c+dZiOjZvRN92PxhXEQm6ao/4OeemAmMob+C4DXgKeBq41X/uTP97RCSIHvhwBdv2FvHwpL7ERWuKVyRUxERF8PQlJ5AUH8u1L2eyfW+R15EkDAS0c4dz7gvn3ACgFTAYGAK0cs4NcM5Nq4mAIuHsy1U7eD1zI9eNTKV/+6ZexxGRIGseH8vzl2dQUFTKtS9lUlRS5nUkCXHHtGWbc267c26Oc+5b59z2YIcSEcgvLJ/iTU+J57ZT07yOIyI1pHurxjzyo34s2pTPL99ejHPq9JWaE1DhZ2aJZvagmS01s/1mts///R8ra/oQkWN334fLydt3kIcn9SU2SlO8IqFsTM+W/OL0dKYs3MLT07O9jiMhrNqFn5m1ARZQvp4fwFTgP5Q3eNwFzDez1kFPKBKGPlu+nbfnb+KGk1Pp01b/phIJBzee0oVxfVvzl09W8dlyTaZJzQhkxO9PQApwtnOul3PuAufcJOdcb2Cs/9qfaiKkSDjZvf8gd7+7hG4tE7h5lKZ4RcKFmfGX8/vQu00Tbvv3AlZu2+t1JAlBgRR+ZwCPOuc+OvyCc+4/wOPAmcEKJhKufvfBMnbvP8j/XdCXmKhjegxXROqpuOhIJl+aQaPYKK7+ZyY79xV7HUlCTCC/VRKAqlaZ3OR/j4gco4+XbmXKwi3cPCqNnq2beB1HRDzQskkcky/LYEdBMT/913wOlvq8jiQhJJDCbxVwvpn94B4ziwTORzt3iByznfuK+fW7S+nZujE3nJLqdRwR8VC/don85fw+zF27i9++v0ydvhI0gezc8RjwPPCFmT3Cf4u8bpQv4jwcuCq48UTCx2/eX8beohL+dc0goiM1xSsS7sb3a8OqbQU89VU23VomcPmQjl5HkhBQ7cLPOfeCmbUAfgu8U+GSAcXA3c65F4MbTyQ8fLh4C1MXb+WOMV3p1rKx13FEpI74xeldWb19H/d9uJzU5HiGpSV5HUnqOQt0+NjMmgOnAR38p9YBnznndgU3Wt2RkZHhMjMzvY4hISq3oJjTH5lOu2YNeeenQ4jSaJ+IVLCvuJTznprN1vwDTLlpGJ2SGnkdSeo4M5vnnMuo7FrAv2Gcczudc/92zj3kP14PdtFnZpPMbJmZ+cwso8L508xsnpkt8X8dVeHaV2a2yswW+o8WlXxuczP70r/w9BPBzCxyLJxz3PPeEvYXl/F/k/qq6BORH4iPjeL5yzOIjDCu+ud35B8o8TqS1GMB/5YxszPM7Akzm2pmH/q/Pz3IuZYC5wIzDjufB4zzrx14OfDyYdcvds718x87KvncIuBe4BdBzityTN5ftIVPlm3n56enk5aipngRqVy7Zg15+pIBbNhZyC2vLaDMp2YPOTaB7NwRY2ZTKN+x4wZgIDDI//1/zOw9M4sJRijn3Arn3A86hJ1zC5xzW/wvlwFxZhYbwOfud859TXkBKOKpHXuL+M2UZfRvn8g1wzt7HUdE6riTOjfnvvG9mL46lz99tMLrOFJPBTLi91tgHPB/QAvnXAvnXDKQDDwMnEP5aFptOQ9Y4JyruLrlP/zTvPeamdViFpGAOOe4+90lFJWU8fCkvkRG6K+riBzdRYPac8WQjjz/9Vre+G6j13GkHgqk8LsIeMU5d6dzLu/QSf8zf78EXgEuqe6Hmdk0M1tayTG+Gvf2BB4Crqtw+mL/FPBw/3FpdbMc4c+41swyzSwzNzf3eD5K5Afemb+ZaSt2cMeYrqQmx3sdR0TqkXvGdmdYlyR+/d4SMteFbF+l1JBACr/WwOwqrn8DtKruhznnTvXv+Xv4MaWq+8ysLfAucJlzLrvC5232fy0AXgVOrG6WI+Sb7JzLcM5lJCcnH89HifyPbflF/O6DZQzs2JSfDO3kdRwRqWeiIiN44qL+tElswPWvzGPT7kKvI0k9EkjhtwU4qYrrg4CtxxenamaWSPkzhr9yzs2qcD7KzJL830cDZ1PeICJSpzjnuOudxZSU+fjL+ZriFZFjk9gwhucvH0hxqY9rXprH/uJSryNJPRFI4fcacKmZ3W9mTQ+dNLOmZvYHyqdWXw1GKDObaGabgMHAVDP7xH/pJqALcO9hy7bEAp+Y2WJgIeV7Cj/n/6xzzOy+Cp+9DvgrcIWZbTKzHsHILFIdb2Zu4qtVudx1Rjc6ai0uETkOXVrE8/iP+7Nq215uf2MRPnX6SjVUewFnf/fsu8AZgAMOPfiWTPnuHR8DEw9rtggJWsBZgmHzngOc8cgMerRuzGvXnESERvtEJAien5nD/VNXcMvoNH5+WrrXcaQOqGoB50C2bCsGzjKzs4GxQEf/pXXAB865j44zp0jIcs5x19uLKXOOv5zfV0WfiATNVcM6sWpbAY99nkV6Sjxn92ntdSSpw6pd+B3inPsQ+LAGsoiErNfmbmRmVh5/mNCL9s0beh1HREKImXH/xF7k5O3nF28uomPzRvRq08TrWFJHaX8okRq2cVchD0xdztAuzbn4xPZexxGREBQbFckzlwygWcMYrnkpkx17tU+BVC6gws+/V+7rZvadmWWbWc5hR/bRP0UkfPh8jjvfWgzAQ+f10RSviNSY5IRYnrs8gz2FJVz78jyKSsq8jiR1UCBbtv2M8gaOkZR3zc4Aph92HL63rkhYe2XOer7J2ck9Z/egbVNN8YpIzerZugmP/KgvCzfu4e53llDdBk4JH4E84/czyou7M5xzB2soj0jIWL9zP3/6aCXD05K4cGA7r+OISJg4o1crfn5aOn/9bDVdWyZw3chUryNJHRLIVG8S8LqKPpGj8/kcd7y5mKgI46Hz+qCto0WkNt08qgtj+7TiwY9XMjNL247KfwVS+M0HOtdUEJFQ8uLsdcxdt4vfjOtB68QGXscRkTBjZjx8fl86JTXi1+8u1fN+8r1ACr/bKN+547SaCiMSCnJy9/HnT1YyqlsLzh/Q1us4IhKmGsREcv/4XmzYVchTX67xOo7UEUd8xs/MPq3kdAHwsX/bs3XA4f+EcM65MUFLJ1LPlPkcd7y1mJjICP50bm9N8YqIp4Z0SWJCv9Y8Mz2HCf3b0Dk53utI4rGqRvzSgbTDjhhgg/++zpVc114xEtZe+Hot89bv5vfje5LSOM7rOCIi3D22O7HREdw7Zam6fOXII37OuY61mEOk3luzo4C/fLqK03qkMKFfG6/jiIgA0CIhjjvGdOU3U5bx/qItjNf/n8Kadu4QCYLSMh+3v7mYhjGRPDCxl6Z4RaROuXhQB/q0bcL9U1ewt6jE6zjiIRV+IkEweWYOizbu4b7xvWiRoCleEalbIiOM+yf0Im9fMX/9dLXXccRDRyz8zMxnZqVmFlPhddlRjtLaiy5SN6zaVsDfPsvizF4tGdenlddxREQq1adtIpee1IGXvlnH0s35XscRj1S1c8d9gANKD3stIn4lZT5+8eYi4uOi+MMETfGKSN12++ld+WjJNn797hLeuWEokdo/POxU1dzxu6peiwg881U2Szbn89TFJ5AUH+t1HBGRKjVpEM09Y7tz2+sLeXXuBi49qYPXkaSW6Rk/kWO0fMteHvsii3F9W3NWb03xikj9ML5fa4akNufPH68kt6DY6zhSy6pawHnEsQq+RrIAACAASURBVHygc27GsccRqR8OlpZP8TZpEMN95/T0Oo6ISLWZGfeN78WZj87gTx+t4K8/6ud1JKlFVT3j9xWBPdNn/vdHHk8gkfrgyS/XsHzrXiZfOoCmjWK8jiMiEpAuLeK5bkQqT3y5hkkZ7Ric2tzrSFJLqir8Tqm1FCL1yNLN+Tz55Rom9m/D6T1beh1HROSY3HhKF95buJl7pyzlo1uGExOlp7/CQVXNHdNrM4hIfVBcWsbtbyyiWaMYfjuuh9dxRESOWYOYSO4b35MrX8zkuZk53HhKF68jSS04pvLezNLMbKiZNQl2IJG67LHPs1i1vYA/ndubxIaa4hWR+m1UtxTG9Ezh8S+y2Lir0Os4UgsCKvzM7Edmth5YCcwABvjPJ5lZlplNqoGMInXCoo17ePqrbM4f0JbR3VO8jiMiEhS/HdeTCDN+/8Eyr6NILah24Wdm44HXgA3AvZQ3cwDgnMsDVgCXBjugSF1QVFLG7W8uokVCHPeerSleEQkdrRMbcNupaUxbsYNPl23zOo7UsEBG/O4BZjjnhgPPVnJ9DtA3KKlE6phHpq1mzY59PHheb5o0iPY6johIUP1kaCe6piTw+w+WU3hQu6+GskAKv57AG1Vc3wZo/ktCzrz1u3luRg4XDmzHyV1beB1HRCTooiMjuH9iLzbvOcCjn2d5HUdqUCCFXxEQV8X1jsCe40ojUscUlZRxx5uLaNWkAb8e293rOCIiNWZgx2ZMGtCWv89cy6ptBV7HkRoSSOH3NfDjyi74u3uvBL4IRiiRuuK1uRvIydvPn87tTUKcpnhFJLT96qzuxMdFce97S3EukD0cpL4IpPD7HdDTzL4EzvWfyzCzm4CFQGPgD8GNJ+KdkjIfz89cS0aHpoxIT/Y6johIjWvWKIa7zujG3HW7eHv+Zq/jSA2oduHnnJsPjAFa8t/mjgeBx4CDwBjn3IqgJxTxyEdLtrJ5zwGuG5nqdRQRkVpzQUY7TmifyB8/WsGewoNex5EgC2gdP+fcTOdcd6A/8CPKp34HAt2cc7PNLL4GMorUOuccz0zPITW5EaO7qaFDRMJHRIRx/4Te5B8o4aGPV3kdR4IskHX8njj0vXNukXPuTefc6865ec45Z2aJwLQaSSlSy2Zm5bFi616uG5FKRIQd/QYRkRDSo3VjrhjSkdfmbmD+ht1ex5EgCmTE7zoz+2NlF8ysGeWNHV2DkkrEY5Nn5NAiIZbx/Vt7HUVExBM/Oy2dlo3juOfdpZSW+byOI0ESSOF3FfBLM/tVxZNm1oLy7ds6AKcHMZuIJ5ZuzufrNXlcOawTsVGRXscREfFEfGwUvxnXg+Vb9/LSN+u9jiNBEkhzx0vALcADZnYjgJm1BWYCLYBRzrnvaiSlSC16dkYO8bFRXDSovddRREQ8dWavloxMT+avn61m+94ir+NIEATa3PEk8GvgUTO7C5gOJAAnO+cW1UA+kVq1cVchUxdv4aJB7WmsdftEJMyZGb8/pycHy3zc9+Fyr+NIEARU+AE45/4E/Bn4IxAFDHfO6W+DhITnZ+YQGWH8ZGhHr6OIiNQJHZMacePJXZi6eCszVud6HUeOU9SRLpjZ5KPcWwCspPy5v0PnnHPuuiBlE6lVu/Yf5PXMjYzv14ZWTRp4HUdEpM64/uTOvLdwM7+ZspSPbxtBXLSef66vjlj4AVdX4/7TDnvtABV+Ui+9/M16ikp8XDuis9dRRETqlNioSP4wvheX/H0Oz0zP5rZT072OJMfoiFO9zrmIYzj0TwCplw4cLOOf36xjdLcWpKckeB1HRKTOGZaWxLi+rXnqq2zW5e33Oo4co4Cf8RMJRW/N28iu/Qc12iciUoV7xnYnJjKCe6csxTnndRw5Bir8JOyVlvl4buZa+rVL5MROzbyOIyJSZ6U0juP209OZmZXH1CVbvY4jx6Cq5o61gI/yfXhL/K+PVt4755x2tJd65eNl29iwq5C7z+pGhUYlERGpxKUndeCteZu474PljExPJkFLX9UrVTV3TKe80PMd9lokZDjneHZ6Dp2SGnFaj5ZexxERqfOiIiN4YGJvJj41i0c+y+I343p4HUkCcMTCzzl3RVWvRULBNzk7WbI5nz9O7E1khEb7RESqo1+7RC46sT0vzl7LeQPa0LN1E68jSTXpGT8Ja89OzyEpPoZzT2jjdRQRkXrlzjHdaNowhnveW4rPpwnB+qKqZ/yOaaNS59yGY48jUntWbN3L9NW5/OL0dC1GKiISoCYNo7n7rO7c/uYi/v3dRu1vXk9U9YzfOo7tmT79BpV6YfKMHBrGRHLJSR28jiIiUi+de0Ib3sjcyEMfr2RMzxSax8d6HUmOoqrC70rUzCEhavOeA7y/aAuXD+5IYsMYr+OIiNRLZsb9E3px5qMz+dN/VvLwpL5eR5KjqKq548VazCFSq174ei0AVw3v5HESEZH6LS0lgWtGdObpr7K5IKOd1kOt49TcIWEnv7CE1+Zu4Jy+rWmT2MDrOCIi9d7No7rQJrEB97y3hJIy39FvEM+o8JOw88qc9RQeLOOa4dqeTUQkGBrGRPG7c3qyevs+/u6fUZG6SYWfhJWikjL+MWstI9KT6dG6sddxRERCxmk9Uji1ewqPTsti854DXseRI1DhJ2Hlnfmbydt3kOtHaLRPRCTYfndO+S4ev39/mcdJ5EhU+EnYKPM5npuZQ+82TRic2tzrOCIiIadt04bcMjqNT5dv5/MV272OI5VQ4Sdh47Pl21mbt5/rRnbGTNuziYjUhKuGdSKtRTy/fX8ZBw6WeR1HDlPtws/M2h/laGdmyabfqFIHOed4Zno27Zs15IyeLb2OIyISsmKiIvjDhF5s2n2Ax7/I8jqOHCaQEb91wNoqjnXANmCfmX1kZicFNanIcfhu3W4WbtzD1cM7ERWpgW4RkZp0UufmnHtCG56bmcOaHQVex5EKAvkNeBWwCNgDPAncBvwMeNp/bgFwK/A8MAiYbmbDg5pW5Bg9Oz2bpg2jmTSgnddRRETCwt1ndadhTBT3vLcU57QRWF0RSOHXEmgIpDnnbnHOPe6ce8w5dxPQFUgAYp1ztwLdgB3Ab4OeWCRAWdsL+HzlDi4f0pEGMdpKWkSkNiTFx3LnGV35NmcX7y3c7HUc8Quk8LseeM45t+vwC865PMpH+m7yv84F/g4MDEZIkeMxeUYOcdERXDa4o9dRRETCyo8Htqdvu0QemLqC/MISr+MIgRV+LYDoKq5HUT4qeMgmqtgLWKQ2bMsv4r2Fm/lRRjuaNYrxOo6ISFiJiDAemNCLXfsP8pdPV3odRwis8FsM3Ghm7Q+/YGYdgBspfwbwkDRg6/HFEzk+/5i1ljKf42ptzyYi4olebZpw2eCO/GvOBhZt3ON1nLAXSOF3O9AMWGVmb5nZn/zHW8BKoDnwCwAziwMuBj4PdmCR6tpbVMK/5mzgrN6taNesoddxRETC1u2np5McH8s97y2lzKdGDy9Vu/Bzzn0NDAE+Bc4Efuk/zgQ+AU7yvwfnXJFzrq1z7rrgRxapnlfnbGBfcSnXjUj1OoqISFhLiIvm3rN7sGRzPq98u97rOGEtoAXNnHOLnHPjKe/gbe0/EpxzE5xzi6q+u/rMbJKZLTMzn5llVDh/mpnNM7Ml/q+jKlz7ysxWmdlC/9Giks894v0SWopLy3jh67UM7dKc3m2beB1HRCTsnd2nFcPTknj4k1Xs2FvkdZywdUwr2TrnfEAZUOr/PtiWAucCMw47nweMc871Bi4HXj7s+sXOuX7+Y0cln3u0+yVETFm4hR0FxRrtExGpI8yM35/Tk+JSH/dPXeF1nLAVUOFnZp3N7N9mlk/5Lh3bzSzfzF41s6A9Pe+cW+GcW1XJ+QXOuS3+l8uAODOLDeBzj+t+qR98PsfkGTl0b9WY4WlJXscRERG/zsnxXH9yKu8v2sKsNXlexwlLgezV2w34DjgfmAn8H/BX//eTgLn+99SW84AFzrniCuf+4Z/mvbcaewZXdr+EgC9W7mDNjn1cP7Iz2jpaRKRuueHkVDo0b8i97y2luLTM6zhhJ5ARvwcBH9DfOXe2c+5O59wdzrmzgf6AA/5Y3Q8zs2lmtrSSY3w17u0JPARUbB652D+FO9x/XBrg/Ye/51ozyzSzzNzc3Or+WFIHPDsjmzaJDTirdyuvo4iIyGHioiO5b3wvcvL2M3l6jtdxwk4ghd9I4HHn3JLDLzjnlgJPAKdU98Occ6c653pVckyp6j4zawu8C1zmnMuu8Hmb/V8LgFeBEwO5v5J8k51zGc65jOTk5Or+WOKxeet389263Vw1rBPRkcf0CKuIiNSwkenJjO3diie+XMOGnYVexwkrgfxmjAH2VnE93/+eGmNmicBU4FfOuVkVzkeZWZL/+2jgbMobRKp1v4SOyTOyadIgmh8NbOd1FBERqcK9Z/cgKsL4zftLcU5r+9WWQHfuuNzMGhx+wX/ucv97jpuZTTSzTcBgYKqZfeK/dBPQBbj3sGVbYoFPzGwxsBDYDDzn/6xzzOy+o9wvISA7dx+fLt/OZYM70ChWuwWKiNRlLZvE8bPT0vlqVS6fLNvmdZywYdWtss1sHPAekAU8DRzquu0GXE95QTXBOfdhDeT0VEZGhsvMzPQ6hhzFr95ZzNvzNzP7rlEkxatZW0Skrist8zHuiVnsKTzItJ+P1D/ag8TM5jnnMiq7FsjOHR8AlwCNgUconzKdSnlnb2PgklAs+qR+2FFQxNvzNnP+gLYq+kRE6omoyAjun9CLrflF/G3aaq/jhIWASmvn3Gtm9iYwAOjoP70OyHTOqSdbPPPirHWU+HxcMzxoy0mKiEgtGNChKT8+sR0vzFrHuSe0pXurxl5HCmkBtz0650qdc3Occ6/7jzkq+sRL+4pLefnb9ZzRsyWdkhp5HUdERAJ055huNGkQzT3vLcXnU6NHTTriiJ+ZtT+WD3TObTj2OCKB+/fcDRQUlXLtCI32iYjUR00bxXDXmd24863FvDVvExdoZYYaU9VU7zrKF2UOVOSxRREJXEmZj79/vZZBnZrRv31Tr+OIiMgxOv+EtryZuZE//WcFp/VIoWmjGl0hLmxVVfhdybEVfiK15oNFW9iaX8QfJ/b2OoqIiByHiAjj/gm9GfvYTB78z0oeOr+P15FC0hELP+fci7WYQyRgzjmenZ5Deko8J3fV7ioiIvVd15YJXDWsE8/OyGFSRlsyOjbzOlLI0Z5WUm99tTqXVdsLuHZEKmbmdRwREQmCW0an0bpJHPe8t5TSMp/XcUKOCj+pt56dnk3LxnGc07e111FERCRIGsVG8ZtxPVm5rYAXZ6/zOk7IUeEn9dKijXv4NmcXVw3rREyU/hqLiISSMT1TGNWtBY98tpqt+Qe8jhNS9BtT6qXJM3JIiIviwhPV8i8iEmrMjN+f05My57jvg+VexwkpKvyk3lmXt5//LN3KJSd1ICEu2us4IiJSA9o1a8jNo9L4z9JtfLlqh9dxQoYKP6l3nv86h6iICH4ypKPXUUREpAZdPbwTnZMb8dspyygq0SZhwRBQ4WdmMWZ2lZn9y8w+M7P+/vNNzewyM2tbMzFFyuXtK+bNzE1M7N+GFo3jvI4jIiI1KDYqkvvH92LDrkKemZ7tdZyQUO3Cz8yaAt8CzwFnAaOAQ1sl5AN/AG4KdkCRil76Zj3FpT6u0fZsIiJhYUiXJM7s1ZK/z1xLfmGJ13HqvUBG/B4E0oHTgTTg+4XTnHM+4B3gjKCmE6mg8GApL32zjtN6pNClRbzXcUREpJbcMjqNguJSXpi11uso9V4ghd85wGPOuWlUvpXbGqBDUFKJVOKN7zayp7CE60dqtE9EJJx0b9WYMT1TeGHWWvIPaNTveARS+DUFcqq4HgloR2WpEaVlPp6buZaMDk0Z0EFb+IiIhJtbRqdRUFTKP7Wo83EJpPBbC/Su4vpIYNXxxRGp3NQlW9m85wDX6tk+EZGw1LN1E07tnsLfv15LQZFG/Y5VIIXfK8DVZjaqwjkHYGY/AyYALwYvmkg55xzPTs+hc3IjTu2e4nUcERHxyK2j08g/UKJRv+MQSOH3EPA58Bnl3b0OeMLMtgP/B0wBHg96Qgl7s9bsZPnWvVw3ojMREXb0G0REJCT1btuEUd1a8PzXa9lXXOp1nHqp2oWfc67UOXcOcAmwFFjpv38ucKlz7lznXGVNHyLH5dkZ2SQnxDKhfxuvo4iIiMduHZ3GnsISXvpmnddR6qWoQG9wzr0GvFYDWUR+YOnmfGZm5fHLM7oRGxXpdRwREfFY33aJnNw1mednruXywR1pFBtwKRPWtGWb1GmTZ+TQKCaSiwa19zqKiIjUEbeMTmPX/oO88u16r6PUO0csk83shWP4POecu+o48oh8b+OuQqYu2cqVQzvSpEG013FERKSOOKF9U4anJTF5Rg6XDu5AwxiN+lVXVf+lRlH5Qs1V0TN+EjR//3otBlw5rJPXUUREpI657dQ0znv6G16ds4Grh2upr+o6YuHnnOtYizlE/sfu/Qd5/buNjO/XhlZNGngdR0RE6pgBHZoxtEtznpmew8WDOtAgRs+BV4ee8ZM66eVv13OgpEwLNouIyBHdOjqdvH3FvDZ3g9dR6g0VflLnFJWU8eLsdYzq1oKuLRO8jiMiInXUiZ2acVLnZjwzPZuikjKv49QLRyz8zGytmWWbWXSF1zlHObJrL7qEqjfnbWLX/oMa7RMRkaO6dXQ6OwqK+bdG/aqlquaO6ZQ3a/gOey1SY8p8judn5tC3XSKDOjXzOo6IiNRxg1Obc2KnZjw9PZsLT2xPXLSe9atKVYXfLUChc64MwDl3Ra0kkrD2ybJtrN9ZyF1ndMNM27OJiMjR3To6jYufn8ObmRu5dHBHr+PUaVU947cbmHTohZm9YGaDaj6ShCvnHM9Oz6Zj84ac3rOl13FERKSeGJLanIwOTXnqq2yKS/WsX1WqKvxKgNgKr68AUms0jYS1b3N2sWhTPteM6ExkhEb7RESkesyMW0ansTW/iLfmbfI6Tp1W1VTvGuBiM5sH5PvPNTezKvfOcs7p6Uo5Js/OyCYpPobzTmjrdRQREalnhqcl0b99Ik99mc2kAe2IidLCJZWp6r/KfcBIYCGwlvLGjr/5v6/qEAnYym17+WpVLpcP7qgHc0VEJGBmxq2j09i85wDvzNeo35FUtXPHG/7RvpOBFOB+4C1gUe1Ek3AyeUYODaIjuXRwB6+jiIhIPTUyPZm+bZvwxJdrOG9AW6IjNep3uCp3NXbOZQPZAGZ2DfCKc+792ggm4WPLngO8v3ALlw7uQGLDGK/jiIhIPWVm3HpqGle+mMm78zdzwcB2Xkeqc6pdCjvnOqnok5rwwtdrccBVwzp5HUVEROq5U7q2oHeb8lG/0jLf0W8IMwGNgZpZlJldb2YfmtkyM1vq//5aM6ty9FCkMvmFJbw2dwPj+rSibdOGXscREZF67lCH74Zdhby3cIvXceqcahd+ZtYE+AZ4ChgC7AcOAIOBZ4DZ/veIVNsrc9az/2AZ147QSkEiIhIcp3ZvQY9WjXlSo34/EMiI3wNAf+AmIMU5d6JzbiDljR83+q/dH/yIEqqKSsr4x6x1DE9Lokfrxl7HERGREHFo1G9t3n4+WKxRv4oCKfwmAk87555yzpUcOumcK3XOPQ08C5wb7IASut5bsJm8fcVcP1KjfSIiElyn90ihW8sEHv9iDWU+53WcOiOQwq85sLyK68v87xE5Kp/PMXlGDr3aNGZIqv7aiIhIcEVElI/65eTu50ON+n0vkMJvHTCmiutn+N8jclSfrdhOTt5+rhuRipm2ZxMRkeA7o2dLuqZo1K+iQAq/F4BzzOwVM+trZnH+o5+ZvQScDTxfMzEllDjneGZ6Nu2aNeDMXi29jiMiIiEqIsK4eXQX1uzYx3+WbvU6Tp0QSOH3F2AycBEwn/Ku3v3APOASYLJz7uGgJ5SQk7l+Nws27OHqYZ2J0qrqIiJSg87s1YouLeJ57PMsfBr1C2gBZ+ecux7oDdxNeRE42f99b+fcT2smooSaZ6dn07RhNJMy2nodRUREQlxkhHHzqC6s3r6Pj5dt8zqO56q16LKZNQRmAs85556hvJFDJGBrdhQwbcUObh2dRsMYrfktIiI17+w+rXn08ywe+zyLM3q2JCIifJ8tr9aIn3OuEOgMaBVEOS6TZ+QQFx3BZYM7eB1FRETCxKFRv5XbCvh0+Xav43gqkAesvgROrqEcEga27y3i3QWbuSCjHc3jY72OIyIiYWRcn9Z0SmrEY59n4Vz4PusXSOF3C9DXzB4xs3TtzSuBemHWWsp8jquHdfY6ioiIhJmoyAhuPKULy7fuZdqKHV7H8Uwghd9aIJ3yAnAFUGxmBw87imskpdR7BUUlvPrtBs7s3Yr2zRt6HUdERMLQhH6t6dC8IY9+vjpsR/0CGbX7FxCe/5XkuL02dwMFxaVcN0KjfSIi4o1Do353vrWYL1ftYFS3FK8j1bpqF37OuStqMIeEsIOlPv7+9VqGpDanT9tEr+OIiEgYm9i/DY9/kcWj07I4pWuLsNs9SqvnSo2bsnAz2/cWc93IVK+jiIhImIuOjODGk7uwaFM+X63O9TpOrQuo8DOzRDN7wMwWmVm+/1jkP9e0pkJK/eXzOSbPyKFbywRGpCV5HUdERIRzT2hLm8QGPDot/Dp8q134mVkqsBj4FRAJTAM+93//K2Cx/z0i3/ty1Q6yduzj+pGpYTecLiIidVNMVAQ3nJLKwo17mJmV53WcWhXIiN9jQFPgdOdcL+fcec65c51zvYAxQKL/PSLfe3Z6Dq2bxDG2Tyuvo4iIiHzv/AFtad0kjkfDbF2/QAq/kcDfnHPTDr/gnPuM8qJvZLCCSf03f8Nu5q7bxVXDOxMdqcdJRUSk7oiNiuSnJ6cyb/1uZmfv9DpOrQnkt/E+oKqnILf73yMCwOTpOTRpEM2FA9t5HUVEROQHLhjYjpaN48LqWb9ACr/XgEvMLObwC2YWB1wKvBqsYFK/5eTu45Pl27j0pA40itUmLyIiUvfERkVy/cjOzF23i29zdnkdp1YEUvh9AEQDC83sFjM7w8zGmNmtwDzKmzw+MLMhFY+aCC1133Mz1xIdGcHlQzp6HUVEROSILjyxPS0SYnn089VeR6kVgQzFVHy272/8dxcPO8J7zP+eyGOLJvVVbkExb8/fxPkD2pKcEOt1HBERkSOKi47k+pGp3Pfhcubk7GRQ5+ZeR6pRgRR+P6mxFBJSXvl2PSVlPq4Zru3ZRESk7rtoUHue+iqbx77I4l8q/Mo55/5Zk0EkNDjnmLJwM0NSm9MpqZHXcURERI6qfNSvM/dPXUHmul1kdGzmdaQaUyfX2DCzSWa2zMx8ZpZR4fxpZjbPzJb4v46qcO0rM1tlZgv9R4tKPvfECtcXmdnE2vqZwsXSzXtZt7OQc/q29jqKiIhItV00qD3NG8Xw6OdZXkepUXWy8AOWAucCMw47nweMc871Bi4HXj7s+sXOuX7+Y8cRPjfDOdcPOAN41szUchpEHyzeQnSkMaZnS6+jiIiIVFvDmCiuHdGZmVl5zN+w2+s4NaZOFn7OuRXOuVWVnF/gnNvif7kMiDOzancPOOcKnXOl/pdx/LdBRYLA53N8uGgLw9OSSWz4g1V/RERE6rRLTupAs0YxPDotdEf96mThV03nAQucc8UVzv3DP417rx1hY1gzG2Rmy4AlwPUVCkE5Tgs27mZLfhHj+mp7NhERqX8axUZx9fBOTF+dy8KNe7yOUyM8K/zMbJqZLa3kGF+Ne3sCDwHXVTh9sX8KeLj/uLSye51zc5xzPYGBwK/8i09X9mdca2aZZpaZm1vVhiVyyAeLthIbFcGp3VO8jiIiInJMLhvckcSG0TwWos/6eVb4OedOdc71quSYUtV9ZtYWeBe4zDmXXeHzNvu/FlC+g8iJR/nzVwD7gV5HuD7ZOZfhnMtITk4O7IcLQ2U+x4eLtzKqWwsS4qK9jiMiInJM4mOjuHpYJ75YuYMlm/K9jhN01S78zKynmZ172LlTzOxzf4ft7cGP94MMicBU4FfOuVkVzkeZWZL/+2jgbMobOQ6/v9OhZg4z6wB0BdbVdO5wMCdnJ3n7ijm7j7p5RUSkfrt8SEeaNIgOyQ7fQEb8/gxcfeiFmbUB3gf6AA2AP5vZZcEIZWYTzWwTMBiYamaf+C/dBHQB7j1s2ZZY4BMzWwwsBDYDz/k/6xwzu89//zBgkZktpHzU8AbnXF4wMoe7DxZvpWFMJKO6/WAVHRERkXolIS6aq4Z1YtqK7SzdHFqjfuZc9RpbzWwr8Dfn3EP+13cCvwPSnHObzexDINk5N6imwnolIyPDZWZmeh2jziop8zHwgWmMTE/m0Qv7ex1HRETkuOUfKGHYQ18wJLU5z16acfQb6hAzm+ecqzR0ICN+TYHtFV6fCUw/9Gwd8AGQfmwRpT77ek0eewpLGKdpXhERCRFNGkRz5dBOfLJsOyu27vU6TtAEUvjlAW0BzKwR5dOw0ypcjyKwvX8lRHywaAsJcVEMT0/yOoqIiEjQXDm0EwmxUTz+Reg86xdI4TcT+Km/weNvQDTlz/gd0pXyZ+skjBSVlPHZsu2c0bMlsVGRXscREREJmiYNo7liaEc+WrKNVdsKvI4TFIEUfncDB/6/vTuPk6Ou8z/+emdyh4RAbgIhhHAmJAIB5VIukSPR9VhX5Vh0FVEOUXbdZfEA97fuevxcXdcDRFHwWBU5nIRDkfsQCMJMEiBASAJkJvdJCLnms39UDbRNz9k9U328n49HP7q76lvVn291p+aT71EF3AD8A/D1iHgOQFId8AHg3pJHaGXt3mdXsWnrDmb53rxmZlaF/uHYfRjSv47/rpJWv053zUbEYkkHAgcDGyJiac7qwcCngIYSx2dlrr6hid2H9OfofUdkHYqZOlhIZQAAIABJREFUmVnJDR/cn78/eiI/uHcRz63YxH5jhmYdUlG6dAHniNgREY15SR8RsSkibomIJSWNzsraq9t28KenV3La1LH0ravku/+ZmZm17ePHTWJQvzq+e9fzWYdStC79tZa0u6R/k/SgpOckHZUuHyHpS2mLoNWIPz29ki3bd7qb18zMqtruQ/pzzlETqW9s4vmVr2QdTlG6cueOvUgujvx5YCgwieTCzUTEGuDDwAU9EKOVqfqGJsYMG8ARE3fPOhQzM7Me9Ynj9mFg3zq+d3dlt/p19c4dA4G3ACcCylt/S7rcasDG17Zzz8JVnH7IOOr65P8UzMzMqsuIXQZw9lF7c8uTy3hhVeW2+nUl8Xsn8N8R8TRQ6HYfi4G9ShKVlb0/LljBtp0t7uY1M7Oa8YnjJtG/bx++d/eirEPptq4kfkOAle2s36XIWKyC1Dc2MX74IA7da3jWoZiZmfWKUUMHcOZb9+bmJ5exdM3mrMPplq4kfguBt7Wz/nRgfnHhWCVYu3kbDzy3mlnT90ByN6+ZmdWOT759En37iP+p0Bm+XUn8rgLOkvQxoPUWDSFpqKT/Ao4Hvl/i+KwM3T5/OTtagpnTxmUdipmZWa8aPWwgHz5yAjc+sYyX1r6adThd1unELyJ+AFwNXAO8kC6+AVgHfAb4bkT8vOQRWtmZ3djEpJFDmLLHsKxDMTMz63WfOn5f6vqoImf4dvUCzhcCx5Akf7cBjwI/AI6LiEtKH56Vm5UbX+PhF9Yw0928ZmZWo8YMG8iHjtiLGx5/mZfXVVarX5dvtxARD0fEJRFxRkScFhEXRcSDPRGclZ9b5zUTAbPczWtmZjXsU8fvSx+J799TWTN8fZ8t65L6xmYOHDu04u9VaGZmVoxxuw7ig0fsyW/nvsSy9VuyDqfTunrLto9LekjScklbJW3Le2ztqUAte8vWb+Hxpet87T4zMzPgU8dPBuCHFdTq17ezBSV9E/gs0AQ8DKzvqaCsPM1pbALwbF4zMzNg/PBBfODwvfj1Yy/x6RP2Zdyug7IOqUOdTvyAj5FM6HhPROzsoXisjNU3NDN9z13Ze8SQrEMxMzMrC58+fl9+O/clrrr3Ba5495Ssw+lQV7p6BdQ76atNS1ZvZt6yDcyc5m5eMzOzVnvtPpj3H7Ynv3z0RVZsfC3rcDrUlcTvDuCIngrEytvstJv3DHfzmpmZ/ZULTpjMzpbgh/eW/1i/riR+FwNHSPqypPE9FZCVp/qGZo6YuBt7DC//8QtmZma9acKIwbz30PH88pEXWbmpvFv9unLnjpXA9cCXgBclbfes3tqwcPkmFq7Y5Nm8ZmZmbbjwhMls39nC1fe+0HHhDHVlVu9XgMtJZvXOxbN6a8bsxib6CE6b6m5eMzOzQiaOHMLfvGU8P39kKZ98x76MGjog65AK6sqs3k/iWb01JyKY3djMUfuOKNsfsZmZWTm48MTJ3PzkMq65/wUuO/2grMMpqCtj/AbiWb01Z0HTRhav3swsz+Y1MzNr16RRu/Du6Xtw3cNLWfNKeY5+60ri9wc8q7fm1Dc00bePOHXq2KxDMTMzK3sXnjiZ13bs5JoHFmcdSkFdSfwuBA6XdKVn9daG1m7e4/YbyfDB/bMOx8zMrOxNHj2UmdP24LqHlrBu87asw3mTriR+LwNTgS/gWb014S8vrmfZ+i2ezWtmZtYFF504mVe37+SaB8pvhm9XJnf8AoieCsTKT31DE/379uGdB4/JOhQzM7OKsf+YoZw+dRw/e2gpnzhuUln1mnU68YuIc3swDiszO1uCOfOaOeGAUQwd2C/rcMzMzCrKRSdNZs68Zn7ywGI+d8oBWYfzuq509VoNeWTxGlZt2upuXjMzs244cOwwTp0ylmsfXMKGLduzDud1bbb4SXo7QETcl/u+I63lrbLNbmxmcP86TjxwdNahmJmZVaSLT9qP2xcs59oHF3PJyftnHQ7QflfvPUBIGhQR21rft1Ne6fq6kkVnmdi+s4Xb5jVz8kFjGNy/K8NAzczMrNXBewzjlIPH8JMHFvOxY/dhWBkMnWrvr/oJAGnS9/p7q34PPr+ada9uZ+Y036LNzMysGBeftB9/eGoFP3twCRedtF/W4bSd+EXEve29t+o1u7GZoQP78o4DRmUdipmZWUWbOn5XTj5oNNc8sJhzj5mY+YTJTk/ukHSXpJPaWX+CpLtKE5ZlZeuOndwxfznvmjKWAX3da29mZlasi0/ajw1btnPdw0uzDqVLs3qPB9q7oNto4B1FRWOZu3fhKjZt3eHZvGZmZiUybc/hnHDAKK65/wU2b92RaSxdvZxLe5M79gVeKSIWKwP1jc3sNrgfR+87IutQzMzMqsbFJ+3Hulezb/Vrd8qmpLOBs3MWXSbpowWKDgcOBf5Ywtisl726bQd3PrWC9x42nn51vsSjmZlZqRw6YTfevv8o5jdtyDSOjq7VsTvQOgUlgLHA0LwyAWwmuaXbF0oanfWqu55ZyZbtO5k1zd28ZmZmpfbDsw7L/DJp7X56RHwH+A6ApBbgkoj4ZW8EZr2vvqGJ0UMHcOQ+u2cdipmZWdXJOumDrt2r131/VWzTa9u5e+EqPnLkBOr6KOtwzMzMrAc4mTMA/vjUCrbtaPFsXjMzsyrmxM+ApJt3/PBBHDZheNahmJmZWQ9x4mes27yN+59bzcxp45DczWtmZlatnPgZty9Yzo6WcDevmZlZlXPiZ8xubGKfkUOYssewrEMxMzOzHuTEr8at3PQaDy9awyx385qZmVW9TiV+koZLOlLSvu2U2UfSOaULzXrDbfOW0xIw0928ZmZmVa/DxE/SvwMrgIeBZyU9KumwAkWPBq4tcXzWw+obmjhgzFD2H5N/QxYzMzOrNu0mfpI+CFwG3AtcAHwVmAA8JOmsng/PelLT+i3MXbqOWdPHZR2KmZmZ9YKO7txxCXB3RJzSukDSt4BfAz+VNDIivt2TAVrPmdPYDMBM35vXzMysJnTU1Xsg8LvcBRGxDjgVuA74/5L+rYdisx5W39jEIeN3ZeLIIVmHYmZmZr2goxa/lkILI6IF+JikdcDlknYDHil1cNZzlqzeTOPLG/jX0w/MOhQzMzPrJR0lfs8CxwLfL7QyIi6VtBH4MjCzxLFZD5ozL+nmPcPdvGZmZjWjo67e24D3SBrRVoGIuBL4LLBXKQOznlXf0MSMvXdj/PBBWYdiZmZmvaSjFr+fAGuB0cCatgpFxHckLQWmlzA26yHPrtjEM8s3ccWsg7MOxczMzHpRu4lfRCwDvteZHUXEzcDNpQjKetbshib6CE6f5su4mJmZ1ZJu37JN0kBJ50gaU8qArGdFBLMbm3nbpBGMHjow63DMzMysFxVzr95dSe7UMaVEsVgvWNC0kRdWb2aWb9FmZmZWc4pJ/ABUkiis19Q3NtG3jzh1ytisQzEzM7NeVmziFyWJwnpFRDC7oZlj9xvJbkP6Zx2OmZmZ9TK3+NWQJ15az7L1W5jla/eZmZnVpI4u59KeVcA+wPISxWI9rL6hif59+/DOKZ6PY2ZmVou63eIXES0RsTQitpYyIABJfytpgaQWSTNylr9T0uOS5qXPJ+asu0fSQklPpo/R7ex/gqRXJP1jqWMvVztbgjmNzRy//yiGDeyXdThmZmaWgWJa/HrSfOB9wFV5y1cDsyKiSdJU4A5gfM76MyNibif2/18kdyWpGY8tWcvKTVs9m9fMzKyGlWXiFxFPA0jKX/5EztsFwEBJA7rS6ijpb4AXgM0lCLVi1Dc0MahfHScd1GZDqJmZmVW5Yid3ZOn9wBN5Sd+1aTfvF5WfNQKShgD/DFzZW0GWg+07W7ht/nJOOmg0g/uXZa5vZmZmvSCzLEDSnUChi8ldHhG3dLDtFOBrwCk5i8+MiGWShgK/A84Grsvb9ErgvyLilQJ5Yf5nnAecBzBhwoR2y5a7hxatYe3mbe7mNTMzq3GZJX4RcXJ3tpO0J3ATcE5ELMrZ37L0eZOkXwJH8ubE763AByR9HRgOtEh6LSL+p0B8VwNXA8yYMaOir1c4u6GJoQP68o79R2UdipmZmWWoovr9JA0H5gCXRcSDOcv7AsMjYrWkfsBM4M787SPiuJxtrgBeKZT0VZOtO3Zy+4LlnDJlLAP71WUdjpmZmWWoLMf4SXqvpJeBo4A5ku5IV10ITAa+mHfZlgHAHZIagSeBZcCP0n29W9JXer8W5eG+Z1ez6bUdzJw+LutQzMzMLGOKqOhezF4xY8aMmDu3M1eJKT8X/+oJ7ntuFY9dfjL96soyzzczM7MSkvR4RMwotM6ZQBXbsm0ndz69gtOmjnPSZ2ZmZk78qtldz6zk1W07meVuXjMzM8OJX1Wrb2hi1NABvHWfEVmHYmZmZmXAiV+V2vTadu5auJIzDhlHXZ/2r1loZmZmtcGJX5W68+kVbNvR4m5eMzMze50TvypV39DM+OGDOHSv3bIOxczMzMqEE78qtP7Vbdz37CrOmDaOPu7mNTMzs5QTvyp0+/zl7GgJZk3zvXnNzMzsDU78qtDsxmYmjhjM1PHDsg7FzMzMyogTvyqzatNWHlq0mlnT90ByN6+ZmZm9wYlflbltfjMtATPdzWtmZmZ5nPhVmfqGJvYfswsHjB2adShmZmZWZpz4VZHmDVt4bMk6T+owMzOzgpz4VZE5jc0AzJzuxM/MzMzezIlfFalvaGLq+GHsM3JI1qGYmZlZGXLiVyWWrtlMw8sb3M1rZmZmbXLiVyVmp928Z0zzvXnNzMysMCd+VaK+oYnDJgxnz90GZx2KmZmZlSknflXguRWbeGb5JmZ5UoeZmZm1w4lfFahvbEaCMw5xN6+ZmZm1zYlfhYsIZjc28bZ9RjB62MCswzEzM7My5sSvwj3VvJEXVm1m5nS39pmZmVn7nPhVuPqGZur6iNOmOvEzMzOz9jnxq2Ct3bzHTh7J7kP6Zx2OmZmZlTknfhXsyZfW8/K6LZ7Na2ZmZp3ixK+C1Tc007+uD6dMGZN1KGZmZlYBnPhVqJaWYM68Jt5xwCiGDeyXdThmZmZWAZz4VajHlqxlxcat7uY1MzOzTnPiV6HqG5sY1K+Okw8anXUoZmZmViGc+FWgHTtbuHXeck48aDSD+/fNOhwzMzOrEE78KtBDi9awdvM2Zk1zN6+ZmZl1nhO/CjS7sYldBvTl+ANGZR2KmZmZVRAnfhVm646d3D5/OadMGcPAfnVZh2NmZmYVxIlfhbn/2dVsfG2Hu3nNzMysy5z4VZj6xiaGD+7HMZNHZh2KmZmZVRgnfhVky7ad3PnUCk6bOpb+ff3VmZmZWdc4e6ggdy9cyeZtO5npbl4zMzPrBid+FaS+oYmRuwzgbZNGZB2KmZmZVSAnfhXila07uOuZlZxxyFjq+ijrcMzMzKwCOfGrEHc+tYKtO1p8b14zMzPrNid+FaK+oYlxuw7ksAm7ZR2KmZmZVSgnfhVg/avbuO+5VcycNo4+7uY1MzOzbnLiVwHuWLCc7TvD3bxmZmZWFCd+FWB2YzN7jxjMIeN3zToUMzMzq2BO/Mrc6le28uDzq5k5bRySu3nNzMys+5z4lbnb5jXTErib18zMzIrmxK/M1Tc0s9/oXThgzNCsQzEzM7MK58SvjDVv2MJjS9cya/oe7uY1MzOzojnxK2NzGpuJgJnTxmUdipmZmVUBJ35lrL6xmSl7DGPSqF2yDsXMzMyqgBO/MvXimldpeGm9J3WYmZlZyTjxK1Oz5zUBcMYh7uY1MzOz0nDiV6bqG5o5dMJw9tp9cNahmJmZWZVw4leGnl/5Ck83b2TWNHfzmpmZWek48StDsxubkOAMz+Y1MzOzEnLiV2YigvqGJo6cuDtjhg3MOhwzMzOrIk78yszTzZtYtGqzZ/OamZlZyTnxKzP1jU3U9RGnTR2bdShmZmZWZZz4lZGIYHZjE8dMHsmIXQZkHY6ZmZlVGSd+ZaTh5Q28tHaLb9FmZmZmPcKJXxmpb2iiX5141xR385qZmVnpOfErEy0twZzGZt6x/2h2HdQv63DMzMysCjnxKxNzl65j+cbXmDXd3bxmZmbWM5z4lYn6hiYG9uvDyQeNyToUMzMzq1JO/MrAjp0t3DqvmZMOHMOQAX2zDsfMzMyqlBO/MrD6lW1MGjXEF202MzOzHlWWiZ+kv5W0QFKLpBk5y98p6XFJ89LnE3PW3SNpoaQn08foAvudKGlLTpkf9lad2jN214H89vyjOdUXbTYzM7MeVK79ivOB9wFX5S1fDcyKiCZJU4E7gPE568+MiLkd7HtRRLyldKGamZmZVYayTPwi4mkASfnLn8h5uwAYKGlARGztxfDMzMzMKlJZdvV20vuBJ/KSvmvTLtwvKj9rfMM+kp6QdK+k43ohTjMzM7OykFmLn6Q7gUKD2i6PiFs62HYK8DXglJzFZ0bEMklDgd8BZwPX5W3aDEyIiDWSDgduljQlIjYW+IzzgPMAJkyY0NlqmZmZmZWtzBK/iDi5O9tJ2hO4CTgnIhbl7G9Z+rxJ0i+BI8lL/NLWwa3p68clLQL2B940LjAirgauBpgxY0Z0J1YzMzOzclJRXb2ShgNzgMsi4sGc5X0ljUxf9wNmkkwQyd9+lKS69PUkYD/ghd6I3czMzCxrZZn4SXqvpJeBo4A5ku5IV10ITAa+mHfZlgHAHZIagSeBZcCP0n29W9JX0u3fDjRKagBuAM6PiLW9VzMzMzOz7CjCvZgdmTFjRsyd29FVYszMzMyyJ+nxiJhRaF1ZtviZmZmZWek58TMzMzOrEU78zMzMzGqEEz8zMzOzGuHEz8zMzKxGOPEzMzMzqxFO/MzMzMxqhBM/MzMzsxrhxM/MzMysRjjxMzMzM6sRvmVbJ0haBSxtp8hIYHUvhVNuXPfaVcv1d91rUy3XHWq7/pVW970jYlShFU78SkDS3LbuiVftXPfarDvUdv1dd9e9FtVy/aup7u7qNTMzM6sRTvzMzMzMaoQTv9K4OusAMuS6165arr/rXptque5Q2/Wvmrp7jJ+ZmZlZjXCLn5mZmVmNcOLXDkmnSloo6XlJ/1Jg/fGSNkh6Mn18KWfdcEk3SHpG0tOSjurd6ItXZP0/K2mBpPmSfiVpYO9GX5yO6p6WOT6t9wJJ93Zl23LW3bpL2kvS3envfYGkz/Ru5MUr5ntP19VJekLS7N6JuLSK/N1X9DmvyLpX9flO0j/lnOfnS9opaffObFvuulv3ij7fRYQfBR5AHbAImAT0BxqAg/PKHA/MbmP7nwEfT1/3B4ZnXafeqj8wHlgMDErf/wY4N+s6lbjuw4GngAnp+9Gd3bacH0XWfRxwWPp6KPBsrdQ9Z/3ngF+2dV4o50ex9a/kc16Rv/uqP9/llZ8F3NWdbcvtUWTdK/Z85xa/th0JPB8RL0TENuB/gfd0ZkNJw4C3Az8GiIhtEbG+xyLtGd2uf6ovMEhSX2Aw0NQDMfaUztT9I8CNEfEiQESs7MK25azbdY+I5oj4S/p6E/A0yR/FSlHM946kPYEzgGt6Kd5S63b9q+CcV9R3T/Wf73J9GPhVN7ctN92ueyWf75z4tW088FLO+5cp/KUeJalB0m2SpqTLJgGrgGvTbp9rJA3p4XhLrdv1j4hlwDeBF4FmYENE/KGnAy6hztR9f2A3SfdIelzSOV3YtpwVU/fXSZoIHAo80kNx9oRi6/5t4PNAS8+G2WOKqX+ln/O6XfcaOd8BIGkwcCrwu65uW6aKqXvuuolU0PnOiV/bVGBZ/hTov5DcFmU68F3g5nR5X+Aw4AcRcSiwGai0sQ/drr+k3Uj+17QPsAcwRNJZPRhrqXWm7n2Bw0laeN4FfFHS/p3ctpwVU/dkB9IuJCfHSyJiY08F2gO6XXdJM4GVEfF4D8fYk4r57iv9nFfMd18L57tWs4AHI2JtN7YtR8XUPdlBBZ7vnPi17WVgr5z3e5LXfB8RGyPilfT1rUA/SSPTbV+OiNbs/waSk2IlKab+JwOLI2JVRGwHbgSO7p2wS6LDuqdlbo+IzRGxGrgPmN7JbctZMXVHUj+Sk+AvIuLGXoi3lIqp+zHAuyUtIekuOlHSz3s+5JIq9ndfyee8YupeC+e7Vh/ijW7erm5bjoqpe+We77IeZFiuD5L/3b1A8r+41kGfU/LKjOWNayEeSdLU3/r+fuCA9PUVwDeyrlNv1R94K7CAZKyLSAZ9X5R1nUpc94OAP6VlBwPzgamd2bacH0XWXcB1wLezrkdv1z2vzPFU5uSOoupfyee8In/3VX++S8vtCqwFhnR123J9FFn3ij3f9cUKiogdki4E7iCZ+fOTiFgg6fx0/Q+BDwCfkrQD2AJ8KNJfBHAR8AtJ/Ul+WB/t9UoUocj6PyLpBpKu4B3AE1TQVc87U/eIeFrS7UAjyZiuayJiPkChbTOpSDcUU3dJxwJnA/MkPZnu8l8jaQ0ue8V+75WuBPWv2HNeCf7NV/X5Li36XuAPEbG5o217twbdV0zdSVr5K/J85zt3mJmZmdUIj/EzMzMzqxFO/MzMzMxqhBM/MzMzsxrhxM/MzMysRjjxMzMzM6sRTvzMrCjpLazuyTqO3ibpg5IWSNoqKSQNzzqm9kiamMZ5btaxdIWkc9O4J2Ydi1k1cOJnVkUk3Shpu6RR7ZT5TPqHdFZvxlZNJE0GfgEsBz5Fcj2vze1uZGZWBnwBZ7Pqcj3JxUY/RHL/5ELOAlYDt5foM08p0X4qyfEk589LI+LJDsqamZUNt/iZVZc5JLcWKniTeEn7AzOA/43kvqLdJmkwQERsi4htxeyrAo1On9d3Z2NJgyT5/Gtmvc4nHrMqkiZgvwGOlLRfgSJnp8/XS+ov6UpJj0paK2mLpCcLjQGTtETSnZLeLukhSVuAr6br3jTGT9Klku6XtCodA/eMpH+UpLxy90h6XtJkSXdI2ixppaT/LJQYSXq/pAckbZK0UdJcSf+QV+YwSb+XtC6t01xJf9OZ4yepj6TPS1qYxt0k6Xu54/ckLQH+PX27OO02/2k7+2wdo3aypG9JaiLpFh4maXdJX5fUkNZnc3p8ZxbYT0i6RtJpkp6Q9Fp67D5SoOw4Sb9Nj9NaST8GhrUR3zGS/pSWfSV9fVQbdThR0tckLU/L/1bScEl9JX01PV5b0uPf5nCDnP0OSfe3KK3PGkl/lvSBTmw7Ky37qqT1km6RdFBemSvSuKdK+ln6m9go6VeSRhfYZ7d/O2aVwomfWfW5Pn0+s8C6jwDPRsSjJInA+cCfgS8Cl5G0YF0r6RMFtt0HuAV4GLgYuLudGD4HPA38v/T1M8A30vf5hgJ3ktzf9VLgQeCfgb+KQdK/ADcAg0mSzn8BHgdm5ZQ5Lt1+PEly9k/Aq8BNkj7cTrytvg98DXg2jfsmkmP0JyX3oAW4JI0D4LMkyfRVndj3t4Gj0v1fDmwDJgEfBP4AfB74MjAI+L2kQl3oRwA/BX4P/CPwCkkS/3rCI2kg8CfgPcA16T4nkdxQ/q9IejtwV7r+qyTHbF/gbknHFPj8bwJHknyPPwPeD/wY+AHwtnQf1wAzge904ph8n+QY1gMXpp//LPDW9jZKv8tbSI7VF0iO7bHAQ0rGX+b7ObAnye/8WpL7jP8h5zstxW/HrDJEhB9++FFlD+B54Lm8ZccAAXwhfV8HDCiw7R8LbLsk3fZ9BcrfA9yTt2xwgXLXkCQqA/K2DeD8vLJPAo/lvN8H2EGSIPXLK9t6z3GRJJv3AXW564EHgJday7ZxzKamsfwqb/kF6fJP5yz7QrpsYie+i3PTso8XiH1Abqzpsv7AAuCPecsjPQZTc5aNAbYC38hZdlFa9qM5y+rS4xLAuTnL5wLrgDE5y8YBG4BHC9ThPqBPzvJfAy0k/wnIX74d2KWDY7MO+F4nj9/E9H0/oBl4Lnf/wDRgJ/CbnGVXpNvelfeb+ES6/JOl+O344UclPdziZ1adfg5MlvS2nGVnkfyx+wVAROyMiK0Akvql3Y4jSf5ITpa0a94+l5O0gHUoIl5N99tX0m7pfu8BhgAH5BXfTpIU5rqXpBWq1ftIkpcrIm9sYkRE+nI6cGBav90kjUw/dwRwK0mLz/7thN3avfqNvOU/ImkJfVP3axf9qEDsWyNiJ4CkAZJGkLTE3gccXmAf90XE/JztV5C0puYeq5nAGt5o+SX9jL+a7CNpbPoZ16f7aS3bTPL7OULSmAJ1aMl5/zBJcvSTAsv7AhMK1CHXeuCtkvbqoFyuw4GxwA8i4pWcuBtJJiydVmCYwHdbj3Pqp/z1d1rsb8esYjjxM6tOP0+fzwJIu7Q+CDwQEYtbC0n6e0mNwGskycIq0rF7QH7i90JOktUuSadL+jOwhWSyySreSETyr3e3LCJ25C1bB+ye8761+25eOx/b+of5h+nn5T5ax+S9aVxXjonp8zO5CyMZN/k8SatjMRblL1DiUknPknwHq9N4z+fNxwlgaYFl+cdqb5LvKv+YLsx7PzF9foY3eyqvTKsX896v72D5bgX2netS4CBgqZLxpd+QVCjhzdUaU1tx7wLkjy/8q7qnCfjinH0V+9sxqxi+nItZFYqI5yU9DPydpEuA00mSg9dbgST9HUnLxxzgW8AKkta300nGXeX/x3BLZz5b0tEkY7YeBj4NLCMZz3YYyfi2/P3upGPquMjr+/1X4LE2ysxvY3lnPr9TSW87Ch2/zwP/SfK9XEmS+O0EPkoyHjNfW8dKea8LxdqZY5hfNn8/bX1+Z+J6k4i4UdIDJOM0TwY+Blwq6fKI+I/OBlvg8/Lj7uh49ORvx6ysOPEzq17XkwyefxdJy99W4Lc56z9E0uoxK7clT9KJRX7u35IkeidHxGs5+53U9iYdei59PgR4qI0yz6fPmyPizm58xpL0+UDgL60LJfUj6Ur9czf22ZEPkYyPPCd3oaSPFbHPJcDhkvrmtfrld1UuSZ8PLLCP1mWFWhhLKiJ4jK71AAADLElEQVRWkkwQ+bGSSwTNAa6U9M38rvHUkpwYb81bdyDJONLVBZY/3fom/U4nkozfg+J/O2YVw129ZtXr1yQJ2AUkY5nqIyL3unOtY7JePw+kY8yKSTpa99tCMiavdb8DSSYddNeNJK1KV6Z/tF8nvX6JmL+QJIiXqsDt0zpxeZHZ6fPn8pZ/nKTLsr6rQXdCC3nnYSWX4XlvEfucQzI2rfXSPUiqI+/4R8RykskdZ+de2iQd+3c2yeSOFfQQSXX540jTsaELSSZwDGlj07kk403Pl/R6GUlTgVOBW/PGGwJclB6DVueSdKXPSd8X+9sxqxhu8TOrUhGxVtKtQOt1yK7PK3ILyaSJ2ZJuJhnDdB7QRDJbtLt+T9JVfKek60ku1/L3JGPYuiUiFkv6MsllRB6R9BuSmadTSWahvi8iWiR9lGTm71OSfkLSOjSW5PIgB5NcqqStz5gv6Srgk5KGkUwUOJhkvN1fSFqlSu0WkmT2FySTX/Ym6R5/BnhLN/f5o3QfV0k6hOQyOe+n8HX8LiWZxf3ntO4CPgkM5M0JcKkNBZZJugloIBkLeihJon1b3n9SXhcROyR9jmQixoOSfkZSt4uATSSXysm3O8nlW24iGS96AUnX7bXpPov67ZhVEid+ZtXtepLEbw1wW+6KiLgubeG7gOSaay+SXKdtA+kfxO6IiHslnU0yXupbwEqSsYT3k/xh7e5+/13SIpLr6H2J5NImC0m6s1vLPCjpSJLrtZ1H0qqzgiSxKJQQ5Ps0SaL0cZLWozXA1cDl0TN3J/kPkku6nE2SnD1Lcj27/elm4hcRWySdRPKdnkfS6ntT+r4hr+x9admvkBwzgEeBMyOirS71UnkV+B+SsX1nkByHF0kmF329vQ0j4leSNpN8p18lqeM9wGUR8XyBTc4iGU/5byR/924EPtM6qz3dZ7G/HbOKoOjcJD0zM7OKIukKkgtY7xURL2ccjllZ8Bg/MzMzsxrhxM/MzMysRjjxMzMzM6sRHuNnZmZmViPc4mdmZmZWI5z4mZmZmdUIJ35mZmZmNcKJn5mZmVmNcOJnZmZmViOc+JmZmZnViP8DO+V1qo/A6kcAAAAASUVORK5CYII=\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",
      "/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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3yV5f3G8c83m7ACIayEvTeEjeIWUQQVpGLdW6utWq3WUcWtVX+OWgdV66pKFawiKooDVJbsvcIOAknYhJB1//7IwSIlIYeck+ecnOv9ep1Xc57nnHDRFnLxnOf+3uacQ0RERETCV5TXAURERESkYlToRERERMKcCp2IiIhImFOhExEREQlzKnQiIiIiYU6FTkRERCTMxXgdwGv16tVzzZs39zqGiIiIyFHNmTMn2zmXcvjxiC90zZs3Z/bs2V7HEBERETkqM1t/pOP6yFVEREQkzKnQiYiIiIQ5FToRERGRMKdCJyIiIhLmVOhEREREwpwKnYiIiEiYU6ETERERCXMqdCIiIiJhToVOREREJMyp0ImIiIiEORU6ERERkTCnQiciIiIS5lToRERERMKcCp2IiIhImFOhExEREQlzKnQiIiIiYU6FTkRCinOOwqJir2OIiISVGK8DHImZjQRGAx2APs652b7jpwOPA3FAPvAn59w3vnPfAY2A/b5vM8g5t61yk4vIsdidV8C01TlMXZXF1JVZbN2dx6ntGzCiZxontUshNlr/9hQRKUtIFjpgMTAceOWw49nAUOfcZjPrDEwCUg85f9HB8icioauo2LE4cxdTV2YxdVUWczfspKjYUT0umv6t6nFSuxQ+X7SFL5ZsoW71OIZ1a8yI9DQ6p9bCzLyOLyISckKy0DnnlgH/8xe3c27eIU+XAAlmFu+cO1CJ8UTkGGzdnceUlVl8vyqbH1ZlsSO3AIDOqbW47oSWnNA2hfSmdYiLKbkad//QTkxdmcX4uZm8O3MDb0xbR9sGNRiRnsa5PVJpUCvBy9+OiEhIMeec1xlK5fsY9fYjXXUzs/OB651zpx3y2mSgCBgHPOxK+c2Z2bXAtQBNmzbtuX79+qDkF4lkeQVF/LRue8lVuJXZrNi6B4B6NeI5oW09TmiTwvFt6lGvRvxRv9eu3AImLNzMuLmbmLdhJ1EGx7dJYUR6KoM6NqRaXHSwfzsiIiHBzOY453r9z3GvCp2ZTQYaHuHUPc65j32v+Y4jFDoz6wR8Qsl9chm+Y6nOuUwzq0lJoXvHOffW0XL06tXLzZ6tT2lFKso5R0bWXqaszGbqyixmrs0hr6CYuOgoejWvwwltUzihTQodGtWs0Mema7L2Mn5uJh/NyyRz535qxMcwpEsjhqen0qdFXX0kKyJVWsgVuvI4UqEzszTgG+AK59yPpbzvcqCXc+6mo/0aKnQix25XbgE/rC4pcN+vymLzrjwAWtarXlLg2tajX8tkEuMCf3dHcbFjxtocxs/N5LNFP5ObX0STutUY3iON4empNEuuHvBfU0TEa1Wi0JlZEjAFeNA5N+6Q18UASc65bDOLBd4DJjvnXj7ar6FCJ1J+hUXFLNj038UMCzbupNhBzfgYjmtdjxPapjCwTT2a1E2s1Fy5+YV8sXgL4+dm8mNGNs5B7+Z1GJ6expCujaiVEFupeUREgiWsCp2ZnQf8DUgBdgLznXNnmNm9wF3AqkNePgjYB0wFYoFoYDLwR+dc0dF+LRU6kbJt3rn/lwL3w6psducVYgZd05I4sU1JieveJImYEBktsnnnfj6al8m4uZtYk7WP+JgoBnVqyPD0VAa2rhcyOUVEjkVYFbrKpEIn8mv784uYsTbHt5ghi4ysfQA0qBXPCW1SOKFtCse3rked6nEeJy2bc44Fm3Yxfu4mPlmwmZ25BaTUjOe8HqkMT0+lfcNaXkcUEfGbCl0pVOgk0jnnWLF1zy+rUWet205+YTHxMVH0aVGXE9uWlLg29WuE7YKDA4VFfLt8G+PmZvLt8m0UFjs6Na7FiPQ0hnVvXK6VtiIioUCFrhQqdBKJdu0vYIrvCtz3q7LYurtklGOb+jV8ixlS6NuiLgmxVW8cSM7eA0xYsJlxczNZlLmLmCjjpHYpDE9P49QO9YmPqXq/ZxGpOlToSqFCJ5Fm2c+7uez1WWzbc4Da1WI5vnU9Tmhbj4FtUmicVM3reJVq5dY9jJu7if/My2Tr7pL/PoZ2a8Tw9DR6NEkK2yuSIlJ1qdCVQoVOIsnMNTlc/dZsqsfF8Oyo7vRuXpfoKJWWomLHD6uzGT93E5OWbCGvoJiW9aozPD2V89LTSI2woisioUuFrhQqdBIpvlyyhZvem0danWq8fVVflZRS7Mkr4PNFW/hw7iZmrd2OGfRvmczw9DTO7NyQ6vEhuWOiiEQIFbpSqNBJJPj3Txv58/iFdElL4p+X96ZuiK9QDRUbcnL5aF4m4+dtYn1OLolx0Qzu3JAR6Wn0b5lMlK5uikglU6ErhQqdVGXOOV6esoYnvljOwDb1ePninrrCdAycc8xZv4Nxczfx6YKf2XOgkMa1EzgvPZXh6Wm0SqnhdUQRiRAqdKVQoZOqqrjY8chny3jth7UM69aYp0Z2Iy5GQ3UrKq+giK+WbmXc3E1MXZlFsYPuTZIYkZ7K0G6NSUrU1U8RCR4VulKo0ElVVFBUzB0fLuSjeZlcPqA5953dUR8PBsG23Xl8PH8z4+ZuYvmWPcRFR3FK+/qM6JnGSe1SiNWuFCISYCp0pVChk6omN7+Q3/1rLt+tyOL2QW258eTWGr8RZM45lv68m3FzMvl4fiY5+/JJrh7HsO6NGZGeRqfGtfS/gYgEhApdKVTopCrZmZvPFW/8xIKNO3nkvC5c2Kep15EiTkFRMVNXZjFu7iYmL91GflEx7RrULBmB0iOV+rUSvI4oImFMha4UKnRSVfy8az+XvjaL9dtzeX5UdwZ3buR1pIi3K7eACQtLPpKdt2EnUQYD26QwPD2VMzo1rJI7cYhIcKnQlUKFTqqC1dv2culrM9mdV8g/Lu1F/1bJXkeSw6zJ2sv4uZl8NC+TzJ37qRkfw1ldGjGiZxq9m9fRR7IiUi4qdKVQoZNwN2/DDq584yeio6J444redE6t7XUkKUNxsWPG2hzGzcnk88U/k5tfRJO61RjeI43ze6bRpG6i1xFFJISp0JVChU7C2ZSVWVz/9hxSasbz9lV9aJZc3etI4ofc/EK+WLyF8XMz+TEjm9ioKF68KJ3TOjbwOpqIhKjSCp3W1IuEqY/nZ3LVGz/RvF51Pryhv8pcGEqMi2F4ehrvXN2XH+48hfaNanL9O3P4fNHPXkcTkTCjQicShv7541pufn8+6c3qMPa6ftSvqZWT4S41qRrvXN2Xrmm1uem9eXw8P9PrSCISRlToRMKIc46nJq3ggQlLGdSxAW9d2YdaCbFex5IAqZUQy1tX9aVXszrcOnY+H87Z5HUkEQkTKnQiYaKo2HH3R4t44dvVjOrdhBcvStfYiyqoRnwMb1zRhwGt6vGnDxfw7swNXkcSkTCgQicSBvIKivjdv+bw3qyN3HRyax4b3oUYbStVZVWLi+bVy3pxUtsU7v5oEW9OW+d1JBEJcfqJIBLiducVcNnrs5i0ZCv3D+3I7We008yyCJAQG83Ll/RkUMcG3P/JEv4xdY3XkUQkhKnQiYSwbXvyGPXKDOas38Fzo7pzxXEtvI4klSg+Jpq/X5TOkK6NeOSzZbzwzSqvI4lIiIrxOoCIHNn6nH1c8tossvYcKPn4rV19ryOJB2Kjo3jugu7ER0fx1JcrOVBYzB9Pb6urtCLyKyp0IiFoyeZdXPb6TxQWF/PuNX3p0bSO15HEQzHRUTw5shux0VH87ZvV5BcW8+cz26vUicgvVOhEQsyMNTlc8+ZsaibE8P61/Wldv6bXkSQEREcZjw3vQlxMFK9MXcOBwmLuH9pRpU5EABU6kZDyxeIt/OH9eTStm8hbV/ahcVI1ryNJCImKMh48pxNxMVG89sNa8ouKeficzkRFqdSJRDoVOpEQ8f6sDdz90SK6NUni9ct6U6d6nNeRJASZGfcO6UBcTBQvfZdBfmExT4zoSrRKnUhEU6ET8Zhzjhe/y+DJSSs4sW0KL12cTmKc/mhK6cyMO85oR3xMFM9OXkVBUTFPj+ym2YQiEUw/NUQ8VFzsePDTpbwxbR3ndm/8y43vIkdjZtxyWlviYqL46xcryC8s5rlRPYiL0f9/RCKRCp2IR/ILi7n9gwV8smAzVx3fgnvO6qB7ocRvvzupNXHRUTw8cRkF/5rD3y9KJz5GW8KJRBr9U07EA/sOFHL1W7P5ZMFm7hzcnnuHqMzJsbt6YEseOqcTk5dt45q35pBXUOR1JBGpZCp0IpVs+758fvvqTH5YlcUTI7pww0mtNHpCKuyS/s15YkQXvl+VxZVv/ERufqHXkUSkEqnQiVSizJ37GfnyNJb9vJuXL+7JBb2beh1JqpALejfl/37TjRlrcrjs9VnsySvwOpKIVBIVOpFKsmrrHs5/aRrb9hzg7Sv7MKhTQ68jSRV0Xo80nr+wB3M37OSS12axa79KnUgkUKETqQRz1u9g5CvTKSx2jL22P31bJnsdSaqws7s25sWL0lmyeRcXvTqDHfvyvY4kIkGmQicSZN+u2MbFr84kqVos464fQMfGtbyOJBHgjE4NGXNJL1Zu3cuF/5hB9t4DXkcSkSBSoRMJov/My+SaN2fTMqU6H1w/gKbJiV5Hkghycvv6vHZZL9bl7GPUmBls253ndSQRCRIVOpEgee2Htdwydj69m9fl/Wv7kVIz3utIEoEGtknhjSv6sHnnfi4YM4Ofd+33OpKIBIEKnUgQPDt5JQ99upQzOzfkn1f0pmZCrNeRJIL1a5nM21f1IXvPAX7zynQ2bs/1OpKIBJgKnUiALdi4k+e+XsXwHqm88Nt0EmI1tV+817NZXd65ui+7cgu44JXprMve53UkEQkgFTqRACoudtz/yRLq1YjngXM6Ea3dHySEdGuSxHvX9mN/QRG/eWU6q7ft9TqSiASICp1IAI2bu4n5G3fy58Ht9TGrhKROjWvz/rX9KXYwasx0VmzZ43UkEQkAFTqRANmdV8ATXyynR9MkzuuR6nUckVK1a1iTsdf1IzrKGDVmOoszd3kdSUQqSIVOJECem7yKnH35PDisM1H6qFVCXKuUGvz7uv4kxsXw23/MYP7GnV5HEpEKUKETCYBVW/fw5rR1jOrdhC5ptb2OI1IuzZKrM/a6fiQlxnHxqzOZvW6715FE5Bip0IlUkHOO0ROWkBgXze2D2nkdR8QvaXUSGXtdP+rXjOfS12cxPSPH60gicgxU6EQqaNKSLfy4OofbBrUjuYaGB0v4aVS7Gu9f24/UpGpc/s9ZTF2Z5XUkEfGTCp1IBezPL+KhT5fRvmFNLurb1Os4Isesfq0E3r+2Hy1TanD1m7P5ZvlWryOJiB9U6EQq4OUpGWTu3M/oYZ2IidYfJwlvyTXiee+avrRrWJPr3p7DF4u3eB1JRMpJP4FEjtHG7bm8PCWDs7s2ol/LZK/jiAREUmIc/7qmL51Ta3Pju3OZsGCz15FEpBxU6ESO0cMTlxJlxj1DOngdRSSgaiXE8vZVfenZtA43vz+PcXM2eR1JRI5ChU7kGHy/KotJS7Zy0ymtaVS7mtdxRAKuRnwMb1zZm/6tkrn9wwW8P2uD15FEpAwqdCJ+Kigq5oEJS2mWnMjVA1t4HUckaBLjYnjtst6c0CaFP49fxFvT13kdSURKoUIn4qc3p61j9ba93Hd2R+Jjor2OIxJUCbHRjLm0J6d1aMB9Hy/h1e/XeB1JRI5AhU7ED9v25PHs5FWc1C6FU9rX9zqOSKWIj4nmpYvTGdKlEQ9PXMbfv13tdSQROUyM1wFEwskTn6/gQGER953dETPt1yqRIzY6iudGdSc22nhy0goOFBZz62lt9OdAJESo0ImU05z1Oxg3dxPXn9iKlik1vI4jUulioqN4+jfdiY2O4vmvV5FfWMydg9up1ImEABU6kXIoLnaM/mQJDWrF8/tTWnsdR8Qz0VHGEyO6EhcTxctTMnTFWiREqNCJlMO/Z29kUeYunhvVnerx+mMjkS0qynj43M7ExUTxzx/XkV9YzEPndCYqSqVOxCv6ySRyFLtyC/jrpBX0bl6HYd0aex1HJCSY2S8rvV+ekkFBUTGPDe9KtEqdiCdU6ESO4pnJK9mZm8/oYX30sZLIIcyMOwe3Iy7mv/fUPTWym/Y1FvGACp1IGZZv2c3bM9bz275N6dS4ttdxREKOmfHH09sSHxPFk5NWkF9UzHOjehCrUidSqVToRErhnOP+j5dQMyGG205v53UckZB248mtiY+J4uGJy8gvnMvfL+qhwdsilUj/hBIpxcRFPzNz7XZuH9SOOtXjvI4jEvKuHtiSB8/pxORlW7nu7TnkFRR5HUkkYqjQiRxBbn4hj0xcRqfGtbiwT1Ov44iEjUv7N+ex4V2YsjKLq978idz8Qq8jiUQEFTqRI3jx2wx+3pXHA8M6adWeiJ8u7NOUp87vxvSMHC5//Sf2HlCpEwk2FTqRw6zP2ceYqWs4t3tjejWv63UckbA0omcaz43qwZwNO7jktZns2l/gdSSRKk2FTuQwD326lNho466zOngdRSSsDe3WmL//Np3Fmbu4+NWZ7MzN9zqSSJWlQidyiG9XbGPysm38/tQ2NKiV4HUckbA3uHNDXrmkJyu27mHUmBnk7D3gdSSRKqnUQmdmxWZW5O8jEKHMbKSZLfFl6HXI8dPNbI6ZLfL95ymHnIszszFmttLMlpvZiEBkkciRX1jMQxOW0rJeda48roXXcUSqjFPaN+DVS3uxNnsfo8bMYNvuPK8jiVQ5Zc2hexBwhx07F+gMTAJWAAa0AwYBi4CPA5RrMTAceOWw49nAUOfcZjM7mCPVd+4eYJtzrq2ZRQG6+Un88vqPa1mTvY83ruhNXIwuXosE0gltU3jjij5c9eZPjBozg39d05dGtat5HUukyii10DnnRh/63MyuBBoCnZ1zKw471wH4FtgQiFDOuWW+73v48XmHPF0CJJhZvHPuAHAl0N73umJKyp9IuWzdncffvl7FaR0acFK7+l7HEamS+rdK5q0r+3D5P3/iN69M592r+9GkbqLXsUSqBH8uQ9wBvHB4mYNfCtjfgTsDFawcRgDznHMHzCzJd+whM5trZh+YWYNKzCJh7rHPllFQ7PjL2VoIIRJMvZrX5Z2r+7Irt4BRY2awPmef15FEqgR/Cl0zYH8Z53N9rykXM5tsZouP8DinHO/tBDwBXOc7FAOkAT8659KB6cBTZbz/WjObbWazs7KyyhtZqqif1m3nP/M3c+3AljRLru51HJEqr3uTJN69ph+5+YVc8IoWSogEgj+FbiVw7SFXw35hZnWAaym5r65cnHOnOec6H+FR5n14ZpYGfARc6pzL8B3OoaRQfuR7/gGQXsavPcY518s51yslJaW8kaUKKiou2a+1ce0EfndyK6/jiESMzqm1efuqvmzfl8+d4xbh3OG3bIuIP/wpdHcDrYBVZvZ/Zna9mV1nZs9QUvZaUrIwIWh8ZXIicJdz7seDx13J3wQTgJN8h04FlgYzi1QN783awNKfd3P3kA4kxpW1RkhEAq1zam3uPLM9k5dt5d1ZAbkFWyRilbvQOecmAmdQsvDhFuBF4CXgZt+xM32vqTAzO8/MNgH9gYlmNsl36iagNfAXM5vvexy8g/1OYLSZLQQuAW4LRBapunbsy+epL1fQv2UyQ7o08jqOSES6YkBzBrapx0OfLmX1tr1exxEJW3Ysl7l9Cw6aUzK2ZK1zbmuAc1WaXr16udmzZ3sdQzxw738W8d6sjUz8w/G0b1jL6zgiEWvb7jzOeHYqjZOq8dHvjtPYIJEymNkc51yvw48f058a59xW59xM59yMcC5zErmWbN7FuzM3cEm/ZipzIh6rXyuBJ0Z0Zcnm3Tz9VblvxRaRQ/hV6Mwsycwe961G3Wdme31fP3qkxRIiocg5x+hPlpCUGMetp7X1Oo6IAIM6NeS3fZsyZuoapq3WGFERf5W70JlZKjCPknl0ULI44XNKdpP4MzDXzBoHPKFIgH2yYDM/rdvBHWe0o3ZirNdxRMTn3iEdaFGvOn/89wJ25uZ7HUckrPhzhe4xoAFwtm+8yG+ccyOdc12AIb5zjwUjpEig7DtQyKOfLaNrWm1+06uJ13FE5BCJcTE8P6oHOfsOcNd4jTIR8Yc/hW4w8Jxz7rPDTzjnPgf+BpwZqGAiwfC3b1azdfcBRg/rRFSUHf0NIlKpOqfW5rZB7fh88RY+mLPJ6zgiYcOfQlcTyCzj/Cbfa0RC0pqsvbz2wxpGpKeR3rSO13FEpBTXDmxJ/5bJjP5kCeuytTWYSHn4U+hWAOeb2f+8x8yigfPxY6cIkcrknOPBT5cSHxPNnWe28zqOiJQhKsp4+jfdiI2O4uax8ykoKvY6kkjI86fQPQ+cAHxjZueYWXvf41xgMjAQeDYYIUUq6utl2/huRRa3nNaG+jUTvI4jIkfROKkaj57XhQUbd/L816u8jiMS8sq915Fz7nXfrgz3A+MPOWXAAeBu59wbgY0nUnF5BUU8+OlSWtevwWUDmnsdR0TKaUjXRny3Io2/f7uagW1S6NOirteRREKWX3PonHOPA2nARZTs7Xo3cCGQ6px7IvDxRCrutR/WsmF7LqOHdiI2WhPoRcLJ/cM60aRuIreOnc+u/QVexxEJWX7/dHPO5Tjn3nfOPeF7jHXObQ9GOJGK2rxzPy98s5rBnRpyfJt6XscRET/ViI/h2Qu6s2V3Hvd9vNjrOCIhy+9CZ2aDzewFM5toZp/6vh4UjHAiFfXoZ8sodo57hnTwOoqIHKMeTetwy6lt+Hj+Zv4zr6xhCyKRy5+dIuLM7GNKdoj4HdAb6Ov7+nMz+4+ZxQUnpoj/pmfk8OnCn7n+xFY0qZvodRwRqYDfndya3s3r8Jf/LGbj9lyv44iEHH+u0N0PDAWeBuo75+o751KAFOApYBjwl8BHFPFfYVExD0xYQmpSNW44qZXXcUSkgqKjjP/7TXcAbh07n0KNMhH5FX8K3W+Bd5xzdzjnftk52XdP3Z3AO8DFgQ4ociz+NXMDy7fs4S9ndyAhNtrrOCISAE3qJvLweZ2ZvX4HL32X4XUckZDiT6FrDEwr4/x0oFHF4ohUXM7eAzz95QqOb12PMzo19DqOiATQOd1TObd7Y579ehXzNuzwOo5IyPCn0G0G+pVxvi/wc8XiiFTcU1+uIDe/iNHDOmKm/VpFqpoHz+1Mw1oJ3DJ2PnsPFHodRyQk+FPo3gMuMbOHzeyXjTDNrI6ZPQRcArwb6IAi/li4aSfv/7SRywY0p3V9bS0sUhXVSojl2VHd2bg9lwc+WeJ1HJGQ4E+hewCYRMkw4Wwz22JmW4Bs4B7fuQcDH1GkfIqLHfd/soTk6vHcfFobr+OISBD1bl6Xm05uzQdzNjFxoT4cEvFn668DwFlmdjYwBGjuO7UOmOCc+yzg6UT8MH5eJvM27OTJ87tSKyHW6zgiEmS/P7UNU1dlc9f4hfRomkTjpGpeRxLxzLHsFPGpc+4G59yZvscNKnPitT15BTz++XJ6NE1iRHqa13FEpBLERkfx3KjuFBU7/vjv+RQVO68jiXhGG1tKlfD816vI2XeAB4Z1IipKCyFEIkWz5OqMHtaJGWu284/v13gdR8Qz5f7IFcDMTgeuBloCdYHDf3I655ymuEqlWr1tD//8cR0X9GpC17Qkr+OISCU7v2ca363I4ukvV3Bcq3p0SavtdSSRSufP1l+3Al8AJwKZwFRgymGPqUHIKFIq5xwPTFhKtbhobj+jnddxRMQDZsYj53WmXo14bh47j9x8jTKRyOPPFbpbKSltg51z+UHKI+KXSUu28v2qbO4f2pF6NeK9jiMiHklKjOPp33Tjoldn8vDEZTx6XhevI4lUKn/uoasHjFWZk1CRV1DEwxOX0q5BTS7p18zrOCLisQGt6nHdCa14d+YGvlyyxes4IpXKn0I3l5J750RCwitT1rBpx35GD+tETLTW94gI/PH0tnROrcWd4xaybXee13FEKo0/PwVvoWSniNODFUakvDbtyOXF71YzpGsj+rdK9jqOiISIuJgonhvVg/0FRdz2wQKKNcpEIkSp99CZ2ZdHOLwH+MLM1lEyULjosPPOOXdGwNKJlOKRicswg3vO6uB1FBEJMa1SavCXsztyz0eL+ee0dVx1fAuvI4kEXVmLItoCR/qnzQZKruzp41fxxI+rs/l88RZuO72tJsOLyBH9tk9Tvl2exROfL2dAq2Q6NKrldSSRoCq10DnnmldiDpFye+arlaTVqcY1J+jfFCJyZGbGEyO6MPi577n5/Xl8ctPxJMRGex1LJGh0J7mElbXZ+5i9fgcX92umv5xFpEzJNeJ5amQ3Vm7dy+OfL/c6jkhQqdBJWBk3ZxNRBuf1SPU6ioiEgRPbpnDlcS14Y9o6vl2+zes4IkFTaqEzs2IzKzSzuEOeFx3lofHcEjRFxY5xczdxQtsUGtRK8DqOiISJOwa3o33DmvzpwwVk7z3gdRyRoChrUcSDlCyKKDzsuYgnpmfk8POuPO4ZopWtIlJ+CbHRPDeqB0Nf+IE7PlzIa5f1wuzwrchFwltZiyJGl/VcpLJ9OGcjtRJiOK1DA6+jiEiYadewJnef2Z7RE5byzoz1XNK/udeRRAJK99BJWNidV8AXS7YwrHtjLYYQkWNy2YDmnNg2hYcnLmPV1j1exxEJqLIGC59wLN/QOTf12OOIHNlnC38mr6CY83s28TqKiIQpM+PJkV0589nv+cP78/nPjQOIj9E/EKVqKOseuu/w7545871efzok4D6cs4nW9WvQLa2211FEJIzVr5nAX8/vylVvzuapSSu4Z0hHryOJBERZhe7kSkshUoaDs+f+fGZ73cgsIhV2aocGXNKvGf/4fi0ntq3P8W3qeR1JpMLKWhQxpTKDiJRGs+dEJNDuPqsD09fk8Md/z2fSLSdQp3qc15FEKuSYFkWYWRszO87M9PmXBJVmz4lIMFSLi+a5Ud3ZkZvPn4C3jpYAACAASURBVMcvxDlN5ZLw5lehM7MLzGw9sByYCvT0Ha9nZqvMbGQQMkoEOzh7bkR6mtdRRKSK6dS4Nnec0Z5JS7Yy9qeNXscRqZByFzozOwd4D9gA/IWSRRAAOOeygWXAJYEOKJHtwzkbqZkQw+kdNXtORALvquNbcFzrZB6YsJQ1WXu9jiNyzPy5QncvMNU5NxB45QjnZwLdApJKhENmz3XT7DkRCY6oKOPpkd2Jj43i5vfnk19Y7HUkkWPiT6HrBPy7jPNbAF1GkYD57+w5fdwqIsHTsHYCjw/vyqLMXTw7eaXXcUSOiT+FLg8o66705sDOCqUROcSHczbRKqU63ZskeR1FRKq4wZ0bMqp3E16aksGMNTlexxHxmz+F7gfgwiOd8K12vRL4JhChRA7Onju/ZxPNnhORSvGXszvSPLk6fxw7n125BV7HEfGLP4VuNNDJzL4FhvuO9TKzm4D5QC3gocDGk0il2XMiUtmqx8fw7AXd2bbnAHf/Z5FGmUhYKXehc87NBc4AGvLfRRGPA88D+cAZzrllAU8oEefg7LmBbVJoWFuz50Sk8nRrksStp7dl4sKfGT830+s4IuXm1xw659z3zrkOQA/gAko+gu0NtHfOTTOzGkHIKBHm4Ow5LYYQES9cf2Ir+rSoy30fL2Z9zj6v44iUiz9z6F44+LVzboFz7gPn3Fjn3BznnDOzJGByUFJKRNHsORHxUnSU8cwF3YmKMm4dO5/CIo0ykdDnzxW668zs0SOdMLO6lCyIaBeQVBKxNHtOREJBalI1Hj2vC3M37OSFb1d7HUfkqPwpdFcBd5rZXYceNLP6lGwD1gwYFMBsEoE0e05EQsXQbo0Znp7K81+vYs767V7HESmTP4si3gL+ADxiZjcCmFka8D1QHzjFOfdTUFJKxNDsOREJJQ8M60RqnWrcMnY+e/I0ykRCl7+LIv4O3AM8Z2Z/BqYANYGTnHMLgpBPIohmz4lIqKmZEMuzF/Rg88487v9kiddxRErlV6EDcM49BvwVeBSIAQY655YGOphEHs2eE5FQ1LNZHX5/SmvGz83kkwWbvY4jckQxpZ0wszFHee8eYDkl99UdPOacc9cFKJtEEM2eE5FQdtPJrZm6Mot7PlpEz2Z1SE2q5nUkkV8ptdABV5fj/acf9twBKnTit4Oz5+4+q4PXUURE/kdMdBTPXtCDs57/nlvHzue9a/oRHaVbQyR0lPqRq3Mu6hgemjMhx0Sz50Qk1DVNTuTBczoxa+12Xp6S4XUckV/x+x46kUDT7DkRCRfn9UhlaLfGPPPVShZs3Ol1HJFfqNCJ5zR7TkTChZnx8LmdqV8znlvGzmffgUKvI4kAZRQ6M1trZhlmFnvI8zVHeegatPhNs+dEJJzUrhbL/13QnXU5+3joUw15kNBQ1qKIKZQscig+7LlIwBycPXfn4PaaPSciYaNfy2RuOLEVL36XwUntUhjcuZHXkSTClVronHOXl/VcJBA0e05EwtUtp7Xlh9XZ/Hn8Iro3qaORS+Ip3UMnntHsOREJZ3ExUTx7QXcOFBRz2wfzKS7Wh1jinbLuoWt6LI/KDC/h7eDsOS2GEJFw1TKlBvcP7ciPq3N47Ye1XseRCFbWPXTrOLZ75jR3QspFs+dEpCq4oHcTvl2xjb9OWs6A1sl0alzb60gSgcoqdFeiRRASJAdnz41IT9PsOREJa2bG48O7Mvi5qdz8/nwm3HQ81eL095pUrrIWRbxRiTkkwmj2nIhUJXWqx/H0yO5c/NpMHv1sGQ+d29nrSBJhQnJRhJmNNLMlZlZsZr0OOX66mc0xs0W+/zzFd7ymmc0/5JFtZs969zuQo9HsORGpao5vU49rBrbg7Rnr+XrZVq/jSIQJyUIHLAaGA1MPO54NDHXOdQEuA94GcM7tcc51P/gA1gPjKzOwlN/B2XPn92yi2XMiUqXcfkY7OjSqxR0fLmTbnjyv40gECclC55xb5pxbcYTj85xzm31PlwAJZhZ/6GvMrA1QH/g++EnlWGj2nIhUVfEx0Tw/qjt7DxTypw8W4pxuRZfKEZKFrpxGAPOccwcOO34hMNaV8afIzK41s9lmNjsrKyuoIeXXNHtORKq6Ng1qcu+QDkxZmcWb09Z5HUcihGeFzswmm9niIzzOKcd7OwFPANcd4fQo4L2y3u+cG+Oc6+Wc65WSknJsvwE5Jpo9JyKR4OJ+zTilfX0e/Xw5K7bs8TqORADPCp1z7jTnXOcjPD4u631mlgZ8BFzqnMs47Fw3IMY5NyeI0aUCNHtORCKBmfHX87tSKyGGm9+fR15BkdeRpIord6Erxy4RTcwsxYJ4l7uZJQETgbuccz8e4SUXcpSrc+KdPb7Zc8O6NdbsORGp8urViOfJkd1YvmUPf/3if24LFwkof67QrQPWlvFYB2wB9prZZ2bW71hDmdl5ZrYJ6A9MNLNJvlM3Aa2BvxwyoqT+IW/9DSp0IeuzRZo9JyKR5eR29bl8QHNe/3EtU1bqnm0JHivvChwzuwL4PdAM+BewCjCgLSVXxtYCbwBtgIuBGsBpzrmQXm3aq1cvN3v2bK9jRISRL09j+758Jv/xRI0rEZGIkVdQxLAXfmBHbgFf3DyQ5BrxR3+TSCnMbI5zrtfhx/25QtcQSATaOOf+4Jz7m3PueefcTUA7oCYQ75y7GWgPbAPuD0B2qQLWZe/jp3WaPScikSchNprnRvVgV24Bd45bpFEmEhT+FLrrgX8457YffsI5lw28SslHojjnsoDXgN6BCCnhb9xczZ4TkcjVoVEt7jyzPZOXbeXj+ZuP/gYRP/lT6OoDsWWcj6HkKt5Bmyhjr1iJHMXFjnFzNHtORCLbFQOa0zWtNo99vox9Bwq9jiNVjD+FbiFwo5k1PfyEmTUDbgQWHHK4DfBzxeJJVTB9TQ6bNXtORCJcVJRx/9BObN19gJe+yzj6G0T84M8VtNuAScAKM5tIyaIIKCluQ3xfjwIwswTgIkpGjEiE+3DOJs2eExEBejarw3k9Uhnz/Rou6N2EJnUTvY4kVUS5r9A5534ABgBfAmcCd/oeZ1JS9Pr5XoNzLs85l+acO9JODhJB9uQV8PninzV7TkTE587B7YmJMh6ZuMzrKFKF+LVThHNugXPuHEpWtDb2PWo65851zi0o+90SiTR7TkTk1xrWTuDGk1vzxZItTFud7XUcqSKOaesv51wxUAQU+r4WOaIP52yiVUp1ujdJ8jqKiEjIuOr4FqTVqcYDE5ZSWKQfo1JxfhU6M2tpZu+b2S5KdoXYama7zOxdM2sZnIgSrjR7TkTkyBJio7l3SAdWbN3De7M2eB1HqoByL4ows/bAj0Bt4AtgKSU7RXQARgKDzOx459zyYASV8KPZcyIipTujU0MGtErm6a9WMrRbY5IS47yOJGHMnyt0jwPFQA/n3NnOuTucc39yzp0N9AAc8GgwQkr4OTh77njNnhMROSIz476hHdm9v4BnvlrpdRwJc/4UuhOBvznnFh1+wjm3GHgBODlQwSS8afaciMjRtW9Yi4v7NeOdmRtYsWWP13EkjPlT6OKA3WWc3+V7jcgvs+cGafaciEiZbj2tLTXiY3jw0yXa51WOmb87RVxmZtUOP+E7dpnvNRLhDs6eG6rZcyIiR1Wnehy3DWrLj6tz+HLpVq/jSJjyp9A9CnQF5pnZzWY22Pe4BZgHdAEeCUZICS+aPSci4p/f9mlKuwY1eXjiUvIKiryOI2HIn50iJgAXA7WAZyjZ1msi8H++Yxc75z4NRkgJLx/O2UTLlOr00Ow5EZFyiYmO4r6hHdm4fT+v/bDW6zgShvzdKeI9oCnQH/it79EfaOKcez/w8STc/Hf2XJpmz4mI+OG41vU4o1MD/v7tarbsyvM6joQZv3eKcM4VOudmOufG+h4znXO6PizAf2fPDe+hj1tFRPx1z1kdKSx2/PULjXQV/5Q6WNjMmh7LN3TOaeR1hNLsORGRimmanMg1A1vw928zuLh/M9Kb1vE6koSJsq7QrQPWHsNDIpRmz4mIVNzvTmpNg1rxPPDJEoqLNcZEyqesrb+upGT3B5Fy0ew5EZGKqx4fw5/PbM+tYxcwbu4mRvZq4nUkCQOlFjrn3BuVmEPC3MHZc8PT0zR7TkSkgs7plspb09fzxBcrGNy5ITUTYr2OJCHO70URIkei2XMiIoETFWWMHtqJ7L0HeOHb1V7HkTCgQicBodlzIiKB1a1JEuf3TOP1H9ayNnuf13EkxKnQSYVp9pyISHDcMbgdcdFRPDJxmddRJMSp0EmFafaciEhw1K+ZwO9PbcPkZVuZujLL6zgSwlTopEI0e05EJLiuOK45zZMTefDTpRQUFXsdR0KUCp1UiGbPiYgEV3xMNPcO6cjqbXt5e/p6r+NIiPKr0JlZnJldZWb/MrOvzKyH73gdM7vUzPRTPcJo9pyISPCd2qE+J7RN4ZnJK8nZe8DrOBKCyl3ozKwOMAP4B3AWcApwcE+SXcBDwE2BDiih6+DsuaHdGmv2nIhIEJkZ953dgf35RTz91Uqv40gI8ucK3eNAW2AQ0Ab4ZTmjc64YGA8MDmg6CWmaPSciUnla16/Jpf2b896sDSzZvMvrOBJi/Cl0w4DnnXOTOfKWYKuBZgFJJWFBs+dERCrXzae1oU5iHA9MWIpz2p1T/sufQlcHWFPG+WggrmJxJFxo9pyISOWrXS2W2we1Y9ba7Xy2aIvXcSSE+FPo1gJdyjh/IrCiYnEkXGj2nIiINy7o3YQOjWrx6GfL2J9f5HUcCRH+FLp3gKvN7JRDjjkAM7sVOBd4I3DRJFRp9pyIiHeio4zRQzuSuXM/Y6aW9cGZRBJ/Ct0TwNfAV5SsdnXAC2a2FXga+Bj4W8ATSsjR7DkREW/1bZnMkK6NeGnKajJ37vc6joSAchc651yhc24YcDGwGFjue/8s4BLn3HCnOzQjgmbPiYh4764z2+McPP75cq+jSAiI8fcNzrn3gPeCkEXCwMHZc8PT0zR7TkTEQ2l1Ern+xFY89/UqLunXjD4t6nodSTykrb/EL5o9JyISOq4/sRWNayfwwIQlFBXrQ7JIVuoVOjN7/Ri+n3POXVWBPBLiNHtORCR0VIuL5q6zOvD79+bxweyNjOrT1OtI4pGyPnI9hSMPEC6L/nlQhR2cPXfH4HaaPSciEiLO7tqIt6ev58lJKzizSyNqV4v1OpJ4oNRC55xrXok5JAxo9pyISOgxM+4b2pGhL/zA375exb1nd/Q6knhA99BJuWj2nIhI6OqcWptRvZvwxrR1rN621+s44gEVOikXzZ4TEQlttw9qR7W4aB6euNTrKOKBshZFrAWKgfbOuQLf86PdI+ecc60CGVBCg2bPiYiEtuQa8dx8ahsenriMb5Zv5ZT2+vs6kpS1KGIKJQWu+LDnEmE0e05EJDxc2r85787awEOfLuP41inExeiDuEhRVqH7A5DrnCsCcM5dXimJJORo9pyISHiIi4nivrM7cvk/f+KNaWu59gR9aBYpyqruO4CRB5+Y2etm1jf4kSTUaPaciEj4OKldfU5pX5/nv15N1p4DXseRSlJWoSsA4g95fjmgqh9hDs6eO79nmmbPiYiEiXuHdOBAYRFPTtI+r5GirEK3GrjIzLqY2cHR08lm1rSsRyVklkqk2XMiIuGnZUoNrjiuBR/M2cTCTTu9jiOVoKxC9yBwIjAfOLjC9Vnf12U9pIooLnaMn5up2XMiImHo96e0Jrl6HA9MWIpzWtNY1ZW1U8S/zWwOcBLQAHgY+BBYUDnRxGsz1uSQuXM/d57Z3usoIiLip5oJsdxxRnvuGLeQTxZs5pzuqV5HkiAqa5UrzrkMIAPAzK4B3nHOfVIZwcR7mj0nIhLezu+Zxtsz1vPYZ8s5vWMDEuPK/LEvYazcA2qccy1U5iLHnrwCPlv8M0O7NdbsORGRMBUVZYwe1pEtu/N46bsMr+NIEPk1cdDMYszsejP71MyWmNli39fXmplqfxXy+aItmj0nIlIF9GxWl3O7N+aVqWvYuD3X6zgSJOUudGZWG5gOvAgMAPYB+4H+wMvANN9rpArQ7DkRkarjzjPbE23Go58t8zqKBIk/V+geAXoANwENnHN9nHO9KVkwcaPv3MOBjyiVbV32Pmat267ZcyIiVUSj2tW48eRWfL54C9Mysr2OI0HgT6E7D3jJOfeic67g4EHnXKFz7iXgFWB4oANK5Rs/dxNmcF4PrYgSEakqrh7YkrQ61XhwwlIKi4qP/gYJK/4UumRgaRnnl/heI2GsuNgxbm4mx7euR6Pa1byOIyIiAZIQG809Z3Vg+ZY9vPfTRq/jSID5U+jWAWeUcX6w7zUSxg7OntNiCBGRqmdw54b0b5nM01+uYGduvtdxJID8KXSvA8PM7B0z62ZmCb5HdzN7CzgbeDU4MaWyfDhnEzXjYzijU0Ovo4iISICZGfcN7cju/QU8O3mV13EkgPwpdE8CY4DfAnMpWeW6D5gDXAyMcc49FfCEUmkOzp47W7PnRESqrA6NanFR32a8PWM9K7bs8TqOBIg/g4Wdc+56oAtwNyXlbozv6y7OuRuCE1Eqi2bPiYhEhj+e3pYa8TE8+OkS7fNaRZRrGLCZJQLfA/9wzr1MyQIIqWK+Wb6N1KRqpDfV7DkRkaqsTvU4/nh6W+7/ZAlfLt2q22yqgHJdoXPO5QItAa1zrqKKih3T1+QwoFWyZs+JiESAi/o2pW2DGjwycRl5BUVex5EK8uceum+Bk4KUQzy27Ofd7NpfwHGt63kdRUREKkFMdBT3D+3Ehu25vPbDWq/jSAX5U+j+AHQzs2fMrK32bq1aDk4O799KowRFRCLFca3rMahjA/7+7Wq27s7zOo5UgD+Fbi3QlpJitww4YGb5hz0OBCWlBN20jBxapVSnQa0Er6OIiEglundIRwqLHE98sdzrKFIB/lxl+xegpTBVUEFRMbPWbmdEula3iohEmqbJiVw9sAUvfpfBJf2a0aNpHa8jyTEod6Fzzl0exBy/YmYjgdFAB6CPc2627/jpwONAHJAP/Mk5943v3IWUjFBxwGbgYuecdiAuh4WbdpKbX8QAfdwqIhKRfndyaz6cs4nRE5by0Q0DiIrS4rhw489HrpVpMTAcmHrY8WxgqHOuC3AZ8DaA736+54CTnXNdgYXATZUXN7z9uDoHgH4tVehERCJRjfgY/nxmexZs3Mn4eZlex5Fj4FehM7MkM3vEzBaY2S7fY4HvWMCu0TrnljnnVhzh+Dzn3Gbf0yVAgpnFA+Z7VLeSmRu1KLlKJ+UwLSObjo1qUad6nNdRRETEI+d2T6V7kySe+GI5ew8Ueh1H/FTuQmdmrSi58nUXEA1MBr72fX0XsND3msoyApjnnDvgnCsAbgAWUVLkOgKvVWKWsJVXUMTc9Ts5rrWuzomIRLKoKGP0sE5k7TnAC9+s9jqO+MmfK3TPA3WAQc65zs65Ec654c65zsAZQJLvNeViZpPNbPERHueU472dgCeA63zPYykpdD2Axvy3eJb2/mvNbLaZzc7Kyipv5Cppzvod5BcVM6CV5s+JiES67k2SGJGexus/rGVd9j6v44gf/Cl0JwLPOucmH37COfcVJWXuxPJ+M+fcab5iePjj47LeZ2ZpwEfApc65DN/h7r7vmeFKNqX7NzCgjF97jHOul3OuV0pKSnkjV0nTMrKJjjJ6t6jrdRQREQkBdw5uR2y08fDEZV5HET/4U+j2AmVdztrqe03QmFkSMBG4yzn34yGnMoGOZnawnZ1Oyaw8OYppGTl0S6tNjXjNiRYREahfK4GbTmnD5GVbmboysj/FCif+FLr3gIvN7H/unDezBOAS4N1AhDKz88xsE9AfmGhmk3ynbgJaA38xs/m+R33fQokHgKlmtpCSK3aPBiJLVbYnr4CFm3bp41YREfmVK49vTrPkRB78dCkFRdrGPRz4c1lmAiV7uc43s5eBlZTMfGsPXAscACaY2a8+6nTOTfM3lHPuI0o+Vj38+MPAw6W852XgZX9/rUg2a+12ioqd5s+JiMivxMdEc++Qjlzz1mzembGeK45r4XUkOQp/Ct2h9849y393jbBSXmO+10QfWzQJtmkZOcTFRJHeTFPBRUTk107rUJ+BberxzFcrGdatMck14r2OJGXwp9BdEbQU4olpGTn0alaHhFh1bhER+TUz476zOzL4ue/5v69W8sh5XbyOJGXwZ+uvN4MZRCrX9n35LPt5N7cPaut1FBERCVFtGtTk0v7NeHPaOi7q24yOjWt5HUlKEapbf0mQzVhTst1Xfy2IEBGRMtxyaltqV4vlgQlLKJkMJqFIhS5CTcvIpnpcNF3TansdRUREQljtxFhuG9SOmWu38/niLV7HkVKo0EWoaatz6NOiLrHR+r+AiIiU7cI+TWnfsCaPTFxGXkGR13HkCPTTPAL9vGs/a7L3cVxrfdwqIiJHF+3b5zVz537GTF3jdRw5AhW6CDQ94+D9c5o/JyIi5dOvZTJDujTixe9Ws3nnfq/jyGFU6CLQtIwckhJj6dBQq5VERKT87jqrPc7B458v9zqKHKbchc7MOpnZ8MOOnWxmX5vZHDO7LfDxJNCcc0zPyKF/y2SiouzobxAREfFJq5PIdSe24pMFm/lp3Xav48gh/LlC91fg6oNPzCwV+AToClQD/mpmlwY2ngTahu25ZO7cr+2+RETkmFx/Yksa1U5g9CdLKCrWGJNQ4U+hSwemHPL8Ikq29erunOsIfA7cGMBsEgQ/rtb8OREROXaJcTHcdVYHlmzezQezN3odR3z8KXR1gK2HPD8TmOKcy/Q9nwBo24EQNy0jmwa14mmVUt3rKCIiEqaGdm1E7+Z1eHLSCnbnFXgdR/Cv0GUDaQBmVh3oD0w+5HwM/u0NK5Xs4P1zA1rVw0z3z4mIyLExM+4f2ontufk8P3mV13EE/wrd98ANvoURzwKxlNxDd1A7IPNIb5TQsHLrXnL25WtciYiIVFjn1Npc0KsJb0xbR0bWXq/jRDx/Ct3dwH7gQ+Aq4K/OuVUAZhYNnM+v77GTEDMtIxtACyJERCQgbj+jHdVio3no06VeR4l45S50zrm1QHugO9DCOXfXIacTgRuAxwIbTwJpWkYOTesmklYn0esoIiJSBdSrEc/Np7XhuxVZfLt8m9dxIppfg4Wdc4XOuYXOufWHHd/jnPvYObcuoOkkYAqLipmxJkdX50REJKAu7d+clinVeejTpeQXFnsdJ2L5VejMrK6ZPWRmP5rZKjPr7zuebGb3mVn74MSUilqyeTd78gp1/5yIiARUXEwUfzm7I2uy9/HmtHVex4lY/uwU0QSYD9wB1ARaUjJQGOdcDnAhmkMXsqb59m8doPlzIiISYCe3q8/J7VJ4/utVZO054HWciOTvThEJlNxDdwpw+NyLj33HJQRNy8imbYMapNSM9zqKiIhUQX85uyP7C4p4atIKr6NEJH8K3enA8865ZcCR9vpYCzQJSCoJqPzCYn5at11X50REJGhaptTgiuOa8+85G1m0aZfXcSKOP4WuOlDWEpYaFcwiQTJ/407yCop1/5yIiATV709tQ3L1OB6YsATntM9rZfKn0K0A+pVx/ixgccXiSDBMy8jGDPq1UKETEZHgqZUQy5/OaMfs9Tv4ZMFmr+NEFH8K3SvAxWZ2JRDtO+bMrKaZPQOcBLwY4HwSANNW59C5cW1qJ8Z6HUVERKq4kT2b0CW1No99tpzc/EKv40QMfwYLvwSMAV4F1vgOfwjsAG4G/uaceyfgCaVCcvMLmbdxBwNa6+qciIgEX1SUcf/QjmzZncfL32V4HSdi+DtY+CbgOEpK3efALOAlYKBz7pbAx5OKmr1uBwVFTgsiRESk0vRqXpdh3RrzytQ1bNud53WciOBXoQNwzk13zt3inBvinDvTOfd759yPwQgnFTctI4eYKKN38zpeRxERkQhy26C2FBY7XtRVukrhd6GT8DI9I5seTZNIjIvxOoqIiESQZsnVGdkzjXdnbmDzzv1ex6ny/N3662ozm2ZmW8zsgJnlH/bQeOgQsmt/AYsyd9FfH7eKiIgHbjqlNQ7HC9+u9jpKlVfuyzZm9hRwK7AZmA7sDFYoCYyZa3IodjBA8+dERMQDaXUSGdW7Ke/N2sANJ7aiSd1EryNVWf58DnclJQshznHOFQUpjwTQtIwcEmKj6NE0yesoIiISoW48uTVjZ2/k+a9X8eTIbl7HqbL8+cjVgAkqc+FjekYOvZvXJT4m+ugvFhERCYKGtRO4uG8zxs/LZG32Pq/jVFn+FLpJQO9gBZHAytpzgBVb92i7LxER8dwNJ7UiLjqK5yav9DpKleVPofsD0NvM7jez1GAFksCYsSYHQPPnRETEcyk147l0QDM+XrCZVVv3eB2nSvJnp4htwNvAfcAGMyvQKtfQNS0jh5rxMXRuXMvrKCIiIlx3QisSY6N5dvIqr6NUSf6scn0QuIeSVa6z0SrXkDYtI5u+LesSE61RgyIi4r261eO48vgW/O2b1dy4eTcddcEhoPxZ5XodWuUaFjbtyGV9Ti6X9W/udRQREZFfXH18S96Yto5nJq/kH5f28jpOleLP5ZsEtMo1LEzP8N0/11oLIkREJHTUTozlmoEt+f/27jzOrrq84/jnOzPZF7KRhJB9kyUIQoAkILKKUqmKuwKCWkCRgmBtrdZqW63WiuC+ICpo3SqICigJkXUSICwhREjITQIJgZCZ7Htm5ukf50x6GWa5s91z7+T7fr3u6849231+59yZ+8xvO3P/up4n17qhryu1J6G7C49yLQsLcrUMH9Cb6SMHZR2KmZnZK1x80kSG9O/FtXM94rUrtSeh+zhwnKQveJRr6YoIqnO1zJoynIoKZR2OmZnZKwzq24tLT5nCPcs28Ohzm7IOp8doT0K3FpgBfBaPci1Zq2p28NLW3b7dl5mZlawPzpnALLzPwwAAIABJREFUiIG9uXbusqxD6THaMyji50B0VyDWNR7Mef45MzMrbf17V3HZG6bwH7c/zcKVtcya7EqIzio4oYuIi7oxDusiC3I1HHJQXyYO9w2QzcysdJ0/awI/vH8l1961nF9dOgvJ3YQ6w5OU9SANDcGCXC1zpozwL4aZmZW0vr0qufy0qTy8eiMPrKjJOpyy12INnaRTACLivvzXbWnc3orvmZe2sWnnPvefMzOzsvCe48fxvXtyfO2u5Zw81ZURndFak+s9QEjqFxF7G1+3sr3S9ZVdFp21S3Uu+Q9nthM6MzMrA32qKrnijGl8+pYl/GXZy5x+2KisQypbrSV0pwGkydz+11a6FuRqmTRiAGOG9Ms6FDMzs4K887ixfPeeHNfOXc5prxnpWroOajGhi4h7W3ttpaWuvoGHVm3kb48Zk3UoZmZmBetVWcHfnzGNT/5mMX9eup43zRiddUhlqeBBEZLmSzqjlfWnSZrfNWFZez35wha276lz/zkzMys7bztmDJNHDODrc5fT0OAZ0jqiPaNcTwVaa9weCbyhU9FYhzXev3W25/IxM7MyU1VZwZVnTmPZ+m3cvuTFrMMpS+2dtqS1tHkKsL0TsVgnVOdqOGz0IIYP7JN1KGZmZu127mvHMH3UQK6bt5x619K1W6sTC0u6ALggb9GnJV3czKZDgNcBc7swNivQ7n31LFq9iQ+cOCHrUMzMzDqkokJ84szpfPTnj3HbEy9w3rFjsw6prLR1p4hhwLT05wBGA4OabBPADpJbg322S6Ozgjz+/Gb21DW4/5yZmZW1s48czRGHDOb6u5/l3KPH0KvS9z8oVKtnKiKuj4hJETGJZJ65qxpf5z0mR8RREXFRRKwtTtiWb0GuhgrBCZOHZR2KmZlZh1VUiKvPms5ztTu55TGnFO1RcOobERUR8T/dGYx1zIO5Wo4aO4TBfXtlHYqZmVmnnHH4SI4eN4Rv3L2CvXUNWYdTNlyXWea276lj8ZrNnOTmVjMz6wGkpJbuhc27+NWiNVmHUzac0JW5R1ZvpK4hmDNlRNahmJmZdYlTpo1g5oShfHv+Cnbvq886nLLghK7MLcjV0ruyguMmDM06FDMzsy4hiavfOJ2Xtu7mfx56PutwyoITujJXnavhdeOH0K93ZdahmJmZdZk5U0Ywe/JwvnNPjl17XUvXFid0ZWzzzr0sXbfVza1mZtYjXfPG6dRs38NNC1ZnHUrJc0JXxhaurCUC5kz1gAgzM+t5Zk4cxinTD+Z79+bYvqcu63BKWkEJnaQhkk6QNKWVbSZJurDrQrO2VOdq6d+7kqPHDsk6FDMzs25x9VnT2bRzHz95cFXWoZS0NhM6SV8E1gMLgOWSHpZ0bDObzgF+3MXxWSuqc7UcP3EYvatc0WpmZj3TMeOGcObhI/nBfSvZsmtf1uGUrFYzAUnvBj4N3AtcDnwJGA9USzq/+8Ozlry8dTcrXt7u232ZmVmP94mzprN1dx0/esC1dC1pq2rnKuAvEfHGiPheRPwLcDhwH/ATSVd1e4TWrAUrawE8IMLMzHq8I8ccxJtnjObGB1axacferMMpSW0ldIcBv81fEBGbgDcBNwFfk/Tv3RSbtaJ6RS2D+1ZxxJjBWYdiZmbW7T5x1nR27K3jB/evzDqUktRWQtfsTdQioiEiPgRcB3xG0rcKOJZ1oeqVNcyaPJzKCmUdipmZWbebPmoQ5752DD95cDU12/dkHU7JaSsJWw6c3NLKiLgG+ALwMcA1dUWyZuNO1mzc5f5zZmZ2QLnyzGnsqavne/fksg6l5LSV0N0JvFVSi5lDRHwB+AQwrisDs5ZV52oAOGmq+8+ZmdmBY8rBA3n768Zy88LnWL91d9bhlJS2ErobgU8BI1vbKCKuB94B/FtXBCXpXZKWSmqQNDNv+VmSHpW0JH0+PW/deyQ9me73X10RR6mqztUyYmAfpo4cmHUoZmZmRXXlGdOobwi+85cVWYdSUlpN6CLihYj4dkQ83daBIuJ3aW1dV3gKOI9kNG2+GuDciDgK+CBwM0Bag/hV4IyIOBIYJemMLoqlpEQE1bla5kwZjuT+c2ZmdmAZP7w/75o5ll88vIYXNu/KOpyS0eGBDJL6SrpQ0qiuDAggIp6OiGXNLH88ItalL5cCfSX1ASYDyyNiQ7puHkmNYY+T27CdDdv2uP+cmZkdsD5++jQAvjXftXSNOjMy9SCSO0Mc2UWxtNc7gMcjYg+wAjhM0kRJVcDb6KF9+qpznn/OzMwObIcO6cd7TxjHbxat4fnanVmHUxI6O9VIh9v8JM2T9FQzj7cWsO+RwFeAS2H/3HgfBX4F3A+sBlq8i6+kSyQtkrRow4YNLW1Wkh5cUcOhQ/oxbli/rEMxMzPLzOWnTaWyQnxj/rNZh1ISOpvQRYd3jDgzImY087ittf0kjQVuBS6MiP3jliPiDxFxYkTMBpYBLV7hiPhBRMyMiJkHH3xwR4tQdPUNwcKVGzlpqvvPmZnZgW3U4L6cP2sCtzy2lpUbtmcdTuYyq6Hr0JtJQ4DbgU9HxINN1o1Mn4eSzIt3QzFjK4anX9zKll373NxqZmYGfPTUKfSpquT6u11L15mEbgMwCXiwrQ3bS9LbJa0FZgO3S/pzuurjwFTgXyQ9kT4ap1S5XtJf03i+HBHLuzqurDXOPzfbAyLMzMwYMbAPH5wzkd8vXsfy9duyDidTHU7o0tt/PZcOSuhSEXFrRIyNiD4RMSoizk6X/0dEDIiIY/IeL6fr3hcRR6SPX3Z1TKWgOlfLlIMHMGpw36xDMTMzKwmXnjKZAb2ruG5ej6vHaRfff7VM7Ktv4OFVG93camZmlmfogN586KSJ3LHkJZau25J1OJlxQlcmFq/ZzM699Z5/zszMrIkPv34yg/tW8fW5B25fOid0ZaI6V4sEsyY7oTMzM8t3UL9e/N3rJzPv6fUsXrM563Ay4YSuTFTnajjikMEMHdA761DMzMxKzsUnT2Jo/15cO/fA7EvnhK4M7N5Xz2PPbXZzq5mZWQsG9qni0jdM4d7lG3j0uY1Zh1N0TujKwKPPbWJvfYMHRJiZmbXiwtkTGDGwN1+768CrpXNCVwaqczVUVojjJw3LOhQzM7OS1b93FR89dSrVuVoWpPc+P1A4oSsD1blajh57EAP7VGUdipmZWUn7wInjGTW4D9fOXUZEh+9QWnac0JW4bbv38eTaLZw01c2tZmZmbenbq5KPnzaVR1Zv4v5na7IOp2ic0JW4h1dtpL4hfLsvMzOzAr37+HEcOqQfX5u7/ICppXNCV+Kqc7X0rqrg2PFDsw7FzMysLPSpquSK06eyeM1m5j/zctbhFIUTuhJXnatl5oSh9O1VmXUoZmZmZeMdx41l/LD+XHuA1NI5oSthG3fs5ekXt3r+OTMzs3bqVVnBlWdMY+m6rfx56UtZh9PtnNCVsIUrkyHXsz3/nJmZWbu97XWHMvngAXx97rM0NPTsWjondCXswRU1DOhdyWvHHpR1KGZmZmWnskJcdeZ0lq3fxh+XvJh1ON3KCV0JW5Cr5cTJw+lV6ctkZmbWEW856hBeM2oQ181bTl19Q9bhdBtnCiXqxS27WFmzw/3nzMzMOqGiQnzirGms3LCD255Yl3U43cYJXYlqvGWJ558zMzPrnLOPHM2RYwZz/d3Psq+H1tI5oStR1blahvTvxeGjB2cdipmZWVmTxNVnTef5jTv57aNrsw6nWzihK0ERwYJcLbMnD6eiQlmHY2ZmVvZOP2wkx4wbwjfnr2BPXX3W4XQ5J3Ql6LnanbyweZf7z5mZmXWRxlq6Fzbv4tePrMk6nC7nhK4EVaf95+ZM9fxzZmZmXeX100Zw/MShfOsvK9i9r2fV0jmhK0HVuRpGDe7D5BEDsg7FzMysx0hq6V7D+q17+PlDz2cdTpdyQldiGvvPzZkyAsn958zMzLrS7CnDmTNlON+9ZwU799ZlHU6XcUJXYpav307tjr2ersTMzKybXPPG6dRs38tNC57LOpQu44SuxFTnagA8IMLMzKybHDdhGG+YfjDfvzfH9j09o5bOCV2Jqc7VMn5Yf8YO7Z91KGZmZj3W1WdNZ9POffz4gVVZh9IlnNCVkLr6BhaurOWkqa6dMzMz605HjxvCmYeP4of3r2TLrn1Zh9NpTuhKyNJ1W9m2u47ZUzxdiZmZWXe7+qzpbN1dx4/uX5l1KJ3mhK6ENM4/N3uya+jMzMy62xFjBnPOUaO58cHVbNqxN+twOsUJXQmpztUwfdRADh7UJ+tQzMzMDghXnTmdHXvr+P595V1L54SuROyta+CR1RuZ4+ZWMzOzopk+ahB/e/QYflq9mg3b9mQdToc5oSsRT6zZzO59DZ5/zszMrMiuPGMae+rq+d69uaxD6TAndCXiwRU1VAhmuf+cmZlZUU0+eCDnHTuWny18jvVbd2cdToc4oSsRC3K1zDj0IA7q1yvrUMzMzA44V54xjfqG4Nt/WZF1KB3ihK4E7Nxbx+NrNrm51czMLCPjhvXnXTPH8cuH1/DC5l1Zh9NuTuhKwKLVm9hXHx4QYWZmlqErTp8KwLfmP5txJO3nhK4EVOdqqaoQx08cmnUoZmZmB6wxQ/rxvhPG8ZtFa3m+dmfW4bSLE7oSsCBXw+vGD6F/76qsQzEzMzugXX7aVCorxPV3l1ctnRO6jG3ZtY8lL2xxc6uZmVkJGDm4LxfMmsCtj68lt2F71uEUzAldxh5aWUtDwBwPiDAzMysJl506hb69Krl+XvnU0jmhy1h1rpa+vSo4ZvyQrEMxMzMzYMTAPnxwzkT+8OQ6lr20LetwCuKELmMLcrUcP3EYfaoqsw7FzMzMUpe8fjIDeldx3bzlWYdSECd0GdqwbQ/L1m/z/HNmZmYlZuiA3nzo5Enc+dRLLF23Jetw2uSELkMLV9YCeECEmZlZCfrwyZMY3LeKr88t/b50TugyVJ2rZVCfKmaMGZx1KGZmZtbEQf16cckpk5n39HoWr9mcdTitckKXoepcDSdOHk5VpS+DmZlZKbropEkM7d+La+eWdl86ZxIZWbtpJ8/V7vR0JWZmZiVsYJ8qLnvDFO5dvoFHn9uYdTgtckKXkQW5tP/cVCd0ZmZmpezC2RMZMbAPX7urdGvpnNBlZEGuluEDejN95KCsQzEzM7NW9OtdycdOnUJ1rnZ/hUypcUKXgYigOlfLrCnDqahQ1uGYmZlZG95/4nhGD+7LtXOXERFZh/MqTugysKpmBy9t3e3+c2ZmZmWib69KLj99Ko+s3sT9z9ZkHc6rOKHLwINpde1Jnn/OzMysbLxn5jgOHdKPr81dXnK1dE7oMrAgV8OYg/oyYXj/rEMxMzOzAvWuquDvz5jK4jWbmf/My1mH8wpO6IqsoSFYkKtl9pQRSO4/Z2ZmVk7OO3YsE4b359oSq6VzQldkz7y0jU0797n/nJmZWRnqVVnBlWdMY+m6rfx56UtZh7OfE7oiq84lHSlnO6EzMzMrS2895lCmHDyAr899loaG0qilc0JXZAtytUwaMYAxQ/plHYqZmZl1QGWFuOrM6Sxbv40/Lnkx63AAJ3RFVVffwEOrNrq51czMrMz9zVGHcNjoQVw3bzl19Q1Zh+OErpiefGEL2/fUMcfTlZiZmZW1irSWbuWGHdz2xLqsw3FCV0yNtwuZNXlYxpGYmZlZZ5195ChmHDqY6+9+ln0Z19I5oSui6lwNh40exPCBfbIOxczMzDpJElefNZ3nN+7kt4+uzTQWJ3RFsntfPYtWb3Jzq5mZWQ9y2mtGcsy4IXxz/gr21NVnFocTuiJ5/PnN7Klr8IAIMzOzHkQS17xxOsMH9mbDtj2ZxVGV2TsfYKpzNVRWiBPdf87MzKxHOXnqCE6emu0doJzQFUl1rpajDj2IQX17ZR2KmZmZdaFSuJWnm1yLYPueOhav2ezmVjMzM+sWTuiK4JHVG6lrCA+IMDMzs27hhK4IFuRq6V1ZwXEThmYdipmZmfVAJZnQSXqXpKWSGiTNzFt+gqQn0sdiSW/PW3ecpCWSVkj6hkqhQTtVnavhdeOH0K93ZdahmJmZWQ9Ukgkd8BRwHnBfM8tnRsQxwJuA70tqHNjxXeASYFr6eFORYm3V5p17Wbpuq5tbzczMrNuUZEIXEU9HxLJmlu+MiLr0ZV8gACQdAgyOiAUREcBNwNuKFnArFq6sJQJOmuoBEWZmZtY9SjKha42kEyUtBZYAl6UJ3qFA/j031qbLMledq6V/70peO3ZI1qGYmZlZD5XZPHSS5gGjm1n1mYi4raX9IuIh4EhJhwM/lXQn0Fx/uWjlvS8haZ5l/Pjx7Yq7vfr1quTsI0fTu6rscmczMzMrE5kldBFxZif3f1rSDmAGSY3c2LzVY4F1rez7A+AHADNnzmwx8esKnz7n8O48vJmZmVl5NblKmtQ4CELSBOA1wOqIeBHYJmlWOrr1QqDFWj4zMzOznqQkEzpJb5e0FpgN3C7pz+mqk4HFkp4AbgU+FhE16bqPAjcAK4AccGeRwzYzMzPLhJJBoQeumTNnxqJFi7IOw8zMzKxNkh6NiJlNl5dkDZ2ZmZmZFc4JnZmZmVmZc0JnZmZmVuac0JmZmZmVOSd0ZmZmZmXOCZ2ZmZlZmXNCZ2ZmZlbmnNCZmZmZlTkndGZmZmZlzgmdmZmZWZlzQmdmZmZW5pzQmZmZmZU5J3RmZmZmZc4JnZmZmVmZc0JnZmZmVuYUEVnHkClJG4DnuuBQI4CaLjhOuXL5XX6X/8Dl8rv8Ln/xTIiIg5suPOATuq4iaVFEzMw6jqy4/C6/y+/yZx1HVlx+l78Uyu8mVzMzM7My54TOzMzMrMw5oes6P8g6gIy5/Ac2l//A5vIf2Fz+EuA+dGZmZmZlzjV0ZmZmZmXOCV07SXqTpGWSVkj6p1a2O15SvaR3FjO+7tZW+SWdKmmLpCfSx+eyiLO7FHL903PwhKSlku4tdozdqYDr/w951/6p9HdgWBaxdocCyn+QpD9IWpxe/4uziLO7FFD+oZJulfSkpIclzcgizu4g6UZJL0t6qoX1kvSN9Nw8KenYYsfYnQoo/2GSFkjaI+mTxY6vuxVQ/g+k1/1JSdWSji52jESEHwU+gEogB0wGegOLgSNa2G4+cAfwzqzjLmb5gVOBP2Yda4blHwL8FRifvh6ZddzFLH+T7c8F5mcdd5Gv/z8DX0l/PhjYCPTOOvYilv+rwL+mPx8G3J113F1Y/lOAY4GnWlh/DnAnIGAW8FDWMRe5/COB44EvAp/MOt4Myj8HGJr+/OYsrr9r6NrnBGBFRKyMiL3AL4G3NrPdFcBvgZeLGVwRFFr+nqqQ8r8fuCUingeIiJ70GWjv9X8f8IuiRFYchZQ/gEGSBAwkSejqihtmtymk/EcAdwNExDPAREmjihtm94iI+0iuZ0veCtwUiYXAEEmHFCe67tdW+SPi5Yh4BNhXvKiKp4DyV0fEpvTlQmBsUQLL44SufQ4F1uS9Xpsu20/SocDbge8VMa5iabP8qdlpk9Odko4sTmhFUUj5pwNDJd0j6VFJFxYtuu5X6PVHUn/gTST/2PQUhZT/W8DhwDpgCXBlRDQUJ7xuV0j5FwPnAUg6AZhABl9sGSn498N6vA+T1NYWVVWx37DMqZllTYcJXwf8Y0TUJ/+k9yiFlP8xktuSbJd0DvA7YFq3R1YchZS/CjgOOAPoByyQtDAilnd3cEVQSPkbnQs8GBGt1WiUm0LKfzbwBHA6MAWYK+n+iNja3cEVQSHl/zJwvaQnSBLax+k5NZRtac/vh/VQkk4jSehOLvZ7O6Frn7XAuLzXY0n+E883E/hlmsyNAM6RVBcRvytOiN2qzfLnf3FFxB2SviNpRET0hPv8FXL91wI1EbED2CHpPuBooCckdIWUv9F76VnNrVBY+S8GvhxJR5oVklaR9CV7uDghdqtCf/8vhmSQALAqfRwI2vP7YT2QpNcCNwBvjojaYr+/m1zb5xFgmqRJknqTfGn9Pn+DiJgUERMjYiLwv8DHekgyBwWUX9Lo9A95Y5NLBVD0D3Y3abP8wG3A6yVVpc2OJwJPFznO7lJI+ZF0EPAGknPRkxRS/udJamdJ+469BlhZ1Ci7TyG//0PSdQAfAe7rIbWThfg9cGE62nUWsCUiXsw6KCsOSeOBW4ALsmqRcQ1dO0REnaSPA38mGfF1Y0QslXRZur4n9pvbr8DyvxP4qKQ6YBfw3rS2ouwVUv6IeFrSn4AngQbghohodph7uWnH5//twF1pLWWPUWD5/x34iaQlJE1w/9hDaqcLLf/hwE2S6klGe384s4C7mKRfkIziHyFpLfCvQC/YX/Y7SEa6rgB2ktZU9hRtlV/SaGARMBhokHQVySjoHpHQF3D9PwcMB76T1mnURcTMosbYQ75rzczMzA5YbnI1MzMzK3NO6MzMzMzKnBM6MzMzszLnhM7MzMyszDmhMzMzM+skSTdKellSl8xsIKle0hPp41VTRDXlhM7MWpTewuyerOMoNknvlrRU0h5JIWlI1jG1RtLENM6Lso6lPSRdlMY9MetYzLrAT0huedhVdkXEMenjb9va2AmdWZmQdIukfZIObmWbK9MvyHOLGVtPImkq8HPgJeCjwAVAj5pTz8y6XkTcB7zidoeSpkj6U3pv7/slHdZd7++Ezqx83EwyGfh7W9nmfKAG+FMXvecb08eB5FSS83xNRNwYET+LiH0Zx2Rm5ekHwBURcRzwSeA77di3r6RFkhZKeltbG/tOEWbl43aS//7OB77ZdKWk6ST3Ev5WZxMQSf0jYmdE7O3MccrUyPR5c0d2ltQP2BMRDV0XkpmVG0kDgTnAb9K7RwD0SdedB/xbM7u9EBFnpz+Pj4h1kiYD8yUtiYhcS+/nGjqzMpEmV78GTpA0rZlNLkifb5bUW9IXJD0saaOkXWnH2oua7iRptaR5kk6RVC1pF/CldN2r+tBJuiZtOtiQ9jF7RtInlfcXK2/fFZKmSvqzpB1ph+EvS3rV3x5J75D0gKRtkram/5l+uMk2x0r6vaRNaZkWFfKfa7pvhaRPSVqWxr1O0rfz+8dJWg18MX25Km2+/kkrx2zsA3ampGslrSNpnh0saZik/5K0OC3PjvT8vqWZ44SkGyS9WdLjknan5+79zWx7iKTfpOdpo6Qfkdxuqbn4TpJ0d7rt9vTn2S2U4XRJX5H0Urr9b5Tcm7VK0pfS87UrPf8tNvvnHXdAerxcWp7atKbhnQXse2667U5JmyXdJunwJtt8Po17hqSfpp+JrZJ+IWlkM8fs8GfHrIMqgM15/eCOiYjDASLiloiY0cyjMZkjItalzyuBe4DXtfVmZlY+bk6fP9DMuvcDyyPiYZIv+MuAhcC/AJ8mqXH6saS/a2bfScBtwALg74G/tBLD1cDTwH+kPz8DfDV93dQgYB7JDeqvAR4E/hF4RQyS/gn4X6A/STL5T8CjwLl527w+3f9QkqTrH0jumXmrpPe1Em+j7wBfAZancd9Kco7u1v/fUP6qNA6AT5Akyd8v4NjXAbPT438G2AtMBt4N3AV8iuTej/2A30tqrhn7eJJO1b8naZrZTpKc709kJPUF7gbeCtyQHnMycFPTg0k6BZifrv8SyTmbAvxF0knNvP9/AyeQXMefAu8AfgR8F5iVHuMG4C3A9QWck++QnMM/AB9P3385cGJrO6XX8jaSc/VZknN7MlCtpH9jUz8DxpJ8zn9Mcj/pu/KuaVd8dszaLb2P7SpJ7wJQ4uhC9pU0VFJjbd4I4CSS+yO3+oZ++OFHGT1Ibv79bJNlJwEBfDZ9XQn0aWbfuc3suzrd97xmtr8HuKfJsv7NbHcDSQLSp8m+AVzWZNsngEfyXk8C6kgSn15Ntm2837RIksj7gMr89cADwJrGbVs4ZzPSWH7RZPnl6fKP5S37bLpsYgHX4qJ020ebib1Pfqzpst7AUmBuk+WRnoMZectGAXuAr+YtuyLd9uK8ZZXpeQngorzli4BNwKi8ZYcAW4CHmynDfUBF3vJfAQ0kyX3T5fuAgW2cm03Atws8fxPT172AF4Fn848PvBaoB36dt+zz6b7zm3wm/i5dfmlXfHb88KPQB/CL9PO7D1gLfDj9+/YnYDFJQva5Ao81B1iS7rcE+HBb+7iGzqz8/AyYKmlW3rLzSb7Efg4QEfURsQdAUq+0+W8EyZffVEkHNTnmSyQ1Vm2KiJ3pcavS/yJHkCRvA4DXNNl8H0myl+9eklqjRueRJCWfjyZ9/yL9ywYcDRyWlm+opBHp+w4H7iCpoZneStiNzZxfbbL8hyQ1l69qBm2nHzYT+56IqAeQ1EfScJKa0/uA45o5xn0R8VTe/utJaj/zz9VbgFr+v6aW9D1e0adS0uj0PW5Oj9O47Yskn5/jJY1qpgz5/f4WkCQ9NzazvAoY30wZ8m0GTpQ0ro3t8h0HjAa+GxHb8+J+kuRL8c16dXP9NxvPc+onvPKadvazY1aQiHhfRBwSEb0iYmxE/CgiVkXEmyLi6Ig4IiKa6zfX3LGqI+KodL+jIuJHbe3jhM6s/PwsfT4fIG1aejfwQESsatxI0gclPQnsJkkCNpD2jQOaJnQr85KnVkk6R9JCYBfJII0N/H+C0XS+thcioq7Jsk3AsLzXjc1oS1p528Yv3O+l75f/aOzz9qp+U3kmps/P5C+MpF/iCpL/ojvjVR2V0+aVayQtJ7kGNWm8l/Hq8wTwXDPLmp6rCSTXquk5Xdbk9cT0+Rle7a9Ntmn0fJPXm9tYPrSZY+e7BjgceE5J/82vSmoukc3XGFNLcQ8Emvbfe0XZ08R6Vd6xOvvZMSsLHuVqVmYiYoWkBcB7JF0FnEPypb+/1kbSe0hqKm4HrgXWk9SWnUPSr6npP3O7CnlvSXNI+kQtAD4GvEDSX+xYkv5jTY9bT9vU9ib7j/vPwCMtbNPR2dlFUrvZGc2dv08BXyajBP3IAAAD8UlEQVS5Ll8gSejqgYtJ+js21dK5UpOfm4u1kHPYdNumx2np/QuJ61Ui4hZJD5D0gzwT+BBwjaTPRMR/FhpsM+/XNO62zkd3fnbMSoYTOrPydDNJp/OzSWrq9gC/yVv/XpJainPza94knd7J930XSQJ3ZkTszjvu5JZ3adOz6fNRQHUL26xIn3dExLwOvMfq9Pkw4LHGhZJ6kTRpLuzAMdvyXpL+hxfmL5T0oU4cczVwnKSqJrV0TZsMV6fPzU1i2risuRrBLhURL5MMrPiRpP4k/2B8QdJ/N22iTq3Oi/GOJusOI+mnWdPM8qcbX6TXdCJJ/zjo/GfHrCy4ydWsPP2KJLG6nKSv0B8iIn/etMY+T/t/x9M+XJ1JJhqP20DS563xuH1JOut31C0ktUBfSL+M95P2T4XyGEnid42auQ1XAdNo/DF9vrrJ8o+QNB3+ob1BF6CBJn9jlUw38/ZOHPN2kr5fjVPUIKmSJuc/Il4iGRRxQf4UHmnfugtIBkWsp5tIqmzaTzPte7mMZODDgBZ2XUTSn/MySfu3kTSD5JZKd8Sr5/e7Ij0HjS4iadK+PX3d2c+OWVlwDZ1ZGYqIjZLuABrn0bq5ySa3kQw2+KOk35H0EboEWEcyerKjfk/SZDtP0s0k05J8kKSPWIdExCpJ/0oyXcZDkn5NMhJzBsmozPMiokHSxSQjYf8q6UaS2pzRJNNgHEEyJUdL7/GUpO8Dl0oaTNLB/giS/myPkdQidbXbSJLUn5MMGplA0kz9DHBMB4/5w/QY35d0FMl0MO+g+XnoriEZ1bwwLbuAS4G+vDqx7WqDgBck3UoySm8jyRxaHwHubPLPx34RUSfpapIBDA9K+ilJ2a4AtpFMCdPUMJJpSm4l6Y95OUkT6o/TY3bqs2NWLpzQmZWvm0kSulrgzvwVEXFTWiN3OcmcYc+TzDO2hfSLriMi4l5JF5D0R7oWeJmkr979JF+YHT3uFyXlSOaB+xzJFB7LyLtNTkQ8KOkEkvnGLiGphVlPkjA090Xf1MdIEqCPkNT21JLclucz0T13xPhPkqlLLiBJupaTzMc2nQ4mdBGxS9IZJNf0EpJa2lvT14ubbHtfuu2/kZwzgIeBD0RES03bXWUn8C2SvnN/Q3IenicZlPNfre0YEb+QtIPkmn6JpIz3AJ+OiBXN7HI+SX/Ffyf5TrsFuLJxlHd6zM5+dsxKnqKwgW1mZmYlQ9LnSSZWHhcRazMOxyxz7kNnZmZmVuac0JmZmZmVOSd0ZmZmZmXOfejMzMzMypxr6MzMzMzKnBM6MzMzszLnhM7MzMyszDmhMzMzMytzTujMzMzMypwTOjMzM7My93/pKBOpaS6hIAAAAABJRU5ErkJggg==\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.4rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
