{
 "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').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": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2114: ConvergenceWarning: Retrying MixedLM optimization with lbfgs\n",
      "  ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2114: ConvergenceWarning: Retrying MixedLM optimization with cg\n",
      "  ConvergenceWarning)\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": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2118: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2114: ConvergenceWarning: Retrying MixedLM optimization with lbfgs\n",
      "  ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2114: ConvergenceWarning: Retrying MixedLM optimization with cg\n",
      "  ConvergenceWarning)\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": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2118: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/home/travis/build/statsmodels/statsmodels/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\").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": [
      "/home/travis/build/statsmodels/statsmodels/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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3iV5f3H8fc3mxEIkLB3CHsTQKYDERURcI866q6z1mptq1atttra+nMrWm3Vat2g4kRZgoBhb0LCngmBEAjZ9++PHGyKIeTASZ7knM/rup6LnOc55+ETLi/zzX0/3/s25xwiIiIiEvzCvA4gIiIiItVDhZ+IiIhIiFDhJyIiIhIiVPiJiIiIhAgVfiIiIiIhQoWfiIiISIiI8DpAbRAfH+/at2/vdQwRERGRY1q4cGGmcy6hvGsq/Cqhffv2pKSkeB1DRERE5JjMbNPRrmmqV0RERCREqPATERERCREq/ERERERChAo/ERERkRChwk9EREQkRKjwExEREQkRKvxEREREQoQKPxEREZEQocJPREREJESo8BMREREJESr8REREREKECj8RERGREKHCT0RERCREqPATERERCREq/ERERERChAo/ERERkRChwk9ExEPOOQqLS7yOISIhIsLrACIioSY7t5A5aZnMWpfB7NRMdu7Po0/rhgzrFM/QxHj6t4sjOiLc65giEoRU+ImIVLGi4hKWbNnHrNRMZqdmsHTLPkocxEZHMLRTE8b2bsGCDVk8N309z3y7nuiIMAa2b8zQTk0YlhhPz1YNCQ8zr78NEQkCKvxERKrAlqxcZqVmMGtdBnPX7yEnv4gwg96t47j1tCRGJsXTt00cEeH/feJmf14h89OzmJuWydz1e/jLF2uBtcTGRHBSxyYMS2zC0E7xJDWtj5kKQRHxnwo/EZEAOJBfxLy0PcxKLZ2+3ZB5EICWDWMY27sFI5ISGNapCXF1o456jwYxkYzu3ozR3ZsBkJGTz9y0TL5P28OctEy+XrULgITYaIYmlo4GDklsQpvGdav+GxSRoGDOOa8z1HjJyckuJSXF6xgiUoOUlDhWbt//46jeos17KSx21IkM56SOjRnZOYERSQkkJtQL2Ojclqxc5qZlMmf9Huam7SHzQD4AbRvXZVinJgz1FYLx9aMD8veJSO1kZgudc8nlXlPhd2wq/EQEYNf+vB8bMr5bn0nWwQIAurdowMjOCYxMimdA+0bV0pjhnCN19wHmrC8tBOenl04nA3RtHsvQxHiGJjZhcMfGxMZEVnkeEak5VPidIBV+IqEpr7CYBRuyfiz21u7KASC+fjQjk+IZ0Tme4Z0SSIj1foStqLiEFdv3M2d9JnPTMknZuJf8ohLCw4zerRv+ODXcv10jYiLVMSwSzFT4nSAVfiKhwTnHul0HmLUug1mpGSzYkEV+UQlR4WEM7NCIEUkJjExKoGvzWMJqeJdtXmExizbvZe760ucDl23NprjEER0RRnL7Rj+OCPZq1fB/GkxEpPZT4XeCVPiJBK+sgwXM9jVkzE7NYNf+0ufmOjWtz8ikBEZ0juekDk2oE1W7R8ly8gpZsCHL93xgJmt2lo5exkZHMLhjk9IRwU7xdG6mjmGR2q6iwk9dvSISUgqKSli0eS+zUzOYtS6TFduzcQ4a1olkeFJ86RRuUgIt4+p4HTWgYmMiGdWtGaO6lXYMZx7I5/u00iaRuWmZTFtd2jEcXz+KIYnxDPMVguoYFgkuGvGrBI34idRezjk27sn1PaeXwfdpezhYUEx4mNG/bVzp9G3nBHqF+CLJW/fmMtc3GjgnbQ8ZOaUjn20a12Fox3iG+rqGa8LzjCJSMU31niAVfiK1y/68Quauz2RWaum2aFv3HgJKi5iRvkJvSGITGqjbtVzOOdbvPsDctD3MWZ/JvPQ97M8r7Rju3Kw+QxPjGdYpnsEdG+vfUKQGUuF3glT4idRsxSWOpVv3MXtdJrNSM1iyZR/FJY760REMSWzy4/Rt+/h6XketlYpLHCu2Zf84LfzDxizyCksIM+jVOu7HaeEB6hgWqRFU+J0gFX4iNc/2fYd+7L79LjWT/XlFmEHvVg1/nL7t1zaOSHWsBlx+UTGLN+9j7vrSaeGlW/ZRVOKIigjjouTW3De2uwpAEQ+puUNEgoZzjkmz0nnsizU4B80aRDOmR3NGdk5gWKd4Gtc7+pZoEhjREeGc1LEJJ3Vswq8o3a7uhw1ZfLVqJ2/O28zCTft47rJ+dEyo73VUETmCRvwqQSN+IjVDYXEJD0xZwdsLtjC2VwvuOD2JpKZafqQmmb5mN3e+u4SiYsdj5/finN4tvY4kEnIqGvHTHIiI1Ar78wq55p8/8PaCLdxyaiLPXNqPzs1iVfTVMKd2bcrU20fQuVl9bn1rMfdPXkF+UbHXsUTER4WfiNR4W7JyOf/5uXyftoe/XNCbu8d0rfE7Z4SyVnF1eOfGIVw/ogNvzNvE+S/MZfOeXK9jiQgq/ESkhlu8eS8Tn5/Drv15vH7NIC5KbuN1JKmEyPAwfj+2O5OuGMDmPbmMfWY2X6zY6XUskZCnwk9EaqzPlu/gkknzqBMVzoc3D2Nop3ivI4mfzujRnKm3j6BjfD1uenMhD32ykoKiEq9jiYQsFX4iUuM453hxZho3/3sRPVo2YPLNw+jUVB2itVWbxnV596YhXD20Pa/N2ciFL33P1r2a+hXxggo/EalRCotL+O2Hy3ns8zWc07sFb11/Ek3qa5uw2i46IpwHz+3B85f3J333AcY+/R3TVu3yOpZIyKmRhZ+ZXWhmK82sxMySy5wfbWYLzWy578/TylybYWZrzWyJ72h6lHv/1szW+947pjq+HxGpnOxDhfz8tR/4zw9buPXUTjx9ST8tBBxkzu7Vgk9uG06ruDpc93oKf/psNYXFmvoVqS41dQHnFcB5wEtHnM8ExjnntptZT+BLoFWZ65c754664J6ZdQcuAXoALYFpZtbZOae1BkQ8tiUrl2v++QMb9xzkrxf05kI1cQSt9vH1+PDmoTwydRWTZqWzcNNenrm0Hy3j6ngdTSTo1cgRP+fcaufc2nLOL3bObfe9XAnEmJk/c0Djgf845/KdcxuA9cCgE08sIififzt3B6voCwExkeE8MqEXT1/ajzU79jP26dlMX7vb61giQa9GFn6VdD6w2DmXX+bca75p3vut/FVdWwFbyrzeyv+OGIpINTvcuVs3KoIPbx7GkMQmXkeSanRun5Z8fNtwmjWI4eev/cBfvlhDkaZ+RaqMZ4WfmU0zsxXlHOMr8dkewOPAjWVOX+6c6wWM8B1XlPfRcs6Vu2edmd1gZilmlpKRkXHsb0hE/OKc44UZpZ27PVs15KObh6pzN0QlJtRn8i3DuGRgG56fkcZlr8xn1/48r2OJBCXPCj/n3OnOuZ7lHFMq+pyZtQY+Aq50zqWVud823585wFuUP4W7FSg7h9Qa2F7O+3DOTXLOJTvnkhMSEvz75kSkQoXFJdz7wXIe/2IN4/q05N/XDVbnboiLiQznsfN78/eL+rB8azZnPzWb71IzvY4lEnRq1VSvmcUBU4HfOufmlDkfYWbxvq8jgXMobRA50sfAJWYWbWYdgCRgQdUnF5HDsg8VcvVrC3gnZQu3ndaJpy7uq85d+dF5/Vvz8a3DaFwviitenc+TX6+juKTciRkROQ41svAzs4lmthUYAkw1sy99l24FOgH3H7FsSzTwpZktA5YA24CXffc618weBnDOrQTeBVYBXwC3qKNXpPpsycrl/BfmsmBDFk9c2Ie7zuiiPXflJ5KaxTLl1mFM7NeKp75J5cpX55ORk3/sD4rIMZlz+k3qWJKTk11KylFXiRGRSli0eS/X/yuFohLHiz8boCYOOSbnHO+lbOX+KStoUCeSpy/pp/9uRCrBzBY655LLu1YjR/xEJLhMXbaDSyfNo35MBB/ePFQ/vKVSzIyLBrZh8i3DiI2O4PJX5vHst6mUaOpX5Lip8BORKuOc47np67nlrUX0atWQj24eRmKCOnfFP91aNODj24ZzTu+WPPHVOq7+5w/sOaCpX5HjocJPRKpEYXEJv/lgGX/9ci3n9mnJm9cNpnG9KK9jSS1VPzqCpy7py6MTezIvfQ9jn/6OHzZmeR1LpNZR4SciAZd9qJCrXl3Auylbuf20Tjx1iTp35cSZGZcPbseHvxhKdGQYl0yax4sz0zT1K+IHFX4iElCHO3d/2JjF3y7sw6/O6EL5G+mIHJ+erRryyW3DGdOjGY99vobrXk9h78ECr2OJ1Aoq/EQkYBZu2suE5+aQkZPPG9cO5vwBrb2OJEGqQUwkz13Wn4fO7cHs1AzGPj2bRZv3eh1LpMZT4SciAfHpsu1c+vJ/O3dP6qjOXalaZsZVQ9vzwS+GEhZmXPTi97wyOx0tUyZydCr8ROSEHO7cvfWtxfRprc5dqX69W8cx9bYRnNq1KY9MXc2NbywkO7fQ61giNZIKPxE5bgVFJdzzfmnn7vi+6twV7zSsG8mkKwZw39hufLtmN2Ofmc2yrfu8jiVS46jwE5Hjkp1buufuewu3cvuoJP7v4r5ER6hzV7xjZlw3oiPv3jSEkhLHBS98z7/mbtTUr0gZKvxExG+b9+Ry3gtz/tu5O7qzOnelxujfthFTbx/B8KR4/vDxSm59azH78zT1KwIq/ETETws37WXi83PIPFCgzl2psRrVi+KVK5O596yufLFyJ+c+8x0rt2d7HUvEcyr8RKTSPlla2rkbGxPBR+rclRouLMy46eRE/nPDSRwqLGbi83P59/xNmvqVkKbCT0SO6XDn7m1vl3bufnjzMDqqc1dqiYHtG/PZ7SMY3KExv/9oBb98ZwkH84u8jiXiCRV+IlKhgqIS7vZ17k5Q567UUk3qR/Ovnw/irtGd+WTpdsY9+x1rdu73OpZItVPhJyJHlZ1buufu+wu3cseoJJ5U567UYmFhxm2jknjzusHk5BUx4bk5vJuyxetYItVKhZ+IlOtw527Kpiz+flEf7lTnrgSJoYnxTL19OP3aNOKe95dx17tLyS3Q1K+EBhV+IvITCzdlMeH5Oew5WMCb1w7mvP7q3JXg0jQ2hjevG8zto5L4cPFWxj87h9RdOV7HEqlyKvxE5H98vHQ7l748nwYxEXx08zAGq3NXglR4mPGr0Z15/ZpBZB0s4Nxn5/Dhoq1exxKpUir8RAQo7dx99ttUbn97MX1bx/HRzcPoEF/P61giVW5EUgKf3TGCXq0b8qt3l/Kb95eRV1jsdSyRKqHCT0QoKCrh1+8t44mv1jGxXyveuG4QjdS5KyGkWYMY3rpuMDefksg7KVuY+Pxc9uUWeB1LJOBU+ImEuOzcQq58dT4fLNrKL09P4u8X9VHnroSkiPAw7jmzK69enUza7gPc8tYiCotLvI4lElAq/ERC2KY9B5n4whwWbdrHkxf34Zenq3NX5LSuzXh0Yk/mrN/DI5+u8jqOSEBFeB1ARLyRsjGLG95YSIlzvHndYAZ1aOx1JJEa48LkNqzblcPLszfQuXkslw9u53UkkYBQ4ScSgqYs2cbd7y+jVVwdXr16oJo4RMpx71ndSN19gD9MWUnH+PoMSVSHu9R+muoVCSHOOZ75JpU7/rOEvm3i+PAXQ1X0iRxFeJjx9KX9aNekLjf/eyGb9+R6HUnkhKnwEwkRxSWOu99fxt++9nXuXqvOXZFjaRATyStXDaTEwXWv/0BOXqHXkUROiAo/kRDxyux03l+4ldtHqXNXxB8d4uvx/OX9Scs4yJ3vLKG4xHkdSeS4qfATCQGpu3L429frGNOjGXeenqTOXRE/DesUzwPndGfa6t088dVar+OIHDc1d4gEuaLiEn793lLqR0fw6MReKvpEjtOVQ9qxdlcOL8xIo3Oz+kzspz2spfbRiJ9IkHtpVjpLt2bzx/E9ia8f7XUckVrLzHjo3B4M7tCY33ywnMWb93odScRvKvxEgtianfv5v2nrGNu7BWN7t/A6jkitFxkexgs/G0CzBtHc8MZCdmQf8jqSiF9U+IkEqcLiEu56dykN60Tyx/E9vY4jEjQa14vilSsHkptfxA2vL+RQQbHXkUQqTYWfSJB6fnoaK7fv55EJvWisZVtEAqpL81ieuqQfK7Znc88Hy3BOnb5SO6jwEwlCK7dn88y3qYzv25Izezb3Oo5IUDq9ezPuHtOFT5Zu57np672OI1Ip6uoVCTIFRaVTvI3qRfHguB5exxEJar84OZF1O3N44qt1dGoaq1+0pMbTiJ9IkHn221TW7MzhTxN7aWcOkSpmZjx2fm/6tInjV+8uYfWO/V5HEqmQCj+RILJ8azbPzUjjvP6tGN29mddxREJCTGQ4L18xgNiYCK77VwqZB/K9jiRyVCr8RIJEflExd723hPj6UfzhHE3xilSnpg1iePnKZDIP5POLNxdSUFTidSSRcqnwEwkST01LZd2uAzx2fm8a1o30Oo5IyOndOo4nLuzDDxv3cv/kFer0lRpJzR0iQWDx5r28ODONi5PbcGqXpl7HEQlZ4/q0ZN2uHJ75dj1dmsdyzfAOXkcS+R8a8ROp5fIKi/n1e0tp3iCG35/Tzes4IiHvztM7M6ZHMx6ZuoqZ6zK8jiPyP1T4idRyf/96HWkZB3n8gt40iNEUr4jXwsKMv1/Ul87NYrn1rUWkZRzwOpLIj1T4idRiCzdl8fLsdC4b3JYRSQlexxERn3rREbxyVTJR4WFc/68UsnMLvY4kAqjwE6m1DhUU8+v3ltGyYR1+d7ameEVqmtaN6vLiFQPYsjeXW99eRFGxOn3Feyr8RGqpv365lg2ZB/nrBb2pH60+LZGaaGD7xjwyoSezUzN59LPVXscRUVevSG00P30Pr83dwJVD2jG0U7zXcUSkAhcPbMvanQd4dc4GujSL5ZJBbb2OJCHsqIWfmZUAfi9C5JwLP6FEIlKh3IIi7n5/GW0a1eU3Z3b1Oo6IVMLvzu7K+owD3D9lBR0T6jOoQ2OvI0mIqmjE72F+WvhNAHoCXwJrAQO6AGcAy4EpVZBRRMp4/PM1bM7K5Z0bTqKepnhFaoWI8DCeubQfE5+bw01vLmTKLcNo07iu17EkBB31p4Zz7sGyr83sGqA50NM5t/aIa92A6cDmKsgoIj5z12fyr+838fNh7RncsYnXcUTEDw3rRPLKVclMeG4O17+ewvu/GKrnc6Xa+dPccQ/w7JFFH4BzbjXwHPCbQAUTkf91IL90irdDfD3uGaMpXpHaqGNCfZ69rD/rduVw5ztLKCnRtm5Svfwp/NoBhyq4nut7j4hUgT99tprt2Yd44sLe1InSo7QitdXIzgncN7Y7X6/axd+/Xud1HAkx/hR+64AbzCzuyAtm1gi4gdLn/kQkwGaty+Ct+Zu5fkRHBrTTQ+Eitd3Ph7XnkoFteHb6eqYs2eZ1HAkh/jxc8DtgMpBqZm9QWgg6oCvwMyCO0uYPEQmg/XmF3PvBMhIT6vGr0Z29jiMiAWBmPDy+J+kZB7nn/WW0b1KPPm1+Mq4iEnCVHvFzzk0FxlDawPFL4HngBeAO37mzfO8RkQB69NPV7NyfxxMX9iEmUlO8IsEiKiKMF37Wn/j60dzwRgq79ud5HUlCgF87dzjnvnXODQBaAEOAoUAL59wA59y0qggoEsqmr93NOylbuPHkRPq1beR1HBEJsCb1o3nlqmRy8oq44fUU8gqLvY4kQe64tmxzzu1yzs13zs1zzu0KdCgRgezc0inezs3q88vTk7yOIyJVpFuLBjx5cV+Wbs3mNx8swzl1+krV8avwM7M4M3vMzFaY2UEzO+D7+k/lNX2IyPF7+NNVZB4o4IkL+xAdoSlekWA2pkdzfn1GZ6Ys2c4LM9O8jiNBrNKFn5m1AhZTup4fwFTgc0obPO4FFplZy4AnFAlBX6/axQeLtnLzKYn0bq3fqURCwS2ndmJcn5b89cu1fL1Kk2lSNfwZ8fsz0Aw4xznX0zl3kXPuQudcL2Cs79qfqyKkSCjZe7CA3320nK7NY7ntNE3xioQKM+OvF/SmV6uG/PI/i1mzc7/XkSQI+VP4nQk85Zz77MgLzrnPgWeAswIVTCRUPfjJSvYeLOBvF/UhKuK4HsMVkVoqJjKcSVckUy86guv+lcKeA/leR5Ig489PlVigolUmt/reIyLH6YsVO5iyZDu3nZZEj5YNvY4jIh5o3jCGSVcmszsnn1/8exEFRSVeR5Ig4k/htxa4wMx+8hkzCwcuQDt3iBy3PQfy+f1HK+jRsgE3n5rodRwR8VDfNnH89YLeLNiQxR8+XqlOXwkYf3bueBp4BfjWzJ7kv0VeV0oXcR4BXBvYeCKh44GPV7I/r5B/Xz+YyHBN8YqEuvF9W7F2Zw7Pz0ija/NYrhra3utIEgQqXfg55141s6bAH4APy1wyIB/4nXPun4GNJxIaPl22nanLdnD3mC50bd7A6zgiUkP8+owurNt1gIc/XUViQn2GJ8V7HUlqOfN3+NjMmgCjgXa+UxuBr51zWYGNVnMkJye7lJQUr2NIkMrIyeeMJ2fSpnFdPvzFUCI02iciZRzIL+L85+eyI/sQU24dTof4el5HkhrOzBY655LLu+b3Txjn3B7n3H+cc4/7jncCXfSZ2YVmttLMSswsucz50Wa20MyW+/48rcy1GWa21syW+I6m5dy3iZlN9y08/WwgM4scD+cc901ezsH8Yv52YR8VfSLyE/WjI3jlqmTCw4xr//UD2YcKvY4ktZjfP2XM7Ewze9bMpprZp76vzwhwrhXAecCsI85nAuN8awdeBbxxxPXLnXN9fcfucu6bB9wP/DrAeUWOy8dLt/Plyl386ozOJDVTU7yIlK9N47q88LMBbN6Ty+1vL6a4RM0ecnz82bkjysymULpjx83AQGCw7+vPzWyymUUFIpRzbrVz7icdws65xc657b6XK4EYM4v2474HnXPfUVoAinhq9/48Hpiykn5t47h+REev44hIDXdSxyY8PL4nM9dl8OfPVnsdR2opf0b8/gCMA/4GNHXONXXOJQAJwBPAuZSOplWX84HFzrmyq1u+5pvmvd/MrBqziPjFOcfvPlpOXmExT1zYh/Aw/ecqIsd22eC2XD20Pa98t4F3f9jidRyphfwp/C4D3nTO3eOcyzx80vfM32+AN4GfVfZmZjbNzFaUc4yvxGd7AI8DN5Y5fblvCniE77iislmO8nfcYGYpZpaSkZFxIrcS+YkPF21j2urd3D2mC4kJ9b2OIyK1yH1juzG8Uzy/n7yclI1B21cpVcSfwq8lMLeC698DLSp7M+fc6b49f488plT0OTNrDXwEXOmcSytzv22+P3OAt4BBlc1ylHyTnHPJzrnkhISEE7mVyP/YmZ3Hg5+sZGD7Rvx8WAev44hILRMRHsazl/WjVVwdbnpzIVv35nodSWoRfwq/7cBJFVwfDOw4sTgVM7M4Sp8x/K1zbk6Z8xFmFu/7OhI4h9IGEZEaxTnHvR8uo7C4hL9eoCleETk+cXWjeOWqgeQXlXD96ws5mF/kdSSpJfwp/N4GrjCzR8ys0eGTZtbIzP5I6dTqW4EIZWYTzWwrMASYamZf+i7dCnQC7j9i2ZZo4EszWwYsoXRP4Zd99zrXzB4uc++NwN+Bq81sq5l1D0Rmkcp4L2UrM9ZmcO+ZXWmvtbhE5AR0alqfZy7tx9qd+7nr3aWUqNNXKqHSCzj7umc/As4EHHD4wbcESnfv+AKYeESzRVDQAs4SCNv2HeLMJ2fRvWUD3r7+JMI02iciAfDK7HQembqa20cl8avRnb2OIzVARQs4+7NlWz5wtpmdA4wF2vsubQQ+cc59doI5RYKWc457P1hGsXP89YI+KvpEJGCuHd6BtTtzePqbVDo3q885vVt6HUlqsEoXfoc55z4FPq2CLCJB6+0FW5idmskfJ/SkbZO6XscRkSBiZjwysSfpmQf59XtLad+kHj1bNfQ6ltRQ2h9KpIptycrl0amrGNapCZcPaut1HBEJQtER4bz4swE0rhvF9a+nsHu/9imQ8vlV+Pn2yn3HzH4wszQzSz/iSDv2XURCR0mJ4573lwHw+Pm9NcUrIlUmITaal69KZl9uITe8sZC8wmKvI0kN5M+WbXdS2sBxMqVds7OAmUccR+6tKxLS3py/ie/T93DfOd1p3UhTvCJStXq0bMiTF/dhyZZ9/O7D5VS2gVNChz/P+N1JaXF3pnOuoIryiASNTXsO8ufP1jAiKZ5LBrbxOo6IhIgze7bgV6M78/ev19GleSw3npzodSSpQfyZ6o0H3lHRJ3JsJSWOu99bRkSY8fj5vdHW0SJSnW47rRNje7fgsS/WMDtV247Kf/lT+C0COlZVEJFg8s+5G1mwMYsHxnWnZVwdr+OISIgxM564oA8d4uvx+49W6Hk/+ZE/hd8vKd25Y3RVhREJBukZB/jLl2s4rWtTLhjQ2us4IhKi6kSF88j4nmzOyuX56eu9jiM1xFGf8TOzr8o5nQN84dv2bCNw5K8Qzjk3JmDpRGqZ4hLH3e8vIyo8jD+f10tTvCLiqaGd4pnQtyUvzkxnQr9WdEyo73Uk8VhFI36dgaQjjihgs+9zHcu5rr1iJKS9+t0GFm7ay0Pje9CsQYzXcURE+N3YbkRHhnH/lBXq8pWjj/g559pXYw6RWm/97hz++tVaRndvxoS+rbyOIyICQNPYGO4e04UHpqzk46XbGa//P4U07dwhEgBFxSXc9d4y6kaF8+jEnpriFZEa5fLB7ejduiGPTF3N/rxCr+OIh1T4iQTApNnpLN2yj4fH96RprKZ4RaRmCQ8zHpnQk8wD+fz9q3VexxEPHbXwM7MSMysys6gyr4uPcRRVX3SRmmHtzhz+7+tUzurZnHG9W3gdR0SkXL1bx3HFSe14/fuNrNiW7XUc8UhFO3c8DDig6IjXIuJTWFzCr99bSv2YCP44QVO8IlKz3XVGFz5bvpPff7ScD28eRrj2Dw85FTV3PFjRaxGBF2eksXxbNs9f3p/4+tFexxERqVDDOpHcN7Ybv3xnCW8t2MwVJ7XzOpJUMz3jJ3KcVm3fz9PfpjKuT0vO7qUpXhGpHcb3bcnQxCb85Ys1ZOTkex1HqllFCziPPJ4bOs8w01cAACAASURBVOdmHX8ckdqhoKh0irdhnSgePreH13FERCrNzHh4fE/OemoWf/5sNX+/uK/XkaQaVfSM3wz8e6bPfO8PP5FAIrXBc9PXs2rHfiZdMYBG9aK8jiMi4pdOTetz48hEnp2+nguT2zAksYnXkaSaVFT4nVptKURqkRXbsnlu+nom9mvFGT2aex1HROS43HJqJyYv2cb9U1bw2e0jiIrQ01+hoKLmjpnVGUSkNsgvKuaud5fSuF4UfxjX3es4IiLHrU5UOA+P78E1/0zh5dnp3HJqJ68jSTU4rvLezJLMbJiZNQx0IJGa7OlvUlm7K4c/n9eLuLqa4hWR2u20rs0Y06MZz3ybypasXK/jSDXwq/Azs4vNbBOwBpgFDPCdjzezVDO7sAoyitQIS7fs44UZaVwwoDWjujXzOo6ISED8YVwPwsx46JOVXkeRalDpws/MxgNvA5uB+ylt5gDAOZcJrAauCHRAkZogr7CYu95bStPYGO4/R1O8IhI8WsbV4ZenJzFt9W6+WrnT6zhSxfwZ8bsPmOWcGwG8VM71+UCfgKQSqWGenLaO9bsP8Nj5vWhYJ9LrOCIiAfXzYR3o0iyWhz5ZRW6Bdl8NZv4Ufj2Adyu4vhPQ/JcEnYWb9vLyrHQuGdiGU7o09TqOiEjARYaH8cjEnmzbd4invkn1Oo5UIX8KvzwgpoLr7YF9J5RGpIbJKyzm7veW0qJhHX4/tpvXcUREqszA9o25cEBr/jF7A2t35ngdR6qIP4Xfd8Cl5V3wdfdeA3wbiFAiNcXbCzaTnnmQP5/Xi9gYTfGKSHD77dndqB8Twf2TV+CcP3s4SG3hT+H3INDDzKYD5/nOJZvZrcASoAHwx8DGE/FOYXEJr8zeQHK7RozsnOB1HBGRKte4XhT3ntmVBRuz+GDRNq/jSBWodOHnnFsEjAGa89/mjseAp4ECYIxzbnXAE4p45LPlO9i27xA3npzodRQRkWpzUXIb+reN40+frWZfboHXcSTA/FrHzzk32znXDegHXEzp1O9AoKtzbm4V5BPxhHOOF2emk5hQj1Fd1dAhIqEjLMx4ZEIvsg8V8vgXa72OIwHmzzp+rQ9/7Zxb6px7zzn3jnNuofM9CGBmJ1dFSJHqNjs1k9U79nPjyETCwuzYHxARCSLdWzbg6qHteXvBZhZt3ut1HAkgf0b8pplZ/NEumtkZwNQTjyTivUmz0mkaG834fi29jiIi4ok7R3emeYMY7vtoBUXFJV7HkQDxp/ALB74ub39eMxsHfAzMC1QwEa+s2JbNd+szuWZ4B6Ijwr2OIyLiifrRETwwrjurduzn9e83eR1HAsSfwu90oAnwuZnVO3zSzC4APgCmA+cENp5I9XtpVjr1oyO4bHBbr6OIiHjqrJ7NOblzAn//eh279ud5HUcCwJ+u3k3AaKAjMMXMos3sCkr37/0MGO+c038VUqttycpl6rLtXDa4LQ20bp+IhDgz46Fze1BQXMLDn67yOo4EgL9dvWspXdKlP6XTuq9ROtp3vnNOPd9S670yO53wMOPnw9p7HUVEpEZoH1+PW07pxNRlO5i1LsPrOHKC/Cr8oLSjFzgLSATeBC51zhUHOphIdcs6WMA7KVsY37cVLRrW8TqOiEiNcdMpHekQX48Hpqwgr1A/8muzoxZ+ZlZoZgXlHcBsoB5wGZBf5lp+dQUXCbQ3vt9EXmEJN4zs6HUUEZEaJToinD+O78nGPbm8ODPN6zhyAiIquPZvQBv1SUg4VFDMv77fyKiuTencLNbrOCIiNc7wpHjG9WnJ8zPSmNC3Fe3j6x37Q1LjHLXwc85dXY05RDz1/sItZB0s0GifiEgF7hvbjelrdnP/lBW8fs0gzLTAfW3j9zN+IsGmqLiEl2dvoG+bOAZ1aOx1HBGRGqtZgxjuOqMzs1Mzmbp8h9dx5DgcdcTPzEYCOOdmlX19LIffL1JbfLFyJ5uzcvnd2V3126uIyDFccVI73l+4lYc/WcXJnROI1dJXtUpFz/jNAJyZ1fEt1TKDip/5M991bXUgtYZzjpdmptMhvh6juzf3Oo6ISI0XER7GoxN7MfH5OTz5dSoPjOvudSTxQ0WF36kAZdbnO7Xq44hUr+/T97B8WzZ/mtiL8DCN9omIVEbfNnFcNqgt/5y7gfMHtKJHy5/s5io1VEXNHTMrei0SDF6amU58/SjO69/K6ygiIrXKPWO68sWKndw3eQUf3DSUMP3yXCuouUNC1uod+5m5LoOrh7YnJlJPKIiI+KNh3Uh+d3Y3Fm/ex39+2OJ1HKmkipo7HjiO+znn3B9PII9ItZk0K526UeH87KR2XkcREamVzuvfindTtvD4F2sY06MZTepHex1JjqGiZ/wePI77OUCFn9R42/Yd4uOl27lqSHvi6kZ5HUdEpFYyMx6Z0JOznprNnz9fwxMX9vE6khzDUad6nXNhx3FovkxqhVe/2wDAtSM6eJxERKR2S2oWy/UjO/L+wq0s2JDldRw5Bj3jJyEnO7eQtxds5tw+LWkVV8frOCIitd5tp3WiVVwd7pu8nMLiEq/jSAVU+EnIeXP+JnILirl+hLZnExEJhLpRETx4bg/W7TrAP3wzKlIzqfCTkJJXWMxrczYwsnMC3Vs28DqOiEjQGN29Gad3a8ZT01LZtu+Q13HkKFT4SUj5cNE2Mg8UcNNIjfaJiATag+eW7uLx0McrPU4iR6PCT0JGcYnj5dnp9GrVkCGJTbyOIyISdFo3qsvto5L4atUuvlm9y+s4Ug4VfhIyvl61iw2ZB7nx5I6YaYV5EZGqcO3wDiQ1rc8fPl7JoYJir+PIEVT4SUhwzvHizDTaNq7LmT2aex1HRCRoRUWE8ccJPdm69xDPfJvqdRw5QkULOP+PSuzk4YA8YCsw0zm3/USCiQTSDxv3smTLPh4e34OIcP2+IyJSlU7q2ITz+rfi5dnpnNe/FZ2axnodSXwqXfhRupOH83195DzZkeeLzexF4HbnnEPEYy/NTKNR3UguHNDG6ygiIiHhd2d345vVu7lv8grevv4kPWJTQ/gz9NEaWAq8CSQDDX3HIODfwBIgCRgA/Ae4GfhNIMOKHI/UXTl8s2Y3Vw1tT50obS4jIlId4utHc8+ZXZiXnsXkJdu8jiM+/hR+zwBpzrmrnHOLnHM5viPFOXclsAH4k3NusXPuCuAb4OoqyCzil0mz0omJDOPKIe29jiIiElIuHdiWPm3ieHTqarJzC72OI/hX+J1OaTF3NN8AY8q8/hRofxyZRAJmZ3Yek5ds4+LkNjSuF+V1HBGRkBIWZjw6oSdZBwv461drvI4j+Ff4lQB9Krjeh/8+63dYrt+JRALotTkbKC5xXKft2UREPNGzVUOuHNKef8/fzNIt+7yOE/L8KfwmA9eb2b1m9mN7jpnFmtlvget87zlsCLAuMDFF/Lc/r5B/z9/M2b1a0KZxXa/jiIiErLvO6ExC/Wjum7yC4hL1fHrJn8LvTuB74E/AXjPbYWY7gL3Ao8B833swsxhgH/D3wMYVqby35m/mQH4RN45M9DqKiEhIi42J5P5zurN8WzZvztvkdZyQVunCzzm3DxgJnA+8SmmH71LgH8B5wHDfe3DO5TnnfuGce/d4QpnZhWa20sxKzCy5zPnRZrbQzJb7/jytzLUZZrbWzJb4jqbl3Peon5fgkl9UzKvfbWBYpyb0at3Q6zgiIiHvnN4tGJEUzxNfrmX3/jyv44Qsv1aydaU+cs7d4Jw703fc6JybHOD1+lZQWkzOOuJ8JjDOOdcLuAp444jrlzvn+vqO3eXc91iflyAxZcl2dufka7RPRKSGMDMeOrcH+UUlPDJ1tddxQpY/CzgDYGaNgFFAB0qbOTYA3xwe7QsE59xq39915PnFZV6uBGLMLNo5l1/J+57Q56V2KClxTJqVTrcWDRiRFO91HBER8emYUJ+bTknk6W9SuXhgG4Z10v+jq5tfI35m9itKt2R7B3gc+AvwHrDNzO4MfLwKnQ8sPqJoe803zXu/HXuJ8PI+L0Hg2zW7Wb/7ADed3FErxYuI1DA3n5JIuyZ1uX/yCvKLir2OE3IqXfiZ2VXAE5Tu0HEx0BPoBVwELAaeMLMr/bjfNDNbUc4xvhKf7UFp4XljmdOX+6ZwR/iOK/z8/JHvucHMUswsJSMjo7LfltQAL81Ko1VcHc7u1cLrKCIicoSYyHAeHt+T9MyDTJqZ7nWckOPPVO+dwHfAqc65siX6SjP7CJgO/Ap4vTI3c86d7sff/SMzaw18BFzpnEsrc79tvj9zzOwtSreS+0mWo32+nHyTgEkAycnJ6j2vJRZu2ssPG/fywDndiQz3a0BbRESqycmdExjbqwXPTl/P+L6taNtES25VF39+MnYB3j2i6APAd+5d33uqjJnFAVOB3zrn5pQ5H2Fm8b6vI4FzKG0QqdTnJXhMmpVGwzqRXDywjddRRESkAvef052IMOOBj1cQ2P5QqYg/hV8O0LKC66187zlhZjbRzLZSugj0VDP70nfpVqATcP8Ry7ZEA1+a2TJKp6K3AS/77nWumT18jM9LEEjLOMBXq3Zx5ZB21Iv2u29JRESqUfOGMdw5ujMz1mbw5cqdXscJGVbZKtvM3gQmABOcc9OOuDaK0l07JjvnjvpsXW2VnJzsUlJSvI4hx/DbD5fxwaJtzL33NOLrR3sdR0REjqGouIRxz85hX24B0351sn5pDxAzW+icSy7vmj8jfvdSukvHl2a21Mz+4zuWAF/5rv32xOOK+G93Th4fLNzGBQNaq+gTEaklIsLDeGRCT3Zk5/F/07TLa3XwZ+eOrUBf4EkgChjvO6KBvwH9fO8RqXb/nLORwpISrh/R0esoIiLihwHtGnHpoDa8Omcjq3fs9zpO0PN35449zrlfO+e6Oefq+I5uzrl7nHN7qiqkSEUO5BfxxrxNnNmjOR3i63kdR0RE/HTPmK40rBPJfZNXUFKiRo+qpPUupNb7z4LN5OQVccNIjfaJiNRGjepFce9ZXVm4aS/vL9TkYVU66lOUZvbAcdzPOef+eAJ5RPxSWFzCP77bwOAOjenXtpHXcURE5Dhd0L8176Vs4c+fr2Z092Y0qhfldaSgVFH7zIPHcT8HqPCTavPJ0u3syM7jTxN7eR1FREROQFiY8ciEXox9ejaPfb6Gxy/o7XWkoHTUqV7nXNhxHOHVGV5Cm3OOl2am07lZfU7pkuB1HBEROUFdmsdy7fAOvJOyhZSNWV7HCUp6xk9qrRnrMli7K4cbRiZiZl7HERGRALh9VBItG8Zw3+QVFBWXeB0n6Kjwk1rrpZlpNG8Qw7l9KtpQRkREapN60RE8MK4Ha3bm8M+5G72OE3RU+EmttHTLPualZ3Ht8A5EReg/YxGRYDKmRzNO69qUJ79ex47sQ17HCSr6iSm10qRZ6cTGRHDJoDZeRxERkQAzMx46twfFzvHwJ6u8jhNUVPhJrbMx8yCfr9jBz05qR2xMpNdxRESkCrRpXJfbTkvi8xU7mb52t9dxgoYKP6l1XvkunYiwMH4+tL3XUUREpApdN6IDHRPq8YcpK8krLPY6TlA4rsLPzGLMrJWZaXVFqVaZB/J5L2UrE/u1ommDGK/jiIhIFYqOCOeR8T3ZnJXLizPTvI4TFPwq/MxsuJnNBnKAzcBw3/l4M/vGzM6ogowiP3r9+03kF5VwvbZnExEJCUM7xXNWz+b8Y/YGsnMLvY5T61W68DOz4cA3QHPgFeDHhdOcc5m+19cGOqDIYbkFRbz+/UZGd29Gp6b1vY4jIiLV5PZRSeTkF/HqnA1eR6n1/BnxewRYBfQE7i/n+kxgYCBCiZTn3R+2sC+3kJtO1mifiEgo6daiAWN6NOPVORvIPqRRvxPhT+GXDPzTOZdP6Z68R9pK6WigSMAVFZfw8uwNJLdrxIB2jb2OIyIi1ez2UUnk5BXxLy3qfEL8KfxKKL/gO6wlkHticUTKN3X5DrbtO8QNerZPRCQk9WjZkNO7NeMf320gJ0+jfsfLn8LvB+Dc8i74unt/BswNRCiRspxzvDQznY4J9Ti9WzOv44iIiEfuGJVE9qFCjfqdAH8Kvz8Bp5jZ65RO+wK0MbNzgFlAB997RAJqzvo9rNqxnxtHdiQszI79ARERCUq9WjfktK5NeeW7DRzIL/I6Tq1U6cLPOfcNcBlwNvCZ7/SrwMdAZ+Ay59y8gCeUkPfSrDQSYqOZ0K+V11FERMRjd4xKYl9uIa9/v9HrKLVShD9vds69a2afAqMpLfbCgPXAl865A1WQT0Lcim3ZzE7N5DdndiU6ItzrOCIi4rE+beI4pUsCr8zewFVD2lMv2q9SJuT5vXOHcy7XOTfFOfdX59zjzrkPVPRJVZk0K516UeFcNrit11FERKSGuH1UElkHC3hz3iavo9Q62qtXaqwtWblMXb6Dywa3pWGdSK/jiIhIDdG/bSNGJMUzaVY6uQV61s8fRy38zKzEzIr9PPSvLwHzj+82YMA1wzt4HUVERGqYX56exJ6DBbw1f7PXUWqViibGH6bidftEqszegwW888MWxvdtRYuGdbyOIyIiNcyAdo0Z1qkJL85M5/LB7agTpefAK+OohZ9z7sFqzCHyP96Yt4lDhcVasFlERI7qjlGdueil73l7wWbNDlWSnvGTGievsJh/zt3IaV2b0qV5rNdxRESkhhrUoTEndWzMizPTyCss9jpOrXDUET8zGwngnJtV9vWxHH6/yPF6b+FWsg4WaLRPRESO6Y5Rnbn05Xn8Z8Fmrh6mUb9jqegZvxmAM7M6zrmCw68reL/5rmuSXY5bcYnjldnp9GkTx+AOjb2OIyIiNdyQxCYM6tCYF2amccmgtsREqgypSEWF36kAvqIP4DTU7CFV7MuVO9m0J5d7z+yKmbZnExGRY7tjVBKXvzKf91K2cMWQ9l7HqdEqKvz6AF8cfuGcm1HlaSSkOed4aWYa7ZvU5Ywezb2OIyIitcTQxCYkt2vE8zPSuGhgG+30VIGKmjueBJIPv/Ct03dZ1UeSUDUvPYulW7O5fmRHwsM02iciIpVjZtw+Kokd2Xm8v3Cr13FqtIoKv31AfJnX+kksVeqlWWnE14/i/P6tvY4iIiK1zIikePq1jeP56WkUFJV4HafGqmiqdx5wn5m1A7J9584zs04VfMY55/4YsHQSMtbs3M+MtRncNbqzHswVERG/mRl3jEri6td+4MNFW7lkkPZ4L09Fhd8twD+BOygdGXTAeb7jaBygwk/8NmlWOnUiw7liSDuvo4iISC11cucE+rRuyLPT13P+gNZEhmu54iMd9V/EObfROXcKEAO0pnSq9zagTQWHymvx2/Z9h/h4yXYuGdSGuLpRXscREZFaysy44/Qktu49xEeLtnkdp0aqaMQPAOdcEbDdzB4CZjrn9C8pAfXqdxtwwLXabkdERE7QqV2a0qtV6ajfef1bEaFRv/9R6X8N59xDzrkVAGaWYGYDzSzZzBKqLp4Eu+zcQt5esJlxvVvQulFdr+OIiEgtd7jDd3NWLpOXbPc6To3jVxlsZkPMbB6wk9Lmj/nATjOba2YnVUVACW5vzt/EwYJibhiZ6HUUEREJEqd3a0r3Fg14bvp6iorV4VtWpQs/X2H3LdAVeIHS5/1u933dHZhuZoOrIqQEp7zCYl6bs5ERSfF0b9nA6zgiIhIkDo/6bcg8yCfLNOpX1jGf8SvjEWAHMNQ5t7PsBTN7BJjre8/owMWTYDZ58TYyD+Rz08l9vY4iIiJB5ozuzejaPJZnvl3PuX1aaWMAH3+megcDLx1Z9AH4zk3yvUfkmEpKHJNmpdOzVQOGJjbxOo6IiASZsLDSUb/0jIN8qlG/H/lT+DnfcTSaRJdK+3r1LtIzD3LjyETM9FuYiIgE3pk9mtOlWemoX3FJRSVM6PCn8PsBuNHM4o+84Dt3I7AgUMEkeDnneHFmGm0a1+Gsns29jiMiIkEqLMy4bVQn1u8+wOcrdngdp0bw5xm/B4BvgLVm9jqw1ne+K3AFUNf3p0iFUjbtZfHmfTx0bg+tryQiIlXqrJ4t6NQ0lae/SeXsni0IC/Fn/fxZx28OcAawgdJt3J73HbcDacAZzrm5VRFSgstLM9NoVDeSC5Nbex1FRESCXHiYcdtpnVi36wBfrPxJm0LI8Wu4xTk3yzmXDLQAhviOFs65Qc652VURUILL+t05TFu9myuHtKdulD8DziIiIsfnnN4t6ZhQj6e/SaUkxJ/1q1ThZ2Z1zWyhmd0E4Jzb5Zyb7zt2VW1ECSaTZqUTExnGlUPaeR1FRERCxOFRvzU7c/hqVWiXLZUq/JxzuUBH1LkrJ2DX/jw+WryNi5Lb0KR+tNdxREQkhIzr3ZIO8aWjfs6F7qifP1O904FTqiiHhIBX52yguMRx3fCOXkcREZEQExEexi2ndmLVjv1MW73b6zie8afwux3oY2ZPmllnM9MDWlJpOXmFvDVvM2f1akHbJnW9jiMiIiFoQt+WtGtSl6e+WReyo37+FH4bgM6UFoCrgXwzKzjiyK+SlFLrvb1gMzn5Rdw4UqN9IiLijcOjfiu27Wf62tAc9fNn1O7fVLxzh0i5CopK+Md3Gxia2ITereO8jiMiIiFsYr9WPPNtKk9NS+XULk1DbveoShd+zrmrqzCHBLEpS7axa38+f7mgj9dRREQkxEWGh3HLKZ2498PlzFiXwaldmnodqVpp2wSpUiUljkmz0unaPJaRST/Z7U9ERKTande/Na3i6vDUtNDr8PWr8DOzODN71MyWmlm271jqO9eoqkJK7TV97W5Sdx/gppMTQ244XUREaqaoiDBuPjWRJVv2MTs10+s41arShZ+ZJQLLgN8C4cA0SvfuDfedW+Z7j8iPXpqZTsuGMYzt3cLrKCIiIj+6YEBrWjaM4akQW9fPnxG/p4FGlO7J29M5d75z7jznXE9gDBDne48IAIs272XBxiyuHdGRyHA9VSAiIjVHdEQ4vzglkYWb9jI3bY/XcaqNPz+NTwb+zzk37cgLzrmvKS36Tg5UMKn9Js1Mp2GdSC4Z2MbrKCIiIj9x0cA2NG8QE1LP+vlT+B0AMiq4vsv3HhHSMw7w5aqdXHFSO+pFa61vERGpeaIjwrnp5I4s2JjFvPQsr+NUC38Kv7eBn5lZ1JEXzCwGuAJ4K1DBpHZ7efYGIsPDuGpoe6+jiIiIHNUlg9rSNDaap75Z53WUauFP4fcJEAksMbPbzexMMxtjZncACylt8vjEzIaWPaoitNRsGTn5fLBoKxcMaE1CbLTXcURERI4qJjKcm05OZF56FvPTg/9ZP3/m4Mo+2/d//HcXDzvKe8z3nvDjiya11ZvzNlFYXML1I7Q9m4iI1HyXDW7L8zPSePrbVP7dsYnXcaqUP4Xfz6sshQQN5xxTlmxjaGITOsTX8zqOiIjIMZWO+nXkkamrSdmYRXL7xl5HqjL+bNn2r6oMUpaZXQg8CHQDBjnnUnznRwOPAVFAAXC3c+5b37UZQAvgkO82Zzjndh9x30HApMMvgQedcx9V6TcTYlZs28/GPbn84hQt6SgiIrXHZYPb8sKMNJ76JpU3rh3sdZwqU1MXV1sBnAfMOuJ8JjDOOdcLuAp444jrlzvn+vqO3fzUCiDZOdcXOBN4yczUchpAnyzbTmS4MaZHc6+jiIiIVFrdqAhuGNmR2amZLNq81+s4VaZGFn7OudXOubXlnF/snNvue7kSiDGzSncPOOdynXNFvpcx/Pc5RQmAkhLHp0u3MyIpgbi6P2n+FhERqdF+dlI7GteL4qlpqV5HqTI1svCrpPOBxc65/DLnXjOzJWZ2vx1lY1gzG2xmK4HlwE1lCkE5QYu37GV7dh7j+mh7NhERqX3qRUdw3YgOzFyXwZIt+7yOUyU8K/zMbJqZrSjnGF+Jz/YAHgduLHP6ct8U8AjfcUV5n3XOzXfO9QAGAr/1rUFY3t9xg5mlmFlKRkZF61bLYZ8s3UF0RBind2vmdRQREZHjcuWQ9sTVjeTpb4Jz1M+zws85d7pvz98jjykVfc7MWgMfAVc659LK3G+b788cSheSHnSMv381cBDoeZTrk5xzyc655ISEBP++uRBUXOL4dNkOTuvalNiYSK/jiIiIHJf60RFcN7wD367ZzfKt2V7HCbhaNdVrZnHAVOC3zrk5Zc5HmFm87+tI4BxKGzmO/HyHw80cZtYO6AJsrIboQW9++h4yD+RzTu+WXkcRERE5IVcNbU/DOpE8FYSjfpUu/Mysh5mdd8S5U83sGzNbaGZ3BSqUmU00s63AEGCqmX3pu3Qr0Am43/cs3xIzawpEA1+a2TJgCbANeNl3r3PN7GHf54cDS81sCaWjhjc75zIDlTuUfbJsB3Wjwjmta1Ovo4iIiJyQ2JhIrh3egWmrd7FiW3CN+plzlWtsNbOpvvef7XvdClgD5AEZlI6e/dw593oVZfVMcnKyS0lJ8TpGjVVYXMLAR6dxcucEnrqkn9dxRERETlj2oUKGP/4tQxOb8NIVyV7H8YuZLXTOlRvan6ne/sDMMq8vp3Q7tr7Oue7A58Atx51Saq3v1meyL7eQcZrmFRGRINGwTiTXDOvAlyt3sXrHfq/jBIw/hV8jYFeZ12cBMw83VQCfAJ0DFUxqj0+Wbic2JoIRneO9jiIiIhIw1wzrQGx0BM98GzzP+vlT+GUCrQHMrB6lz99NK3M9Av/2/pUgkFdYzNcrd3Fmj+ZER4R7HUdERCRgGtaN5Oph7fls+U7W7szxOk5A+FP4zQZ+4Wvw+D8gEv6/vTsPk6q68z/+/tAsDQg2+yoCsqggqEHjEhXFXTDJZDGZqDGJMTHRjJPJLI46WWeyTCarmkj8xbjEJBOTidOgEo2KCpioA92ICMim0M0qKHtD9/n9cW9rpaxeq7pvLZ/X89RTVeeee+t7bhW3v5xzz738b8ryiUSTKqyEzF+5lV0HDjFrqod5CA+agwAAIABJREFUzcys+HzqPWPo3b2MHxVJr19bEr9/BfYBDwCfAr4TQlgFIKkM+CB/fQ6glYDKqhr69+7OaUcNSDoUMzOznKvo1Z2Pnzaah5bWsmpz4ff6tTrxCyGsBY4GjgfGhBBuTFncC7gW+GZuw7N8trfuEH9avoWLJg+la1lBXRLSzMys1a4+Yyw9u5Xx48dfSTqUrLXpr3UI4VAIoTqEsD6tfFcI4cEQwrqcRmd57U/Lt7DvYL2Hec3MrKj1792dK08dTWV1Da9s2Z10OFlpU+Inqb+kr0taIGmVpFPj8gGS/k3S0R0TpuWjyqoahvTtwUmj+ycdipmZWYf69BljKO9axm1PFHavX1vu3HEE0V0x/gnoA4wFegKEELYDH8XX8SsZb+4/yJMrtnLxccMo66KkwzEzM+tQAw7rwRWnHsmDSzayZmvh9vq1pcfvO0A50Tl+5wDpf+0fjMutBDy6bDN19Q0e5jUzs5Lx6TPG0r1rF257YnXSobRbWxK/84AfhRCWA5nu87YWOCInUVneq6yuYURFT044oiLpUMzMzDrFoD49+Ni7j+QPSzayfvuepMNpl7Ykfr2BLc0sPyzLWKxAvL6njmdWbWPW1OFIHuY1M7PS8Zkzx9K1i7i1QGf4tiXxWwGc0szyi4EXswvHCsEjL27iUENg5pRhSYdiZmbWqQb3LeejJ4/i94s38trre5MOp83akvjdAVwu6ZNA4725gqQ+kr4PTAduz3F8lofmVNcwdmBvJg3vm3QoZmZmne7a6UdR1kUFOcO3LRdw/gkwG7gTWBMXPwDsAP4O+HEI4b6cR2h5Zcub+1m0ZjszPcxrZmYlakjfcj5y0hE88MIGNuworF6/tl7A+TrgdKLk72HgL8BPgDNCCDfkPjzLNw8trSUEmOVhXjMzK2HXTj+KLhK3P1lYM3y7tnWFEMIiYFEHxGIFoLK6lqOH9mH8kD5Jh2JmZpaYYYf35MMnjeQ3z73G588ex4iKnkmH1Cq+waq12sad+3hh/Q5fu8/MzAy4dvo4AH5aQL1+bb1l29WSFkraJOmApLq0x4GOCtSSN7e6BsCzec3MzIARFT354LuO4DfPvUbtG/uSDqdVWj3UK+m7wN8DNURDvTs7KijLT5VVtUwdeThHDuiddChmZmZ54XPTj+K3z7/GHfPX8JVLJyUdTovaco7fJ4kmdLw3hFDfQfFYnlq3bQ9LN77BTRcfk3QoZmZmeeOI/r34wIkjuf8vr3Lt9KMY0rc86ZCa1ZahXgGVTvpK05x4mPcSD/OamZn9lc+fPY76hsBP5+f/uX5tSfzmASd1VCCW3yqrajlpdD+GF8isJTMzs84yakAv3n/CCO7/86ts2bU/6XCa1ZbE7wvASZK+LGlERwVk+WfFpl2s2LzLs3nNzMyacN3Z4zhY38Ds+Wtarpygtty5YwtwL/BvwKuSDnpWb2mYU11DF8FFkz3Ma2Zmlsnogb153/EjuO/P69m6K3/TobbM6v0acBPRrN7n8azekhBCYE51LaceNYBBfXokHY6ZmVneuu6ccfxhyUbufHoNN+bpZMi2zOr9DJ7VW3KW1bzJ2m17+MyZY5MOxczMLK+NHXQYl04dzj2L1nPNmWMZcFj+dZi05Ry/cjyrt+RUVtXQtYu4cPLQpEMxMzPLe9edM479h+q585m1SYeSUVsSvz/iWb0lpXGY94zxA6no1T3pcMzMzPLeuMF9mDllOPcsXMeOPXVJh/MObUn8rgPeJemrntVbGv7v1Z1s3LnPs3nNzMza4PpzxrH3YD13PpN/M3zbkvhtACYDN+NZvSWhsqqG7l27cN6xQ5IOxczMrGBMGNKHiycP4+6F69m5N796/doyueOXQOioQCy/1DcE5i6t5eyJg+hT3i3pcMzMzArK9TPGMXdpLT9/Zi1fPH9i0uG8pdWJXwjhqg6Mw/LMn9duZ+uuAx7mNTMza4ejh/blwklDuWvBOj51xlgO75kfnShtGeq1EjKnupZe3cs45+jBSYdiZmZWkL4wYzy7DhzirgX5M8O3yR4/SWcChBCeSn3fksb6VrgO1jfw8NJazj1mCL26t+VsADMzM2t07PC+nH/sEH7+zFo++Z4x9M2DU6ea+6v+JBAk9Qwh1DW+b6a+4uVlOYvOErHglW3s2HuQmVN8izYzM7NsfGHGeP740mbuXrCO62eMTzqcZhO/swHipO+t91b85lTX0qe8K2dNHJR0KGZmZgVt8ojDOfeYwdz5zFquOn104hMmm0z8Qgjzm3tvxenAoXrmvbiJCyYPpUdXd96amZll6wszxnPprQu4Z9F6Pn/2uERjafXkDkmPS5rRzPKzJT2em7AsKfNXbGXXgUOezWtmZpYjU0ZWcPbEQdz59Br2HDiUaCxtmdU7HWjuSr6DgbOyisYSV1ldS79e3TjtqAFJh2JmZlY0vjBjPDv2HuSeResTjaOtl3NpbnLHUcDuLGKxhO2tO8RjL23mouOG0a3MV/oxMzPLlRNG9ePMCYN4seaNRONo9lodkq4ArkgpulHSJzJUrQBOAB7NYWzWyR5/eQv7DtYza4qHec3MzHLtp5efmPhl0lr69P5A49zjAAwF+qTVCcAeolu63ZzT6KxTVVbVMLhPD04e0z/pUMzMzIpO0kkftJD4hRB+CPwQQFIDcEMI4f7OCMw61679B3lixVb+9uRRlHVR0uGYmZlZB2jLvXp90lcRe/SlzdQdavBsXjMzsyLmZM6AaJh3REVPThxVkXQoZmZm1kGc+Bk79tTx9KptzJwyDMnDvGZmZsXKiZ/xyLJNHGoIHuY1MzMrck78jDnVNYwZ2JtJw/smHYqZmZl1ICd+JW7Lrv0sWr2dWR7mNTMzK3pO/Ercw0s30RBgpod5zczMil6rEj9JFZJOlnRUM3XGSLoyd6FZZ6isqmHikD5MGJJ+XW4zMzMrNi0mfpL+HdgMLAJWSvqLpBMzVD0NuCvH8VkHqtm5j+fX72DW1GFJh2JmZmadoNnET9KHgRuB+cDngf8ARgELJV3e8eFZR5pbXQvATN+b18zMrCS0dOeOG4AnQgjnNxZI+h7wG+AXkgaGEH7QkQFax6msruG4EYczemDvpEMxMzOzTtDSUO/RwO9SC0IIO4ALgXuA/5L09Q6KzTrQum17qN7whod5zczMSkhLPX4NmQpDCA3AJyXtAG6S1A/4c66Ds44zd2k0zHuJh3nNzMxKRkuJ30rgPcDtmRaGEP5B0pvAl4GZOY7NOlBlVQ3TjuzHiIqeSYdiZmZmnaSlod6HgfdKGtBUhRDCV4G/B47IZWDWcVZu3sXLm3Yxc4qHec3MzEpJSz1+PwdeBwYD25uqFEL4oaT1wNQcxmYdZE5VDV0EFzvxMzMzKynNJn4hhI3Aba3ZUAjhD8AfchGUdZwQAnOqazll7AAG9ylPOhwzMzPrRO2+ZZukcklXShqSy4CsYy2reZM12/Ywy7doMzMzKznZ3Kv3cKI7dUzKUSzWCSqra+jaRVw4aWjSoZiZmVknyybxA1BOorBOEUJgTlUt7xk/kH69uycdjpmZmXWybBO/kJMorFMsfm0nG3fuY5av3WdmZlaS3ONXQiqraujetQvnTfJpmWZmZqUom8RvKzAGWJCjWN4i6UOSlklqkDQtpfw8SS9IWho/n5Oy7ElJKyQtiR+Dm9n+KEm7JX0p17Hnq/qGwNzqWqZPGETf8m5Jh2NmZmYJaOk6fk2Kb9u2PoexpHoR+BvgjrTybcCsEEKNpMnAPGBEyvKPhRCeb8X2v090ceqS8dy619my64Bn85qZmZWwdid+HSmEsBxAUnr54pS3y4ByST1CCAdau21J7wPWAHtyEGrBqKyqoWe3MmYc02RHqJmZmRW5bM/xS9IHgMVpSd9d8TDvLUrPGgFJvYF/Br7aWUHmg4P1DTz84iZmHDOYXt3zMtc3MzOzTpBYFiDpMSDTxeRuCiE82MK6k4BvA+enFH8shLBRUh/gd8AVwD1pq34V+H4IYXeGvDD9M64BrgEYNWpUs3Xz3cLV23l9T52Hec3MzEpcYolfCOHc9qwnaSTwP8CVIYTVKdvbGD/vknQ/cDLvTPzeDXxQ0neACqBB0v4Qwq0Z4psNzAaYNm1aQV+2Zk5VDX16dOWsCYOSDsXMzMwSVFDjfpIqgLnAjSGEBSnlXYGKEMI2Sd2AmcBj6euHEM5IWecrwO5MSV8xOXConkeWbeL8SUMp71aWdDhmZmaWoLw8x0/S+yVtAE4F5kqaFy+6DhgH3JJ22ZYewDxJ1cASYCPws3hbl0r6Wue3Ij88tXIbu/YfYubUYUmHYmZmZglTCAU9itkppk2bFp5/vjVXick/X/jVYp5atZXnbjqXbmV5meebmZlZDkl6IYQwLdMyZwJFbF9dPY8t38xFk4c56TMzMzMnfsXs8Ze3sLeunlke5jUzMzOc+BW1yqoaBvXpwbvHDEg6FDMzM8sDTvyK1K79B3l8xRYuOW4YZV2av2ahmZmZlQYnfkXqseWbqTvU4GFeMzMze4sTvyJVWVXLiIqenHBEv6RDMTMzszzhxK8I7dxbx1Mrt3LJlGF08TCvmZmZxZz4FaFHXtzEoYbArCm+N6+ZmZm9zYlfEZpTXcvoAb2YPKJv0qGYmZlZHnHiV2S27jrAwtXbmDV1OJKHec3MzOxtTvyKzMMv1tIQYKaHec3MzCyNE78iU1lVw4QhhzFxaJ+kQzEzM7M848SviNS+sY/n1u3wpA4zMzPLyIlfEZlbXQvAzKlO/MzMzOydnPgVkcqqGiaP6MuYgb2TDsXMzMzykBO/IrF++x6qNrzhYV4zMzNrkhO/IjEnHua9ZIrvzWtmZmaZOfErEpVVNZw4qoKR/XolHYqZmZnlKSd+RWDV5l28vGkXszypw8zMzJrhxK8IVFbXIsElx3mY18zMzJrmxK/AhRCYU13DKWMGMLhvedLhmJmZWR5z4lfgXqp9kzVb9zBzqnv7zMzMrHlO/ApcZVUtZV3ERZOd+JmZmVnznPgVsMZh3veMG0j/3t2TDsfMzMzynBO/ArbktZ1s2LHPs3nNzMysVZz4FbDKqlq6l3Xh/ElDkg7FzMzMCoATvwLV0BCYu7SGsyYOom95t6TDMTMzswLgxK9APbfudTa/ecDDvGZmZtZqTvwKVGV1DT27lXHuMYOTDsXMzMwKhBO/AnSovoGHlm7inGMG06t716TDMTMzswLhxK8ALVy9ndf31DFriod5zczMrPWc+BWgOdU1HNajK9MnDko6FDMzMysgTvwKzIFD9Tzy4ibOnzSE8m5lSYdjZmZmBcSJX4F5euU23tx/yMO8ZmZm1mZO/ApMZXUNFb26cfq4gUmHYmZmZgXGiV8B2VdXz2MvbeaiyUPp3tVfnZmZmbWNs4cC8sSKLeypq2emh3nNzMysHZz4FZDKqhoGHtaDU8YOSDoUMzMzK0BO/ArE7gOHePzlLVxy3FDKuijpcMzMzKwAOfErEI+9tJkDhxp8b14zMzNrNyd+BaKyqoZhh5dz4qh+SYdiZmZmBcqJXwHYubeOp1ZtZeaUYXTxMK+ZmZm1kxO/AjBv2SYO1gcP85qZmVlWnPgVgDnVtRw5oBfHjTg86VDMzMysgDnxy3Pbdh9gwSvbmDllGJKHec3MzKz9nPjluYeX1tIQ8DCvmZmZZc2JX56rrKpl/ODDmDikT9KhmJmZWYFz4pfHat/Yx3PrX2fW1OEe5jUzM7OsOfHLY3OrawkBZk4ZlnQoZmZmVgSc+OWxyupaJg3vy9hBhyUdipmZmRUBJ3556tXte6l6bacndZiZmVnOOPHLU3OW1gBwyXEe5jUzM7PccOKXpyqrajlhVAVH9O+VdChmZmZWJJz45aFXtuxmee2bzJriYV4zMzPLHSd+eWhOdQ0SXOLZvGZmZpZDTvzyTAiByqoaTh7dnyF9y5MOx8zMzIqIE788s7x2F6u37vFsXjMzM8s5J355prK6hrIu4qLJQ5MOxczMzIqME788EkJgTnUNp48byIDDeiQdjpmZmRUZJ355pGrDG7z2+j7fos3MzMw6hBO/PFJZVUO3MnHBJA/zmpmZWe458csTDQ2BudW1nDVhMIf37JZ0OGZmZlaEnPjliefX72DTm/uZNdXDvGZmZtYxnPjlicqqGsq7deHcY4YkHYqZmZkVKSd+eeBQfQMPLa1lxtFD6N2ja9LhmJmZWZFy4pcHtu2uY+yg3r5os5mZmXWovEz8JH1I0jJJDZKmpZSfJ+kFSUvj53NSlj0paYWkJfFjcIbtjpa0L6XOTzurTc0Zeng5v/3saVzoizabmZlZB8rXccUXgb8B7kgr3wbMCiHUSJoMzANGpCz/WAjh+Ra2vTqEcHzuQjUzMzMrDHmZ+IUQlgNISi9fnPJ2GVAuqUcI4UAnhmdmZmZWkPJyqLeVPgAsTkv67oqHcG9Retb4tjGSFkuaL+mMTojTzMzMLC8k1uMn6TEg00ltN4UQHmxh3UnAt4HzU4o/FkLYKKkP8DvgCuCetFVrgVEhhO2S3gX8QdKkEMKbGT7jGuAagFGjRrW2WWZmZmZ5K7HEL4RwbnvWkzQS+B/gyhDC6pTtbYyfd0m6HziZtMQv7h08EL9+QdJqYALwjvMCQwizgdkA06ZNC+2J1czMzCyfFNRQr6QKYC5wYwhhQUp5V0kD49fdgJlEE0TS1x8kqSx+PRYYD6zpjNjNzMzMkpaXiZ+k90vaAJwKzJU0L150HTAOuCXtsi09gHmSqoElwEbgZ/G2LpX0tXj9M4FqSVXAA8BnQwivd17LzMzMzJKjEDyK2ZJp06aF559v6SoxZmZmZsmT9EIIYVqmZXnZ42dmZmZmuefEz8zMzKxEOPEzMzMzKxFO/MzMzMxKhBM/MzMzsxLhxM/MzMysRDjxMzMzMysRTvzMzMzMSoQTPzMzM7MS4cTPzMzMrET4lm2tIGkrsL6ZKgOBbZ0UTr5x20tXKbffbS9Npdx2KO32F1rbjwwhDMq0wIlfDkh6vql74hU7t7002w6l3X633W0vRaXc/mJqu4d6zczMzEqEEz8zMzOzEuHELzdmJx1Agtz20lXK7XfbS1Mptx1Ku/1F03af42dmZmZWItzjZ2ZmZlYinPg1Q9KFklZIekXSv2RYfpWkrZKWxI+rU5aNkvRHScslvSRpdGfGngtZtv87kpbF7f+RJHVu9Nlpqe1xnQ/H3+0ySfenlH9c0qr48fHOizo32tt2ScdLWhSXVUu6rHMjz14233u8rK+kjZJu7ZyIcyvL331BH/OybHtRH+8kfT/lOL9S0s6UZUV9vGuq7QV9vAsh+JHhAZQBq4GxQHegCjg2rc5VwK1NrP8kcF78+jCgV9Jt6qz2A6cBC+JtlAGLgOlJtynHbR8PLAb6xe8Hx8/9gTXxc7/4db+k29RJbZ8AjI9fDwdqgYqk29QZbU9Z/kPg/qaOC/n8yLb9hXzMy/J3X/THu7T61wM/j18X/fGumbYX7PHOPX5NOxl4JYSwJoRQB/waeG9rVpR0LNA1hPAoQAhhdwhhb8eF2iHa3X4gAOVE/5B6AN2AzR0SZcdoTds/DdwWQtgBEELYEpdfADwaQng9XvYocGEnxZ0L7W57CGFlCGFV/LoG2AJkvIBonsrme0fSu4AhwB87Kd5ca3f7i+CYl813XwrHu1QfBX4Vvy6F412qt9peyMc7J35NGwG8lvJ+Q1yW7gNxN+8Dko6IyyYAOyX9XtJiSf8pqayjA86xdrc/hLAIeILof0C1wLwQwvKODjiHWtP2CcAESQskPSvpwjasm8+yaftbJJ1M9IdwdYdFmnvtbrukLsB/Af/YKZF2jGy++0I/5rW77SVyvANA0pHAGODxtq6bp7Jpe+qygjreOfFrWqZzNNKnQFcCo0MIU4DHgLvj8q7AGcCXgJOIupGv6pgwO0y72y9pHHAMMJLoH9E5ks7swFhzrTVt70o09DOd6H+Bd0qqaOW6+SybtkcbkIYB9wKfCCE0dFCcHSGbtn8OeCiE8BqFK5v2F/oxr91tL5HjXaOPAA+EEOrbsW4+yqbt0QYK8HjnxK9pG4AjUt6PBGpSK4QQtocQDsRvfwa8K2XdxXH38SHgD8CJHRxvrmXT/vcDz8bDPbuBh4FTOjjeXGqx7XGdB0MIB0MIa4EVRH8UWrNuPsum7UjqC8wFbg4hPNsJ8eZSNm0/FbhO0jrgu8CVkr7V8SHnVLa/+0I+5mXT9lI43jX6CG8P87Z13XyUTdsL93iX9EmG+fog+t/dGqKu3caTPiel1RmW8rrxHz9EJ4xWAYPi93cBn0+6TZ3Y/suIegC7Ep3v8idgVtJtynHbLwTujl8PJBouGEB0kvNaohOd+8Wv+yfdpk5qe/f4u74h6XZ0dtvT6lxFYU7uyOa7L+hjXpZtL/rjXVxvIrCO+Pq/cVnRH++aaXvBHu8SDyCfH8DFwEqicfub4rKvAZfGr78JLIt/LE8AR6esex5QDSwFfgF0T7o9ndX++I/AHcBy4CXge0m3pQPaLuB7cfuWAh9JWfeTwCvx4xNJt6Wz2g5cDhwElqQ8jk+6PZ31vads4yoKMPHLtv2FfszL4ndf9Me7+P1XgG9lWLeoj3dNtb2Qj3e+c4eZmZlZifA5fmZmZmYlwomfmZmZWYlw4mdmZmZWIpz4mZmZmZUIJ35mZmZmJcKJn5llTdLVkoKkkUnH0pkkHSXpIUk74vZfl3RMLZH0jKTHko6jLSR1jffvzUnHYlbonPiZFZH4XqkHJTV5s3BJfxf/EZ3VmbEVqbuBaUTX+boC+GOi0ZiZtcCJn1lxuZfoavQfaabO5cA24JEcfu5dQM8QwoYcbjOvSSoHTgPuCSH8MIRwXwhhZdJxmZk1x4mfWXGZC7xOlNy9g6QJRD1Uvw4hHMz2wyT1Bggh1IcQ9me7vQIzmOhuDjvbs7IivXIbkplZ85z4mRWREEId8N/AyZLGZ6hyRfx8b2OBpPdJqpS0QVJd/HxbfANyUup9Ix4iPlbSLyRtI7p/ZcZz/CSdJek3ktZLOiBps6R7JA1P227jumdI+rakTZL2SZonaXR6AySNk3RfXO+ApLWS7mhMQuM6fSV9V9K6uM6r8bbLW7Mf41iekLRb0i5Jj0p6d+q+ANbHb78ex3+ome01nqP2U0kfklQF7Ce63RWSPinpsZQ2rZH075K6p23nPkn7JQ2V9EAc2464/eVpdSXpRkmvxftzYWob0uoOiGOrjT9/uaQvSlITbbhUUlW83SWSzozrXCJpcRzjcknnt3J/f0jSXyS9IWmPpFWSbm3FekdKul/Stvgzl0i6Mq3OuDjuf5H0OUmr47r/J+ncDNvM6rdjlu+6Jh2AmeXcvcBngY8RnXuW6m+BlSGEv6SUfQqoA34M7ABOAK4GJgNnZdj+fwOvAv8G9M6wvNFlQAUwG9hCdKPza4B3S5qaoYfw+8A+4D+AQcCX4rac0VhB0mTgaaL/tM4GVgEjgQ8Q3SR+j6SewJNEN16fTXQPzuOBLwKTgJnNxIyks4F5wGvAN+LP+iwwX9L0EMKzwG+Jhsu/DzwAPAg0NLfd2BnAh4DbgNuJ7u8KcD3wIvAwsAc4HfiXuG0fT9tGlzi+auAfgVOI9usW4JaUel8DbgYeBf6XaP8/TNRDuSalveVE++sY4CfAy8AlwH8Bo4Ab0j7/VGBWHP8+4J+AuZI+Ha9zO7A3Ln9A0hEhhDea2iFxcvgb4HHgX4FDwFjggqbWidcbDCwk+o39GKghOsXhbkn9Qwg/SFvlI0S9tLcT/d4/E8c9PYSwKN5mVr8ds4KQ9M2C/fDDj9w/iG6Yviqt7HQgADenlffKsP5Vcd13p5R9Iy77XYb6V8fLRraw3elxvcsyrDsf6JJS/qW4fGJK2VNEScWEDNtuvPf4LUQJybFpyz8Xb+/sFvbdEmA7MCilbASwC1iYUjY60/5sYptd47oNZLiRexP76itAPTAspey+eDvfSKs7B6hNeT+YKLl5NG2fNu6Dx1LKbojLPpW6L4Hfx/FOTGtDXer+J0oCA1EP5rgM5Ve3sG9+THR6Qlkr9t/NKWU/iMtmpJR1B/5MlDxXxGXjUuI+KqXuEOBN4OmUsqx+O374UQgPD/WaFaf7gHGSTkkpu5zoj9cvUyuGEPbCW0ODh0saCDwTL35Xhm3/pDUBNG433nafeLsvEiVQmbZ7Rwghtddsfvw8Nt7GEKIes1+EDJMoQgghfnkZsADYImlg4wNovITJOU3FrGioeipwdwhha8q2NwL3A6dKGtBMs1uyMISwJEPsjd9BF0kVcbxPEfXunZBhO7envZ8PDNXb5wxeAHQDfpS2T+8k2v+pZhL1Xv4iJZ4A/CdRAnhJWv0n0vb/osYYQgivZCgfmyH+VDuBPkCrhoVTzAQWhxD+lBJ3HVFC2It3fs+VIYTVKXU3A78CTpdUERe3+7djViic+JkVp/vi58sB4nPFPgw8E0JYm1pR0Tl7lcBuoj/CW4mGUCEaRku3OkPZO0gaKemXkt4g6lnZGj/6NLHd9Wnvd8TP/ePncfHz0hY+egIwI+XzGh8r4uWDm1l3dPz8coZlL6XVaY+M+07SaZKeIOrN3EEUb2NCk76vDoYQatLKGvdVv/j5yPh5RWqlODFal7buaKLe4fq08sb2jkkrfzXt/c4WyvvRvFvjOB+Kz3G8X9Jlklo6FelImv+e0uNekV4xLhNv769sfjtmBcHn+JkVoRDCK5IWAZdJugG4mCiBuje1XtzTMZ9oaOwWoiHivURDZnPJ/J/DfS19vqQyol6SgcC3if4Y7yHqcfxtE9tNTzze2lzac2iiXmr9x4FvNrF8YwvrN7fd1nx+c96x7yQdRZTkrSI6l+xVomHTUcD/4537qrlzCVuzr5ShrDnp22jqe2rp+8u88RA2SzqBKOG6gKjn76PA85LODCG0+Htr4vPS427Nvuio345oHWX4AAADmUlEQVRZ3nDiZ1a87iUaEryAqOfvAFHSlWoGUXL2vhDCgsZCScdm+dnHE00muDyE8NbQsqKZt4e3c5uNw4hTWqi3GugTQmjP3SnWxc9HZ1jWWJbeM5mt9wHlwEXxkDIAki7OYpvr4uejebv3trHn90hgc1rdKZLK0nr9jknbVocJ0aWFHokfSLoe+BHwN6SdmpBiPc1/T+uaKE81gSghbPxOs/ntmBUED/WaFa/fEJ3Q/nmi86EqQwjp15xr7D1KPxb8Y5af3dx229rjBEAIYRPRjN6rFF2P8K+kXHrk18BJki7NUKdcUp9mPmMD0eSOK+NzuxrXG040S3phCGF7e+Jvxjv2Vdxj+sUstvlHotmx10tK/Q6uJhpqT1VJNIv6rcugxPuycXLN3CziaFET50wujp8znRLQqBI4IZ6F3bitbsDfEfVaP55Wf1bcu9pYdwhRz+LClH8X7f7tmBUK9/iZFakQwuuSHiLqUYK0Yd7Y00QzKn8p6cdEQ4yziHoBs7GMqIfuB4quxbeNaEbvybx9Plp7XEcU83OSZgMrgeFEl3O5GNgAfIdoQsLvJd0H/IVo6Hoi0XmO7+XtySuZfJHocinPSvoZUaJ6LdFkiX/IIvamPAx8C3g4blNXokkG3dq7wXj49D+BG4FHJD1I3AMLrE2rPpsoIZwdD7muINqXFwM/DCFkOjcul34Rn3LwONEw92Ciy+fsJkrumvJNou+zMv7t1hDtt1OAv8/wn5yXgKcl3QYcJLqcSznwzyl1sv3tmOU9J35mxe1eosRvO1GC8VdCCNviIcXvAl8mSvweIur9qW3vh4YQ6iTNJLrOXWOy9CTRrMins9hutaKLEH8V+ARwGNEf/HnECWUIYa+k6UR/0C8j6tXZTXTtuh8RJaXNfcYT8YV9v0p0rcIAPEt0CZpn2xt7M5/3sqT3EV0u51vAG0TXSvw5Ue9je91EdF7ltUSzoRcDFxFday/18/fF++vfiZKb/kTJ4ZeA72Xx+a11D9GFrK+JP3sb0czar4cQ0ieMvCWEsEXSaUT77Bqia0quAD4eQrgnwyq/Jppk9A9E10d8CZiZeopDtr8ds0LQeN0rMzOzoiNpHNF5jjeGEL6VdDxmSfM5fmZmZmYlwomfmZmZWYlw4mdmZmZWInyOn5mZmVmJcI+fmZmZWYlw4mdmZmZWIpz4mZmZmZUIJ35mZmZmJcKJn5mZmVmJcOJnZmZmViL+P0ePA6WMdjJPAAAAAElFTkSuQmCC\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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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",
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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": [
      "/home/travis/build/statsmodels/statsmodels/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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yV5f3G8c/3ZLIDIayEvfeeIm5QEVTQirtubW0dte6Bs1r156h1tVqtVkUBq4iIxQEqS6bsEXZYSdiEkHX//sjBIiUhh5yT55yc6/16nVc9z3MSLirIxfM89/c25xwiIiIiErl8XgcQERERkfJRoRMRERGJcCp0IiIiIhFOhU5EREQkwqnQiYiIiEQ4FToRERGRCBfrdQCv1a1b1zVr1szrGCIiIiLHNHfu3CznXMqRx6O+0DVr1ow5c+Z4HUNERETkmMxs/dGO65ariIiISIRToRMRERGJcCp0IiIiIhFOhU5EREQkwqnQiYiIiEQ4FToRERGRCKdCJyIiIhLhVOhEREREIpwKnYiIiEiEU6ETERERiXAqdCIiIiIRToVOREREJMKp0ImIiIhEOBU6ERERkQinQiciIiIS4VToRERERCKcCp2IhBXnHAWFRV7HEBGJKLFeBzgaM7sQGA20B/o45+b4j58BPAnEA3nAH51zX/vPfQs0BA74v81g59z2ik0uIsdjT24+01dnM21VJtNWZrJtTy6ntavPyJ5pnNw2hbgY/d1TRKQ0YVnogMXACOC1I45nAcOcc5vNrBMwGUg97Pylh8qfiISvwiLH4ozdTFuZybRVmczbsIvCIke1+Bj6t6zLyW1TmLRoK18s2UqdavEM79qIkT3S6JRaEzPzOr6ISNgJy0LnnFsG/M9/uJ1z8w97uwRINLME59zBCownIsdh255cpq7M5LtVWXy/KpOdOfkAdEqtyQ2DWjCoTQo9mtQmPrb4atxDwzoybWUm4+dl8N6sDbw1fR1t6ldnZI80zuueSv2aiV7+dEREwoo557zOUCL/bdQ7jnbVzcwuAG50zp1+2GeTgUJgHPCYK+EnZ2bXA9cDNGnSpOf69etDkl8kmuXmF/Ljuh3FV+FWZrFi214A6lZPYFCbugxqncLA1nWpWz3hmN9rd04+E37azLh5m5i/YRc+g4GtUxjZI5XBHRpQJT4m1D8dEZGwYGZznXO9/ue4V4XOzKYADY5y6j7n3Cf+z3zLUQqdmXUEPqX4Obl0/7FU51yGmdWguNC965z757Fy9OrVy82Zo7u0IuXlnCM9cx9TV2YxbWUms9Zmk5tfRHyMj17NajOoTQqDWqfQvmGNct02XZO5j/HzMvh4fgYZuw5QPSGWoZ0bMqJHKn2a19EtWRGp1MKu0JXF0QqdmaUBXwNXOed+KOHrfg30cs7dfKwfQ4VO5Pjtzsnn+9XFBe67VZls3p0LQIu61YoLXJu69GuRTNX44D/dUVTkmLk2m/HzMvh80RZy8gppXKcKI7qnMaJHKk2TqwX9xxQR8VqlKHRmlgRMBR5xzo077HOxQJJzLsvM4oD3gSnOuVeP9WOo0ImUXUFhEQs3/Xcxw8KNuyhyUCMhlhNa1WVQmxRObF2XxnWqVmiunLwCvli8lfHzMvghPQvnoHez2ozokcbQLg2pmRhXoXlEREIlogqdmZ0P/AVIAXYBC5xzQ8zsfuAeYNVhHx8M7AemAXFADDAFuN05V3isH0uFTqR0m3cd+LnAfb8qiz25BZhBl7QkTmpdXOK6NU4iNkxGi2zedYCP52cwbt4m1mTuJyHWx+CODRjRI5UTW9UNm5wiIscjogpdRVKhE/mlA3mFzFyb7V/MkEl65n4A6tdMYFDrFAa1SWFgq7rUrhbvcdLSOedYuGk34+dt4tOFm9mVk09KjQTO757KiB6ptGtQ0+uIIiIBU6ErgQqdRDvnHCu27f15NersdTvIKygiIdZHn+Z1OKlNcYlrXa96xC44OFhQyDfLtzNuXgbfLN9OQZGjY6OajOyRxvBujcq00lZEJByo0JVAhU6i0e4D+Uz1X4H7blUm2/YUj3JsXa+6fzFDCn2b1yExrvKNA8ned5AJCzczbl4GizJ2E+szTm6bwogeaZzWvh4JsZXv5ywilYcKXQlU6CTaLNuyhyvfnM32vQepVSWOga3qMqhNXU5snUKjpCpex6tQK7ftZdy8Tfx7fgbb9hT//zGsa0NG9Eije+OkiL0iKSKVlwpdCVToJJrMWpPNtf+cQ7X4WJ4f1Y3ezeoQ41NpKSxyfL86i/HzNjF5yVZy84toUbcaI3qkcn6PNFKjrOiKSPhSoSuBCp1Eiy+XbOXm9+eTVrsK71zTVyWlBHtz85m0aCtj521i9todmEH/FsmM6JHGWZ0aUC0hLHdMFJEooUJXAhU6iQYf/riRu8f/ROe0JP7x697UCfMVquFiQ3YOH8/PYPz8TazPzqFqfAxndmrAyB5p9G+RjE9XN0WkgqnQlUCFTioz5xyvTl3DU18s58TWdXn1sp66wnQcnHPMXb+TcfM28dnCLew9WECjWomc3yOVET3SaJlS3euIIhIlVOhKoEInlVVRkePxz5fxxvdrGd61Ec9c2JX4WA3VLa/c/EL+s3Qb4+ZtYtrKTIocdGucxMgeqQzr2oikqrr6KSKho0JXAhU6qYzyC4u4c+xPfDw/g18PaMaD53TQ7cEQ2L4nl08WbGbcvE0s37qX+Bgfp7arx8ieaZzcNoU47UohIkGmQlcCFTqpbHLyCvjNv+bx7YpM7hjcht+e0krjN0LMOcfSLXsYNzeDTxZkkL0/j+Rq8Qzv1oiRPdLo2Kim/h2ISFCo0JVAhU4qk105eVz11o8s3LiLx8/vzMV9mngdKerkFxYxbWUm4+ZtYsrS7eQVFtG2fo3iESjdU6lXM9HriCISwVToSqBCJ5XFlt0HuOKN2azfkcOLo7pxZqeGXkeKertz8pnwU/Et2fkbduEzOLF1CiN6pDKkY4NKuROHiISWCl0JVOikMli9fR9XvDGLPbkF/O2KXvRvmex1JDnCmsx9jJ+XwcfzM8jYdYAaCbGc3bkhI3um0btZbd2SFZEyUaErgQqdRLr5G3Zy9Vs/EuPz8dZVvemUWsvrSFKKoiLHzLXZjJubwaTFW8jJK6RxnSqM6J7GBT3TaFynqtcRRSSMqdCVQIVOItnUlZnc+M5cUmok8M41fWiaXM3rSBKAnLwCvli8lfHzMvghPYs4n4+XL+3B6R3qex1NRMJUSYVOa+pFItQnCzK45q0faVa3GmNv6q8yF4Gqxscyokca717bl+/vOpV2DWtw47tzmbRoi9fRRCTCqNCJRKB//LCWWz5YQI+mtRlzQz/q1dDKyUiXmlSFd6/tS5e0Wtz8/nw+WZDhdSQRiSAqdCIRxDnHM5NX8PCEpQzuUJ9/Xt2HmolxXseSIKmZGMc/r+lLr6a1uW3MAsbO3eR1JBGJECp0IhGisMhx78eLeOmb1Yzq3ZiXL+2hsReVUPWEWN66qg8DWtblj2MX8t6sDV5HEpEIoEInEgFy8wv5zb/m8v7sjdx8Siv+NKIzsdpWqtKqEh/D36/sxcltUrj340W8PX2d15FEJMzpTwSRMLcnN58r35zN5CXbeGhYB+4Y0lYzy6JAYlwMr17ek8Ed6vPQp0v427Q1XkcSkTCmQicSxrbvzWXUazOZu34nL4zqxlUnNPc6klSghNgY/nppD4Z2acjjny/jpa9XeR1JRMJUrNcBROTo1mfv5/I3ZpO592Dx7be29byOJB6Ii/HxwkXdSIjx8cyXKzlYUMTtZ7TRVVoR+QUVOpEwtGTzbq5880cKiop477q+dG9S2+tI4qHYGB9PX9iVuBgff/l6NXkFRdx9VjuVOhH5mQqdSJiZuSab696eQ43EWD64vj+t6tXwOpKEgRif8acRnYmP9fHatDUcLCjioWEdVOpEBFChEwkrXyzeyu8/mE+TOlX559V9aJRUxetIEkZ8PuORczsSH+vjje/XkldYxGPndsLnU6kTiXYqdCJh4oPZG7j340V0bZzEm1f2pna1eK8jSRgyM+4f2p74WB+vfJtOXkERT43sQoxKnUhUU6ET8Zhzjpe/TefpySs4qU0Kr1zWg6rx+q0pJTMz7hzSloRYH89PWUV+YRHPXthVswlFopj+1BDxUFGR45HPlvLW9HWc163Rzw++ixyLmXHr6W2Ij/Xx5y9WkFdQxAujuhMfq18/ItFIhU7EI3kFRdzx0UI+XbiZawY2576z2+tZKAnYb05uRXyMj8cmLiP/X3P566U9SIjVlnAi0UZ/lRPxwP6DBVz7zzl8unAzd53ZjvuHqszJ8bv2xBY8em5HpizbznX/nEtufqHXkUSkgqnQiVSwHfvzuOTvs/h+VSZPjezMTSe31OgJKbfL+zfjqZGd+W5VJle/9SM5eQVeRxKRCqRCJ1KBMnYd4MJXp7Nsyx5evawnF/Vu4nUkqUQu6t2E//tVV2auyebKN2ezNzff60giUkFU6EQqyKpte7nglels33uQd67uw+CODbyOJJXQ+d3TePHi7szbsIvL35jN7gMqdSLRQIVOpALMXb+TC1+bQUGRY8z1/enbItnrSFKJndOlES9f2oMlm3dz6d9nsnN/nteRRCTEVOhEQuybFdu57O+zSKoSx7gbB9ChUU2vI0kUGNKxAa9f3ouV2/Zx8d9mkrXvoNeRRCSEVOhEQujf8zO47u05tEipxkc3DqBJclWvI0kUOaVdPd64shfrsvcz6vWZbN+T63UkEQkRFTqREHnj+7XcOmYBvZvV4YPr+5FSI8HrSBKFTmydwltX9WHzrgNc9PpMtuw+4HUkEQkBFTqREHh+ykoe/WwpZ3VqwD+u6k2NxDivI0kU69cimXeu6UPW3oP86rUZbNyR43UkEQkyFTqRIFu4cRcvfLWKEd1TeemSHiTGaWq/eK9n0zq8e21fdufkc9FrM1iXtd/rSCISRCp0IkFUVOR46NMl1K2ewMPndiRGuz9IGOnaOIn3r+/HgfxCfvXaDFZv3+d1JBEJEhU6kSAaN28TCzbu4u4z2+k2q4Sljo1q8cH1/SlyMOr1GazYutfrSCISBCp0IkGyJzefp75YTvcmSZzfPdXrOCIlatugBmNu6EeMzxj1+gwWZ+z2OpKIlJMKnUiQvDBlFdn783hkeCd8utUqYa5lSnU+vKE/VeNjueRvM1mwcZfXkUSkHFToRIJg1ba9vD19HaN6N6ZzWi2v44iUSdPkaoy5oR9JVeO57O+zmLNuh9eRROQ4qdCJlJNzjtETllA1PoY7Brf1Oo5IQNJqV2XMDf2oVyOBK96czYz0bK8jichxUKETKafJS7byw+ps/jC4LcnVNTxYIk/DWlX44Pp+pCZV4df/mM20lZleRxKRAKnQiZTDgbxCHv1sGe0a1ODSvk28jiNy3OrVTOSD6/vRIqU61749h6+Xb/M6kogEQIVOpBxenZpOxq4DjB7ekdgY/XaSyJZcPYH3r+tL2wY1uOGduXyxeKvXkUSkjPQnkMhx2rgjh1enpnNOl4b0a5HsdRyRoEiqGs+/rutLp9Ra/Pa9eUxYuNnrSCJSBip0IsfpsYlL8Zlx39D2XkcRCaqaiXG8c01fejapzS0fzGfc3E1eRxKRY1ChEzkO363KZPKSbdx8aisa1qridRyRoKueEMtbV/emf8tk7hi7kA9mb/A6koiUQoVOJED5hUU8PGEpTZOrcu2Jzb2OIxIyVeNjeePK3gxqncLd4xfxzxnrvI4kIiVQoRMJ0NvT17F6+z4ePKcDCbExXscRCanEuBhev6Inp7evz4OfLOHv363xOpKIHIUKnUgAtu/N5fkpqzi5bQqntqvndRyRCpEQG8Mrl/VgaOeGPDZxGX/9ZrXXkUTkCLFeBxCJJE9NWsHBgkIePKcDZtqvVaJHXIyPF0Z1Iy7GeHryCg4WFHHb6a31+0AkTKjQiZTR3PU7GTdvEzee1JIWKdW9jiNS4WJjfDz7q27Exfh48atV5BUUcdeZbVXqRMKACp1IGRQVOUZ/uoT6NRP43amtvI4j4pkYn/HUyC7Ex/p4dWq6rliLhAkVOpEy+HDORhZl7OaFUd2olqDfNhLdfD7jsfM6ER/r4x8/rCOvoIhHz+2Ez6dSJ+IV/ckkcgy7c/L58+QV9G5Wm+FdG3kdRyQsmNnPK71fnZpOfmERfxrRhRiVOhFPqNCJHMNzU1ayKyeP0cP76LaSyGHMjLvObEt87H+fqXvmwq7a11jEAyp0IqVYvnUP78xczyV9m9CxUS2v44iEHTPj9jPakBDr4+nJK8grLOKFUd2JU6kTqVAqdCIlcM7x0CdLqJEYyx/OaOt1HJGw9ttTWpEQ6+OxicvIK5jHXy/trsHbIhVIf4USKcHERVuYtXYHdwxuS+1q8V7HEQl7157YgkfO7ciUZdu44Z255OYXeh1JJGqo0IkcRU5eAY9PXEbHRjW5uE8Tr+OIRIwr+jfjTyM6M3VlJte8/SM5eQVeRxKJCip0Ikfx8jfpbNmdy8PDO2rVnkiALu7ThGcu6MqM9Gx+/eaP7DuoUicSaip0IkdYn72f16et4bxujejVrI7XcUQi0sieabwwqjtzN+zk8jdmsftAvteRRCo1FTqRIzz62VLiYox7zm7vdRSRiDasayP+ekkPFmfs5rK/z2JXTp7XkUQqLRU6kcN8s2I7U5Zt53entaZ+zUSv44hEvDM7NeC1y3uyYtteRr0+k+x9B72OJFIplVjozKzIzAoDfQUjlJldaGZL/Bl6HXb8DDOba2aL/P976mHn4s3sdTNbaWbLzWxkMLJI9MgrKOLRCUtpUbcaV5/Q3Os4IpXGqe3q8/crerE2az+jXp/J9j25XkcSqXRKm0P3COCOOHYe0AmYDKwADGgLDAYWAZ8EKddiYATw2hHHs4BhzrnNZnYoR6r/3H3AdudcGzPzAXr4SQLy5g9rWZO1n7eu6k18rC5eiwTToDYpvHVVH655+0dGvT6Tf13Xl4a1qngdS6TSKLHQOedGH/7ezK4GGgCdnHMrjjjXHvgG2BCMUM65Zf7ve+Tx+Ye9XQIkmlmCc+4gcDXQzv+5IorLn0iZbNuTy1++WsXp7etzctt6XscRqZT6t0zmn1f34df/+JFfvTaD967tR+M6Vb2OJVIpBHIZ4k7gpSPLHPxcwP4K3BWsYGUwEpjvnDtoZkn+Y4+a2Twz+8jM6ldgFolwf/p8GflFjgfO0UIIkVDq1awO717bl905+Yx6fSbrs/d7HUmkUgik0DUFDpRyPsf/mTIxsylmtvgor3PL8LUdgaeAG/yHYoE04AfnXA9gBvBMKV9/vZnNMbM5mZmZZY0sldSP63bw7wWbuf7EFjRNruZ1HJFKr1vjJN67rh85eQVc9JoWSogEQyCFbiVw/WFXw35mZrWB6yl+rq5MnHOnO+c6HeVV6nN4ZpYGfAxc4ZxL9x/OprhQfux//xHQo5Qf+3XnXC/nXK+UlJSyRpZKqLCoeL/WRrUS+c0pLb2OIxI1OqXW4p1r+rJjfx53jVuEc0c+si0igQik0N0LtARWmdn/mdmNZnaDmT1HcdlrQfHChJDxl8mJwD3OuR8OHXfF/yWYAJzsP3QasDSUWaRyeH/2BpZu2cO9Q9tTNb60NUIiEmydUmtx11ntmLJsG+/NDsoj2CJRq8yFzjk3ERhC8cKHW4GXgVeAW/zHzvJ/ptzM7Hwz2wT0Byaa2WT/qZuBVsADZrbA/zr0BPtdwGgz+wm4HPhDMLJI5bVzfx7PfLmC/i2SGdq5oddxRKLSVQOacWLrujz62VJWb9/ndRyRiGXHc5nbv+CgGcVjS9Y657YFOVeF6dWrl5szZ47XMcQD9/97Ee/P3sjE3w+kXYOaXscRiVrb9+Qy5PlpNEqqwse/OUFjg0RKYWZznXO9jjx+XL9rnHPbnHOznHMzI7nMSfRasnk3783awOX9mqrMiXisXs1EnhrZhSWb9/Dsf8r8KLaIHCagQmdmSWb2pH816n4z2+f/5yeOtlhCJBw55xj96RKSqsZz2+ltvI4jIsDgjg24pG8TXp+2humrNUZUJFBlLnRmlgrMp3geHRQvTphE8W4SdwPzzKxR0BOKBNmnCzfz47qd3DmkLbWqxnkdR0T87h/anuZ1q3H7hwvZlZPndRyRiBLIFbo/AfWBc/zjRX7lnLvQOdcZGOo/96dQhBQJlv0HC3ji82V0SavFr3o19jqOiBymanwsL47qTvb+g9wzXqNMRAIRSKE7E3jBOff5kSecc5OAvwBnBSuYSCj85evVbNtzkNHDO+Lz2bG/QEQqVKfUWvxhcFsmLd7KR3M3eR1HJGIEUuhqABmlnN/k/4xIWFqTuY83vl/DyB5p9GhS2+s4IlKC609sQf8WyYz+dAnrsrQ1mEhZBFLoVgAXmNn/fI2ZxQAXEMBOESIVyTnHI58tJSE2hrvOaut1HBEphc9nPPurrsTF+LhlzALyC4u8jiQS9gIpdC8Cg4CvzexcM2vnf50HTAFOBJ4PRUiR8vpq2Xa+XZHJrae3pl6NRK/jiMgxNEqqwhPnd2bhxl28+NUqr+OIhL0y73XknHvTvyvDQ8D4w04ZcBC41zn3VnDjiZRfbn4hj3y2lFb1qnPlgGZexxGRMhrapSHfrkjjr9+s5sTWKfRpXsfrSCJhK6A5dM65J4E04FKK93a9F7gYSHXOPRX8eCLl98b3a9mwI4fRwzoSF6MJ9CKR5KHhHWlcpyq3jVnA7gP5XscRCVsB/+nmnMt2zn3gnHvK/xrjnNsRinAi5bV51wFe+no1Z3ZswMDWdb2OIyIBqp4Qy/MXdWPrnlwe/GSx13FEwlbAhc7MzjSzl8xsopl95v/nwaEIJ1JeT3y+jCLnuG9oe6+jiMhx6t6kNree1ppPFmzm3/NLG7YgEr0C2Ski3sw+oXiHiN8AvYG+/n+eZGb/NrP40MQUCdyM9Gw++2kLN57UksZ1qnodR0TK4TentKJ3s9o88O/FbNyR43UckbATyBW6h4BhwLNAPedcPedcCpACPAMMBx4IfkSRwBUUFvHwhCWkJlXhppNbeh1HRMopxmf836+6AXDbmAUUaJSJyC8EUuguAd51zt3pnPt552T/M3V3Ae8ClwU7oMjx+NesDSzfupcHzmlPYlyM13FEJAga16nKY+d3Ys76nbzybbrXcUTCSiCFrhEwvZTzM4CG5YsjUn7Z+w7y7JcrGNiqLkM6NvA6jogE0bndUjmvWyOe/2oV8zfs9DqOSNgIpNBtBvqVcr4vsKV8cUTK75kvV5CTV8jo4R0w036tIpXNI+d1okHNRG4ds4B9Bwu8jiMSFgIpdO8Dl5vZY2b280aYZlbbzB4FLgfeC3ZAkUD8tGkXH/y4kSsHNKNVPW0tLFIZ1UyM4/lR3di4I4eHP13idRyRsBBIoXsYmEzxMOEsM9tqZluBLOA+/7lHgh9RpGyKihwPfbqE5GoJ3HJ6a6/jiEgI9W5Wh5tPacVHczcx8SfdHBIJZOuvg8DZZnYOMBRo5j+1DpjgnPs86OlEAjB+fgbzN+zi6Qu6UDMxzus4IhJivzutNdNWZXHP+J/o3iSJRklVvI4k4pnj2SniM+fcTc65s/yvm1TmxGt7c/N5ctJyujdJYmSPNK/jiEgFiIvx8cKobhQWOW7/cAGFRc7rSCKe0caWUim8+NUqsvcf5OHhHfH5tBBCJFo0Ta7G6OEdmblmB3/7bo3XcUQ8U+ZbrgBmdgZwLdACqAMc+Senc85piqtUqNXb9/KPH9ZxUa/GdElL8jqOiFSwC3qm8e2KTJ79cgUntKxL57RaXkcSqXCBbP11G/AFcBKQAUwDph7xmhaCjCIlcs7x8ISlVImP4Y4hbb2OIyIeMDMeP78TdasncMuY+eTkaZSJRJ9ArtDdRnFpO9M5lxeiPCIBmbxkG9+tyuKhYR2oWz3B6zgi4pGkqvE8+6uuXPr3WTw2cRlPnN/Z60giFSqQZ+jqAmNU5iRc5OYX8tjEpbStX4PL+zX1Oo6IeGxAy7rcMKgl783awJdLtnodR6RCBVLo5lH87JxIWHht6ho27TzA6OEdiY3R+h4RgdvPaEOn1JrcNe4ntu/J9TqOSIUJ5E/BWyneKeKMUIURKatNO3N4+dvVDO3SkP4tk72OIyJhIj7WxwujunMgv5A/fLSQIo0ykShR4jN0ZvblUQ7vBb4ws3UUDxQuPOK8c84NCVo6kRI8PnEZZnDf2e29jiIiYaZlSnUeOKcD9328mH9MX8c1A5t7HUkk5EpbFNEGONpfbTZQfGVPt1/FEz+szmLS4q384Yw2mgwvIkd1SZ8mfLM8k6cmLWdAy2TaN6zpdSSRkCqx0DnnmlVgDpEye+4/K0mrXYXrBunvFCJydGbGUyM7c+YL33HLB/P59OaBJMbFeB1LJGT0JLlElLVZ+5mzfieX9Wuq/ziLSKmSqyfwzIVdWbltH09OWu51HJGQUqGTiDJu7iZ8Bud3T/U6iohEgJPapHD1Cc15a/o6vlm+3es4IiFTYqEzsyIzKzCz+MPeFx7jpfHcEjKFRY5x8zYxqE0K9Wsmeh1HRCLEnWe2pV2DGvxx7EKy9h30Oo5ISJS2KOIRihdFFBzxXsQTM9Kz2bI7l/uGamWriJRdYlwML4zqzrCXvufOsT/xxpW9MDtyK3KRyFbaoojRpb0XqWhj526kZmIsp7ev73UUEYkwbRvU4N6z2jF6wlLenbmey/s38zqSSFDpGTqJCHty8/liyVaGd2ukxRAiclyuHNCMk9qk8NjEZazattfrOCJBVdpg4UHH8w2dc9OOP47I0X3+0xZy84u4oGdjr6OISIQyM56+sAtnPf8dv/9gAf/+7QASYvUXRKkcSnuG7lsCe2bO/J/X7w4JurFzN9GqXnW6ptXyOoqIRLB6NRL58wVduObtOTwzeQX3De3gdSSRoCit0J1SYSlESnFo9tzdZ7XTg8wiUm6nta/P5f2a8rfv1nJSm3oMbF3X60gi5VbaooipFRlEpCSaPSciwXbv2e2ZsSab2z9cwORbB1G7WrzXkUTK5bgWRZhZazM7wcx0/0tCSrPnRCQUqsTH8MKobuzMyePu8bo+4MIAACAASURBVD/hnKZySWQLqNCZ2UVmth5YDkwDevqP1zWzVWZ2YQgyShQ7NHtuZI80r6OISCXTsVEt7hzSjslLtjHmx41exxEplzIXOjM7F3gf2AA8QPEiCACcc1nAMuDyYAeU6DZ27kZqJMZyRgfNnhOR4LtmYHNOaJXMwxOWsiZzn9dxRI5bIFfo7gemOedOBF47yvlZQNegpBLhsNlzXTV7TkRCw+cznr2wGwlxPm75YAF5BUVeRxI5LoEUuo7Ah6Wc3wroMooEzX9nz+l2q4iEToNaiTw5oguLMnbz/JSVXscROS6BFLpcoLSn0psBu8qVRuQwY+duomVKNbo1TvI6iohUcmd2asCo3o15ZWo6M9dkex1HJGCBFLrvgYuPdsK/2vVq4OtghBI5NHvugp6NNXtORCrEA+d0oFlyNW4fs4DdOflexxEJSCCFbjTQ0cy+AUb4j/Uys5uBBUBN4NHgxpNopdlzIlLRqiXE8vxF3di+9yD3/nuRRplIRClzoXPOzQOGAA3476KIJ4EXgTxgiHNuWdATStQ5NHvuxNYpNKil2XMiUnG6Nk7itjPaMPGnLYyfl+F1HJEyK23rr//hnPsOaG9mXYE2FBfC1cA8p7/KSJAcmj1379ntvY4iIlHoxpNaMnVlJg9+sphezWrTNLma15FEjimQOXQ/LzV0zi10zn3knBvjnJt7qMyZ2UmhCCnRRbPnRMRLMT7juYu64fMZt41ZQEGhRplI+AvkGbopZlbiDsZmNhiYWP5IEs00e05EwkFqUhWeOL8z8zbs4qVvVnsdR+SYAil0McB/jrZ/q5kNAz4FZgYrmEQnzZ4TkXAxrGsjRvRI5cWvVjF3/Q6v44iUKpBCdzqQDEwys58fKDCzC4BxwDfAOcGNJ9FGs+dEJJw8PLwjqbWrcOuYBezN1SgTCV+BrHJdD5wBtAA+MbMEM7uc4v1dPwfOdc7lhiamRAPNnhORcFMjMY7nL+rO5l25PPTpEq/jiJQokCt0OOdWUDy6pAfFt1f/QfHVuZHOubzgx5NootlzIhKOejatze9ObcX4eRl8unCz13FEjiqgQgfFK1yBs4CWwLvAxc65wmAHk+ii2XMiEs5uPqUVPZokcd/Hi8jYdcDrOCL/o8RCZ2b5ZpZ3tBfwHVANuAQ4eNi5gxUVXCqXQ7PntBhCRMJRbIyP5y/qjnNw25gFFBZp9KqEl9IGC/8L0K9YqRCaPSci4a5JclUeObcjt3+4kFenpvPbU1p5HUnkZyUWOufcryswh0SxQ7PnRvZI0+w5EQlr53dP5ZsVmTz3n5UMbFWXrlqRL2Ei4GfoRIJNs+dEJFKYGY+d14l6NRK4dcwC9h8s8DqSCFDKFTozGwTgnJt2+PtjOfR5kbLS7DkRiSS1qsTxfxd14+K/zeTRz5by5MguXkcSKfUZum8BZ2ZV/CNJvqX0Z+rMf173zKTMDs2eu+vMdpo9JyIRo1+LZG46qSUvf5vOyW1TOLNTQ68jSZQrrdCdAnDYfLlTQh9Hoo1mz4lIpLr19DZ8vzqLu8cvolvj2hq5JJ4qbVHE1NLei5SXZs+JSCSLj/Xx/EXdGPri9/zhowW8c3VffD7daRBvaFGEeEaz50Qk0rVIqc5Dwzrww+ps3vh+rddxJIqVtijiweP4fs4592g58kgU0ew5EakMLurdmG9WbOfPk5czoFUyHRvV8jqSRKHSnqEbfRzfzwEqdHJMmj0nIpWFmfHkiC6c+cI0bvlgARNuHkiVeP13TSpWibdcnXO+43jpV7CUiWbPiUhlUrtaPM9e2I3V2/fxxOfLvI4jUSgsn6EzswvNbImZFZlZr8OOn2Fmc81skf9/T/Ufr2FmCw57ZZnZ8979DORYNHtORCqbga3rct2JzXln5nq+WrbN6zgSZcKy0AGLgRHAkUOKs4BhzrnOwJXAOwDOub3OuW6HXsB6YHxFBpayOzR77oKejTV7TkQqlTuGtKV9w5rcOfYntu/N9TqORJGwLHTOuWXOuRVHOT7fObfZ/3YJkGhmCYd/xsxaA/WA70KfVI6HZs+JSGWVEBvDi6O6se9gAX/86CecK20ev0jwhGWhK6ORwHzn3MEjjl8MjHGl/C4ys+vNbI6ZzcnMzAxpSPklzZ4Tkcqudf0a3D+0PVNXZvL29HVex5Eo4VmhM7MpZrb4KK9zy/C1HYGngBuOcnoU8H5pX++ce90518s51yslJeX4fgJyXDR7TkSiwWX9mnJqu3o8MWk5K7bu9TqORAHPCp1z7nTnXKejvD4p7evMLA34GLjCOZd+xLmuQKxzbm4Io0s5aPaciEQDM+PPF3ShZmIst3wwn9z8Qq8jSSUXUbdczSwJmAjc45z74SgfuZhjXJ0T7+z1z54b3rWRZs+JSKVXt3oCT1/YleVb9/LnL/7nsXCRoCptsPAvlGHnCAfkApuAqYctXgiYmZ0P/AVIASaa2QLn3BDgZqAV8ICZPeD/+GDn3Hb/P/8KOPt4f1wJrc8XafaciESXU9rW49cDmvHmD2s5qW0KJ7XRYz4SGlbWFThmVkRxaQM4ctbEkccLgVeB35e2OCEc9OrVy82ZM8frGFHhwlens2N/HlNuP0njSkQkauTmFzL8pe/ZmZPPF7ecSHL1hGN/kUgJzGyuc67XkccDueWaBiwE3gV6AbX8rz7Av4AFQGugJ/AB8BvgrvLFlspiXdZ+flyn2XMiEn0S42J4YVR3dufkc9e4RRplIiERSKH7C5DunLvSOTfPP8x3r3NujnPuCmAt8IR/VtzlwFfAr0OQWSLQuHmaPSci0at9w5rcdVY7pizbxicLjvuJJJESBVLoTqe4pJXkK2DIYe8/A5odRyapZIqKHOPmavaciES3qwY0o0taLf40aRn7DxZ4HUcqmUAKXRHQtZTzXfnvs3SH5AScSCqdGWuy2azZcyIS5Xw+46FhHdm25yCvfJt+7C8QCUAghe7fwHVmdreZ1Th00MxqmNk9wLX+zxzSH1gZnJgSycbO3aTZcyIiQM+mtTm/eyqvf7eGjTt0zUOCJ5BCdxswA3gC2GlmW8xsC7ATeByY5f8MZpYI7AL+L7hxJdLszc1n0uItmj0nIuJ315ntiPUZj09c5nUUqUTKPIfOObfLzAYB5wFnAU38p9YDk4BPDo0occ7lAjcFOatEIM2eExH5pQa1EvntKa14evIKpq/OYkCrul5HkkqgzIUOwF/YPva/RI5p7NxNtEypRrfGSV5HEREJG9cMbM77szfw8ISlTPz9QGJjImrjJglDAf8KMrPaZnaBmf3RzO4ws5H+LblEfkGz50REji4xLob7h7Znxba9vD97g9dxpBII6Aqdmd0OPAok8svdIg6Y2f3OueeCGU4im2bPiYiUbEjHBgxomcyz/1nJsK6NSKoa73UkiWBlvkJnZlcCz1C8I8RFQCegM8X7p84HnjGzK0IRUiLPodlzAzV7TkTkqMyMB4d1YM+BfJ77j4ZCSPkEusr1e2CQc26sc26pc26Jc24scBLwA3B7KEJK5NHsORGRY2vXoCaX9WvKu7M2sGLrXq/jSAQLpNC1BT50zhUeecJ/7EP/Z0R+nj03WLPnRERKddvpbaieEMsjny3RPq9y3AIpdHuBRqWcT/V/RqLcodlzwzR7TkTkmGpXi+cPg9vww+psvly6zes4EqECKXRfAr83s9OPPGFmpwE3A5ODFUwil2bPiYgE5pI+TWhbvwaPTVxKbv7/3AgTOaZACt3dFO8KMdnMFprZB/7XAorL3k7gnlCElMgydu4mWqRUo7tmz4mIlElsjI8Hh3Vg444DvPH9Wq/jSAQqc6Fzzm0CugHPAfHAuf5XAvAs0N3/GYli/509l6bZcyIiATihVV2GdKzPX79ZzdbduV7HkQgT0GBh51y2c+4O51x751wV/6u9c+5O51x2qEJK5Dg0e25Ed91uFREJ1H1nd6CgyPHnL5Z7HUUijPYakaDR7DkRkfJpklyV605szvj5GczbsNPrOBJBStwpwswePI7v55xzj5Yjj0SwQ7Pn7j67vddRREQi1m9ObsXYuZt4+NMlfPybE/D59PiKHFtpW3+NPo7v5yjeGkyikGbPiYiUX7WEWO4+qx23jVnIuHmbuLBXY68jSQQo8Zarc853HC8NHYtSmj0nIhI853ZNpXuTJJ76YgV7c/O9jiMRQM/QSVBo9pyISPD4fMboYR3J2neQl75Z7XUciQAqdBIUmj0nIhJcXRsncUHPNN78fi1rs/Z7HUfCnAqdlJtmz4mIhMadZ7YlPsbH4xOXeR1FwpwKnZSbZs+JiIRGvRqJ/O601kxZto1pKzO9jiNhTIVOykWz50REQuuqE5rRLLkqj3y2lPzCIq/jSJhSoZNyOTR7ToshRERCIyE2hvuHdmD19n28M2O913EkTB1XoTOzRDNLNbP4YAeSyKLZcyIioXda+3oMapPCc1NWkr3voNdxJAwFVOjMbKCZfQfsBTYAA/3H65rZV2Y2OAQZJUxp9pyISMUwMx48pz0H8gp59j8rvY4jYajMhc7MBgJfAQ2AvwM/L2d0zmX5318T7IASvjR7TkSk4rSqV4Mr+jfj/dkbWLJ5t9dxJMwEcoXuMWAp0Al44CjnpwK9gxFKIoNmz4mIVKxbTm9N7arxPDxhKc45r+NIGAmk0PUC3nLOHaR4z9YjbaL46p1EAc2eExGpeLWqxHHH4LbMXruDzxdt9TqOhJFACl0RRy9yhzQCcsoXRyKFZs+JiHjjot6Nad+wJk98vowDeYVex5EwEUih+xEYfrQT/tWulwHTgxFKwptmz4mIeCfGZ4we1oGMXQd4fdoar+NImAik0D0BnGxm/6T49itAYzM7B5gGNPd/Rio5zZ4TEfFW3xbJDO3SkFemriZj1wGv40gYKHOhc859BVwCnA187j/8JvAp0Aa4xDk3M+gJJexo9pyIiPfuOasdzsGTk5Z7HUXCQEBz6JxzHwJNgBHAXcC9wIVAE+fc2ODHk3Cj2XMiIuEhrXZVbjypJRMWbmb22h1exxGPBbxThHMuxzn3iXPuaefcU865cc65faEIJ+FHs+dERMLHjSe1pFGtRB6esITCIo0xiWbay1UCotlzIiLho0p8DPec3Z4lm/fw0ZyNXscRD5VY6MysyMwKA3wVVGR4qViaPSciEn7O6dKQPs3q8PTkFew+kO91HPFIbCnnHqH0uXMSZTR7TkQk/JgZDw7rwLCXvucvX63i/nM6eB1JPFBioXPOja7AHBLmNHtORCR8dUqtxajejXlr+jpG9WlCq3rVvY4kFUzP0EmZaPaciEh4u2NwW6rEx/DYxKVeRxEPlHiFzswGATjnph3+/lgOfV4qF82eExEJb8nVE7jltNY8NnEZXy/fxqnt9N/raFLaM3TfAs7Mqjjn8g69L+Xz5j+v4WSVzKHZcyN6pGn2nIhIGLuifzPem72BRz9bxsBWKcTH6kZctCit0J0C4C9zAKeiRRJRSbPnREQiQ3ysjwfP6cCv//Ejb01fy/WDWnodSSpIaYWuK/DFoTfOuW9DnkbCkmbPiYhEjpPb1uPUdvV48avVnN89jZQaCV5HkgpQ2rXY54Beh97458xdEvpIEk40e05EJPLcP7Q9BwsKeXqy9nmNFqUVul1A3cPe60/zKKTZcyIikadFSnWuOqE5H83dxE+bdnkdRypAabdcZwL3m1lTYLf/2Agza1XK1zjn3KNBSyeeKipyjJ+XodlzIiIR6HentmL8vE08PGEpY2/sr7sslVxphe63wFvALRRfyXPACP+rJA5QoaskZq7JJmPXAe46q53XUUREJEA1EuO4c0g77hz3E58u3My53VK9jiQhVOItV+fcOufcyUAikEbxLdffAY1LeTUJcV6pQJo9JyIS2S7omUbn1Fr86fPl5ORpu/XK7JgDapxzBc65zcDDwFTnXEZpr9BHloqwNzefzxdvYVjXRpo9JyISoXw+Y/TwDmzdk8sr36Z7HUdCqMwTB51zDzvnFgOYWYqZ9TazXmaWErp44pVJi7Zq9pyISCXQs2kdzuvWiNemrWHjjhyv40iIBDRC2sz6m9lMYCvFiyZmAVvNbLqZ9QtFQPGGZs+JiFQed53Vjhgznvh8mddRJETKXOj8he1roB3wCsXP0/3e/88dgG/MrG8oQkrFWpe1n9nrdmj2nIhIJdGwVhV+e0pLJi3eyvT0LK/jSAgEcoXuMWAL0M45d7Nz7mXn3F+dczdTXPK2+D8jEW78vE2YwfndtSJKRKSyuPbEFqTVrsIjE5ZSUFjkdRwJskAKXV/gNefc1iNP+I+97v+MRLCiIse4eRkMbFWXhrWqeB1HRESCJDEuhvvObs/yrXt5/8eNXseRIAuk0Dn/qySq+5XAodlzWgwhIlL5nNmpAf1bJPPslyvYlZPndRwJokAK3Y/ADWZW98gT/mM3ALODFUy8MXbuJmokxDKkYwOvo4iISJCZGQ8O68CeA/k8P2WV13EkiErbKeJIDwJfASvM7J/ACv/xdsDlQFX//0qEOjR77vzuaZo9JyJSSbVvWJNL+zblnZnrubhPE9o2qOF1JAmCQObQ/QAMBtZSvB3Yy/7X74F0YLBzbnooQkrF0Ow5EZHocPsZbaieEMsjny3BudKeppJIEdAcOufcNOdcL6Ah0N//auic6+Oc+y4UAaXifL18O6lJVejRRLPnREQqs9rV4rn9jDb8sDqbL5du8zqOBEGZCp2ZVTWzuWZ2I4Bzbptzbpb/pV8JlUBhkWPGmmwGtEzW7DkRkShwad8mtKlfnccnLiM3v9DrOFJOZSp0zrkcoAVayVppLduyh90H8jmh1f+seRERkUooNsbHQ8M6smFHDm98v9brOFJOgdxy/QY4OUQ5xGOHJof3b5nscRIREakoJ7Sqy+AO9fnrN6vZtifX6zhSDoEUut8DXc3sOTNrY2aBrJCVMDc9PZuWKdWoXzPR6ygiIlKB7h/agYJCx1NfLPc6ipRDIIVuLdCG4mK3DDhoZnlHvA6GJKWEVH5hEbPX7mBAS91uFRGJNk2Sq3Ltic0ZPy+D+Rt2eh1HjlMgV9n+Rek7RQSNmV0IjAbaA32cc3P8x88AngTigTzgj865r/3nLgbu9WfcDFzmnNMOxGXw06Zd5OQVMkC3W0VEotJvTmnF2LmbGD1hKR/fNACfT4vjIk2ZC51z7tchzHGkxcAI4LUjjmcBw5xzm82sEzAZSPXf/n0B6OCcyzKzPwM3U1wK5Rh+WJ0NQL8WKnQiItGoekIsd5/Vjts/XMj4+RmaRxqBAppDV1Gcc8uccyuOcny+c26z/+0SINHMEgDzv6pZ8cyNmhRfpZMymJ6eRYeGNaldLd7rKCIi4pHzuqXSrXEST32xnH0HC7yOIwEKqNCZWZKZPW5mC81st/+10H+sdqhClmAkMN85d9A5lw/cBCyiuMh1AN6o4DwRKTe/kHnrd3FCK12dExGJZj6fMXp4RzL3HuSlr1d7HUcCVOZCZ2YtgZ+Ae4AYYArFe7vG+I/95P9MWb/fFDNbfJTXuWX42o7AU8AN/vdxFBe67kCjw3KW9PXXm9kcM5uTmZlZ1siV0tz1O8krLNKCCBERoVvjJEb2SOPN79eyLmu/13EkAIFcoXsRqE3xnq2dnHMjnXMjnHOdgCFAkv8zZeKcO93/fY58fVLa15lZGvAxcIVzLt1/uJv/e6a74k3pPgQGlPJjv+6c6+Wc65WSklLWyJXS9PQsYnxG7+Z1vI4iIiJh4K4z2xIXYzw2cZnXUSQAgRS6k4DnnXNTjjzhnPsPxWXupGAFOxozSwImAvc453447FQG0MHMDrWzMygerSLHMD09m65ptaieoLGCIiIC9WomcvOprZmybBvTVkb3XaxIEkih2weU9m92m/8z5WZm55vZJqA/MNHMJvtP3Qy0Ah4wswX+Vz3/QomHgWlm9hPFV+yeCEaWymxvbj4/bdqt260iIvILVw9sRtPkqjzy2VLyC7XrZyQIpNC9D1xmZv+zFNLMEoHLgfeCEco597FzLs05l+Ccq++cG+I//phzrppzrtthr+3+c68659o757o454Y557KDkaUym712B4VFTvPnRETkFxJiY7h/aAdWb9/HuzPXex1HyiCQ+2wTKN7LdYGZvQqspHiIbzvgeuAgMMHMfvHsmnNuenCiSrBNT88mPtZHj6YVvUBZRETC3ent63Fi67o895+VDO/aiOTqCV5HklIEUugOf3buef67a4SV8Bnzfybm+KJJqE1Pz6ZX09okxulfkYiI/JKZ8eA5HTjzhe/4v/+s5PHzO3sdSUoRSKG7KmQppMLt2J/Hsi17uGNwG6+jiIhImGpdvwZX9G/K29PXcWnfpnRoVNPrSFKCQLb+ejuUQaRizVxT/Ihhfy2IEBGRUtx6Whv+PT+Dhycs4YPr+1G8IZOEm7Dc+ktCb3p6FtXiY+iSVsvrKCIiEsZqVY3jD4PbMmvtDiYt3up1HCmBCl2Umr46mz7N6xAXo18CIiJSuov7NKFdgxo8PnEZufmFXseRo9Cf5lFoy+4DrMnazwmtdLtVRESOLca/z2vGrgO8Pm2N13HkKFTootCM9EPPz2n+nIiIlE2/FskM7dyQl79dzeZdB7yOI0dQoYtC09OzSaoaR/sGWq0kIiJld8/Z7XAOnpy03OsocgQVuijjnGNGejb9WyTj82mlkoiIlF1a7arccFJLPl24mR/X7fA6jhymzIXOzDqa2Ygjjp1iZl+Z2Vwz+0Pw40mwbdiRQ8auA9ruS0REjsuNJ7WgYa1ERn+6hMIid+wvkAoRyBW6PwPXHnpjZqnAp0AXoArwZzO7IrjxJNh+WK35cyIicvyqxsdyz9ntWbJ5Dx/N2eh1HPELpND1AKYe9v5Sirf16uac6wBMAn4bxGwSAtPTs6hfM4GWKdW8jiIiIhFqWJeG9G5Wm6cnr2BPbr7XcYTACl1tYNth788CpjrnMvzvJwDaRyqMHXp+bkDLupr0LSIix83MeGhYR3bk5PHilFVexxECK3RZQBqAmVUD+gNTDjsfS2B7w0oFW7ltH9n78zSuREREyq1Tai0u6tWYt6avIz1zn9dxol4ghe474Cb/wojngTiKn6E7pC2QcbQvlPAwPT0LQAsiREQkKO4Y0pYqcTE8+tlSr6NEvUAK3b3AAWAscA3wZ+fcKgAziwEu4JfP2EmYmZ6eTZM6VUmrXdXrKCIiUgnUrZ7ALae35tsVmXyzfLvXcaJamQudc24t0A7oBjR3zt1z2OmqwE3An4IbT4KloLCImWuydXVORESC6or+zWiRUo1HP1tKXkGR13GiVkCDhZ1zBc65n5xz6484vtc594lzbl1Q00nQLNm8h725BXp+TkREgio+1scD53RgTdZ+3p6+zus4USugQmdmdczsUTP7wcxWmVl///FkM3vQzNqFJqaU13T//q0DNH9ORESC7JS29TilbQovfrWKzL0HvY4TlQLZKaIxsAC4E6gBtKB4oDDOuWzgYjSHLmxNT8+iTf3qpNRI8DqKiIhUQg+c04ED+YU8M3mF11GiUqA7RSRS/AzdqcCRg8w+8R+XMJNXUMSP63bo6pyIiIRMi5TqXHVCMz6cu5FFm3Z7HSfqBFLozgBedM4tA462edtaoHFQUklQLdi4i9z8Ij0/JyIiIfW701qTXC2ehycswTnt81qRAil01YDS1iRXL2cWCZHp6VmYQb/mKnQiIhI6NRPj+OOQtsxZv5NPF272Ok5UCaTQrQD6lXL+bGBx+eJIKExfnU2nRrWoVTXO6ygiIlLJXdizMZ1Ta/Gnz5eTk1fgdZyoEUihew24zMyuBmL8x5yZ1TCz54CTgZeDnE/KKSevgPkbdzKgla7OiYhI6Pl8xkPDOrB1Ty6vfpvudZyoUea9V51zr5hZR+DvQI7/8FigFsXF8EXn3LvBjyjlMWfdTvILnRZEiIhIhenVrA7DuzbitWlruKxfU+rVTPQ6UqUX6GDhm4ETKC51k4DZwCvAic65W4MfT8preno2sT6jd7PaXkcREZEo8ofBbSgocrysq3QVosxX6A5xzs0AZoQgi4TAjPQsujdJomp8wP+qRUREjlvT5Gpc2DON92Zt4PpBLWiUVMXrSJVaQFfoJLLsPpDPoozd9NftVhER8cDNp7bC4Xjpm9VeR6n0At3661ozm25mW83soJnlHfHSfh9hZNaabIocDND8ORER8UBa7aqM6t2ED3/cyMYdOcf+AjluZb4PZ2bPALcBmym+5borVKEkOKanZ5MY56N7kySvo4iISJT67SmtGDNnIy9+tYqnL+zqdZxKK5AHq66meCHEuc65whDlkSCakZ5N72Z1SIiNOfaHRUREQqBBrUQu69uUt2es4zentKJ53WpeR6qUArnlasAElbnIkLn3ICu27dV2XyIi4rmbTm5JfIyPF6as9DpKpRVIoZsM9A5VEAmumWuyATR/TkREPJdSI4ErBjTlk4WbWbVtr9dxKqVACt3vgd5m9pCZpYYqkATH9PRsaiTE0qlRTa+jiIiIcMOgllSNi+H5Kau8jlIplbnQOee2A+8ADwIbzCxfq1zD1/T0LPq2qENsjCbTiIiI9+pUi+fqgc2ZuGgLSzfv8TpOpRPIKtdHgPsoXuU6B61yDVubduawPjuHK/s38zqKiIjIz64d2IK3pq/juSkr+dsVvbyOU6kEssr1BrTKNSLMSPc/P9dKCyJE5P/bu/M4Oeo6/+Ovz1yZ3JM7hNyZDJBwBAmQgyMQDkUR1GUBF8FrUQGXQ3dXVt2fuLrr6nKpILC7ioAueC6LiJAQCCGTAAESIIYcPeSGkOncmVwz8/n9UTXYdHrunqmqmffz8ehHd3/rW9/6fKt7Zj7zrW9VicRH/17F/O3p47ltzipe27iD40fqslr50prjcaXoLNdEWJRKM6h3CRVD+0YdioiIyPt8ZuZYynoVc9sc5Pk1MQAAIABJREFUnfGaT61J6J5CZ7nGnrtTmUozbcIgCgos6nBERETep29pMV84YwLPrtzKy+u2Rx1Ol9GahO464CQzu0VnucbXW9V7eWfXft3uS0REYuuqGWMY3KeE2+asjDqULqM1Cd1G4FjgG+gs19hamNL150REJN56lRTxxTMnsHBN+r3rpkr7tOakiF8A3lGBSH4sSlVzRP9Sxg7qFXUoIiIijbpi2hj+c0EVtz21ike+MA0zTRNqjxYndO7+6Q6MQ/Kgvt5ZlEpz9tHD9IMhIiKxVlpcyLVnlfPPjy7n+TXVnD5xSNQhJZquOtuFvPnObrbXHNL8ORERSYRLTx7FiP6l3PrUKtx1ELA9Gh2hM7MzANz9ucz3zWmoL52vMlUNwHQldCIikgA9igr58uyJ3Py713lm5bucffSwqENKrKYOuT4LuJn1dPeDDe+bqG/h8sK8RSetsiiVZtzg3owo6xl1KCIiIi3yVyeN5CfPprhtzirOOmqopgy1UVMJ3VkAYTL33nuJp9q6el54axsfnTIi6lBERERarLiwgL+bPZGv/noZTy7fwgePHR51SInUaELn7vObei/x8tqmnew5UKv5cyIikjgXTxnB3c+s4fY5qzhv0jBdGL8NWnxShJnNM7PZTSw/y8zm5Scsaa2G+7dOH6+ETkREkqWosIDrz5nIyi27efz1t6MOJ5Fac5brLKCp2YpDgTPbFY20WWWqmqOH92VQnx5RhyIiItJqFx4/gophfbhj7irq6nXGa2u19rIlTe3hCcCedsQibbT/UB1L1m7X3SFERCSxCgqMG8+pILV1L48u3RR1OInT5IWFzexTwKcyim42s8/kqFoGnAjMyWNs0kKvrt/Bgdp6zZ8TEZFEO3/ycCYd0Y87n17NhSeMoLhQl8ttqeb21EBgYvhwYHjG+4ZHOdCT4NZgV3dYpNKoRalqCgxOGT8w6lBERETarKDAuOncCtala/jdKxujDidRmhyhc/c7gTsBzKweuMHdf9kZgUnLLUylOW5kGf1Ki6MORUREpF1mHzOUE0aV8cOn1/CxE0dSUqRRupZo8V5y9wIlc/Gz50AtyzbsYKYOt4qISBdgFozSbdqxj0eWbIg6nMRQ2ptwL63dRm2964QIERHpMs6YOJipYwZw17w17D9UF3U4iaCELuEWpdKUFBZw0pgBUYciIiKSF2bGTedV8M6u/fzyhfVRh5MISugSrjJVzYmjy+hZolvoiohI1zFjwmCmjx/E3c+m2HdQo3TNUUKXYDtqDrJ88y4dbhURkS7pK+dVUL3nAA8sWht1KLGnhC7BFlelcYcZ5TohQkREup6pYwdyRsUQ7pmfYs+B2qjDiTUldAlWmUrTq6SQE0aWRR2KiIhIh7jp3Aq21xzi/oVvRR1KrLUooTOzMjM7xcwmNFFnnJldmb/QpDmVqTQnjx2oa/SIiEiXNWVUGeccM5T7nqti575DUYcTW81mAmb2XWALsAhYZWYvmtkHclSdAfwsz/FJI97dtZ817+7R7b5ERKTLu/HcCnbtr+W/n9coXWOaTOjM7K+Bm4H5wLXAvwKjgUozu6Ljw5PGLKpKA+iECBER6fImj+jPh44dzk+ff4vtew9GHU4sNTdCdwPwjLuf5+73uPs3gWOA54D7zeyGDo9Qcqpck6ZfaRGTRvSLOhQREZEOd+O5Few9WMt9C6qiDiWWmkvojgZ+m1ng7tuBDwIPALea2b90UGzShMqqaqaNH0RhgUUdioiISIerGNaXC48fwf0L11K950DU4cROcwldfa5Cd693988CdwBfN7Mft6AtyZMN22rYsG2f5s+JiEi3cv05EzlQW8c9z6aiDiV2mkvCVgGnNbbQ3b8C3AJcA2ikrpNUpqoBmFmu+XMiItJ9TBjSh4+dOJIHF69jy679UYcTK80ldE8AF5lZo0NB7n4LcCMwKl9BmdklZrbczOrNbGpG+blm9rKZvR4+n52x7FIzey1c7/v5iiWOKlNpBvfpQfnQPlGHIiIi0qmunz2Runrn7mfWRB1KrDSX0P0U+AdgaFOV3P1O4BPAt/MU1xvAxwlOvshUDVzo7scBVwEPAoQJ5w+A2e4+GRhmZrPzFEusuDuVqTQzJgzCTPPnRESkexk9qBeXTB3J/7y4gU079kUdTmw0mdC5+yZ3v8vdVzTXkLv/bzha127uvsLdV+Yof9XdN4dvlwOlZtYDGA+scvet4bK5BAlml5Pauoetuw9o/pyIiHRb1509EYAfz9MoXYM2n8hgZqVmdqWZDctnQK3wCeBVdz8ArAGONrOxZlYEXEweDwHHSWVK158TEZHu7ciynlx2yih+vWQD69M1UYcTC+05M7U/wZ0hJrdlZTOba2Zv5Hhc1IJ1JwP/DnwB3ruUypeAR4AFwFqg0bv4mtnVZrbEzJZs3bq1sWqxtHBNNUeW9WTUwJ5RhyIiIhKZa88qp7DA+OG81VGHEgvtvdRImydxufs57n5sjsejTW7QbCTwe+BKd3/vvGV3f8zdT3X36cBKoNFP2N3vc/ep7j51yJAhbe1Cp6urdxZXbWNmuebPiYhI9zasXylXTBvD717ZSNXWPVGHE7n2JnSelyhayMzKgMeBm919YdayoeHzAILLqPxXZ8bWGVa8vYud+w7pcKuIiAjwpVkT6FFUyJ1Pa5QushG6Jhs1+5iZbQSmA4+b2ZPhouuAcuCbZrY0fDScgXunmf0ZWAh8z91XdURsUWq4/tx0nRAhIiLC4D49uGrGWP5v2WZWbdkddTiRak9CtxUYR5BA5ZW7/97dR7p7D3cf5u7nh+Xfcffe7j4l4/FuuOxyd58UPh7Od0xxUJlKM2FIb4b1K406FBERkVj4whnj6V1SxB1zu9w4Tqu0OaELb/+1LjzLVDrYobp6Xnxrmw63ioiIZBjQu4TPzhzLH19/h+Wbd0YdTmR0/9WEWLZhBzUH63T9ORERkSyfO308/UqLuH1O951Lp4QuISpTacxg2ngldCIiIpn69yzmb08fz9wVW1i2YUfU4URCCV1CVKaqmXREPwb0Lok6FBERkdj5zGnjGNCrmNvmdM+5dEroEmD/oTpeWbdDh1tFREQa0adHEV84cwLzV23l5XXbog6n0ymhS4CX123nYF29TogQERFpwpXTxzC4Twm3PtX9RumU0CVAZaqawgLj5HEDow5FREQktnqVFPGlWeVUptIsCu993l0ooUuAylSaE0b2p0+PoqhDERERibW/OXU0w/r14LY5K3Hv1BtaRUoJXczt3n+I1zbuZGa5DreKiIg0p7S4kOvOKueltdtZsLo66nA6jRK6mHvxrW3U1btu9yUiItJCf33yKI4s68mtc1Z1m1E6JXQxV5lKU1JUwAdGD4g6FBERkUToUVTIl88uZ9mGHcx7892ow+kUSuhirjKVZuqYAZQWF0YdioiISGJ84qSRjB7Yi9u6ySidEroY27b3ICve3qXrz4mIiLRScWEB18+eyPLNu3hy+TtRh9PhlNDF2OKq4JTr6br+nIiISKtdfOKRjB/Sm9vnrKa+vmuP0imhi7GFa6rpXVLI8SP7Rx2KiIhI4hQWGDecU8HKLbv5w+tvRx1Oh1JCF2OLUmlOHT+I4kJ9TCIiIm3xkeOO4Khhfblj7ipq6+qjDqfDKFOIqbd37qOqeq/mz4mIiLRDQYFx47kTqdq6l0eXbo46nA6jhC6mGm5ZouvPiYiItM/5k4czeUQ/7nx6NYe66CidErqYqkylKetVzDHD+0UdioiISKKZGTedW8H6bTX89uWNUYfTIZTQxZC7syiVZvr4QRQUWNThiIiIJN7ZRw9lyqgyfjRvDQdq66IOJ++U0MXQunQNm3bs0/w5ERGRPGkYpdu0Yx+/emlD1OHknRK6GKoM58/NKNf150RERPLl9ImDOXnsAH78zBr2H+pao3RK6GKoMlXNsH49GD+4d9ShiIiIdBnBKN1RbNl1gF+8sD7qcPJKCV3MNMyfmzFhMGaaPyciIpJP0ycMYsaEQfzk2TXUHKyNOpy8UUIXM6u27CG996AuVyIiItJBvnJeBdV7DvLAonVRh5I3SuhipjJVDaATIkRERDrISWMGcmbFEO6dn2LPga4xSqeELmYqU2lGD+zFyAG9og5FRESky7rp3Aq21xziZ8+/FXUoeaGELkZq6+pZXJVmZrlG50RERDrSCaPKOOeYYfzngip27jsUdTjtpoQuRpZv3sXu/bVMn6DLlYiIiHS0m86tYNf+Wv57QVXUobSbEroYabj+3PTxGqETERHpaJNG9OOC44bz04Vr2b73YNThtIsSuhipTFVTMawPQ/r2iDoUERGRbuGGcyrYe7CWe59L9iidErqYOFhbz0trtzFDh1tFREQ6TcWwvnz0hBH8vHItW3cfiDqcNlNCFxNLN+xg/6F6XX9ORESkk10/eyIHauu4Z34q6lDaTAldTCxcU02BwTTNnxMREelU44f04eMfGMlDi9exZdf+qMNpEyV0MbEolebYI/vTv2dx1KGIiIh0O9fPnkhdvXPXM2uiDqVNlNDFQM3BWl7dsF2HW0VERCIyamAvLpk6iodf3MCmHfuiDqfVlNDFwJK12zlU5zohQkREJEJfPrscgB/PWx1xJK2nhC4GKlNpigqMk8cOiDoUERGRbmtEWU8uP2UUv16ykfXpmqjDaRUldDGwKFXNiaPL6FVSFHUoIiIi3dq1Z5VTWGDc+XSyRumU0EVs575DvL5ppw63ioiIxMDQfqV8atoYfv/qRlJb90QdTospoYvYC1Vp6h1m6IQIERGRWPjirAmUFhdy59zkjNIpoYtYZSpNaXEBU0aXRR2KiIiIAIP79OCqGWN57LXNrHxnd9ThtIgSuogtSqU5eexAehQVRh2KiIiIhK4+fTy9S4q4Y+6qqENpESV0Edq6+wArt+zW9edERERiZkDvEj572jieeOMdlm/eGXU4zVJCF6HFVWkAnRAhIiISQ587bRz9Sou4fU7859IpoYtQZSpN3x5FHDuiX9ShiIiISJb+PYu5+ozxzF2xhWUbdkQdTpOU0EWoMlXNqeMHUVSoj0FERCSOPj1zHAN6FXPbnHjPpVMmEZGN22tYl67R5UpERERirE+PIr545gTmr9rKy+u2RR1Oo5TQRWRRKpw/V66ETkREJM6unD6WwX16cOtT8R2lU0IXkUWpNIN6l1AxtG/UoYiIiEgTepYUcs2sCVSm0u8NyMSNEroIuDuVqTTTJgyioMCiDkdERESa8clTRzO8Xym3zVmJu0cdzmGU0EXgreq9vLNrv+bPiYiIJERpcSHXnl3OS2u3s2B1ddThHEYJXQQWhsO1M3X9ORERkcS4dOoojizrya1zVsVulE4JXQQWpaoZ0b+UMYN6RR2KiIiItFBJUQF/N7ucZRt2MO/Nd6MO532U0HWy+npnUSrN9AmDMdP8ORERkST5+AdGMmZQL26L2SidErpO9uY7u9lec0jz50RERBKouLCA62dPZPnmXTy5/J2ow3mPErpOVpkKJlJOV0InIiKSSBdNOZIJQ3pz+5zV1NfHY5ROCV0nW5RKM25wb0aU9Yw6FBEREWmDwgLjhnMqWLllN394/e2owwGU0HWq2rp6Xnhrmw63ioiIJNyHjzuCo4f35Y65q6itq486HCV0nem1TTvZc6CWGbpciYiISKIVhKN0VVv38ujSzVGHo4SuMzXcLmTa+IERRyIiIiLtdf7kYRx7ZD/ufHo1hyIepVNC14kqU9UcPbwvg/r0iDoUERERaScz46ZzK1i/rYbfvrwx0liU0HWS/YfqWLJ2uw63ioiIdCFnHTWUKaPK+NG8NRyorYssDiV0neTV9Ts4UFuvEyJERES6EDPjK+dVMKhPCVt3H4gsjqLIttzNVKaqKSwwTtX8ORERkS7ltPLBnFYe7R2glNB1kspUmuOO7E/f0uKoQxEREZE8isOtPHXItRPsOVDLsg07dLhVREREOoQSuk7w0tpt1Na7TogQERGRDqGErhMsSqUpKSzgpDEDog5FREREuqDYJnRmdomZLTezejObmlF+ipktDR/LzOxjGcs+aGYrzWyNmX0tmsgPV5mq5sTRZfQsKYw6FBEREemCYpvQAW8AHweey1E+1d2nAB8E7jWzIjMrBO4CPgRMAi43s0mdGXAuO2oOsnzzLh1uFRERkQ4T27Nc3X0FHH7miLvXZLwtBTx8fQqwxt2rwvUeBi4C/tzhwTZhcVUad5hZrhMiREREpGPEeYSuUWZ2qpktB14HvujutcCRwIaMahvDskhVptL0Kink+JFlUYciIiIiXVSkI3RmNhcYnmPR19390cbWc/cXgMlmdgzwczN7Ash1ERjPUYaZXQ1cDTB69OhWx90aPYsLOX/ycEqKEpk7i4iISAJEmtC5+zntXH+Fme0FjiUYkRuVsXgksLmR9e4D7gOYOnVqzqQvX26+4JiObF5EREQkeYdczWycmRWFr8cARwFrgZeAieHyEuAy4P8iC1RERESkk8Q2oTOzj5nZRmA68LiZPRkuOg1YZmZLgd8D17h7dTiP7jrgSWAF8Ct3Xx5F7CIiIiKdydw79Ihj7E2dOtWXLFkSdRgiIiIizTKzl919anZ5bEfoRERERKRllNCJiIiIJJwSOhEREZGEU0InIiIiknBK6EREREQSTgmdiIiISMIpoRMRERFJOCV0IiIiIgmnhE5EREQk4ZTQiYiIiCScEjoRERGRhFNCJyIiIpJwSuhEREREEk4JnYiIiEjCKaETERERSThz96hjiJSZbQXW5aGpwUB1HtpJKvVf/Vf/uy/1X/1X/zvPGHcfkl3Y7RO6fDGzJe4+Neo4oqL+q//qv/ofdRxRUf/V/zj0X4dcRURERBJOCZ2IiIhIwimhy5/7og4gYup/96b+d2/qf/em/seA5tCJiIiIJJxG6EREREQSrtsldGb2QTNbaWZrzOxrOZb3MLNHwuUvmNnYjGU3h+Urzez85to0s3FhG6vDNkua20a4fLSZ7TGzr3an/pvZWDPbZ2ZLw8c93an/4bLjzWyRmS03s9fNrLS79N/M/ibjs19qZvVmNqUb9b/YzH4efu4rzOzmfPY9Af0vMbOfhf1fZmazumj/zzCzV8ys1sz+Kmv7V4X1V5vZVd2w/38ysx1m9od89z3u/TezKfaX3/2vmdmlre6gu3ebB1AIpIDxQAmwDJiUVeca4J7w9WXAI+HrSWH9HsC4sJ3CptoEfgVcFr6+B/hSU9vIiOG3wK+Br3an/gNjgTe66+cPFAGvASeE7wcBhd2l/1lxHAdUdbPP/5PAw+HrXsBaYGw36v+1wM/C10OBl4GCLtj/scDxwAPAX2VseyBQFT4PCF8P6C79D5fNBi4E/pDPn/0k9B+oACaGr0cAbwNlrepjvndanB/AdODJjPc3Azdn1XkSmB6+LiK4WKBl122o11ib4TrVQFH2thvbRvj+YuAHwLfIf0IX6/7T8Qld3Pt/AfBQd+1/Vhz/Cny3O/UfuBx4LCwbBKwCBnaj/t8FXJHR1tPAKV2t/xl17+f9f9AvB+7NeH8vcHl36X9G+Sw6JqFLRP8zli8jTPBa+uhuh1yPBDZkvN8YluWs4+61wE6CX66NrdtY+SBgR9hG9rZybsPMegP/CNzS5h42Ldb9D5eNM7NXzWy+mZ3etm42Ku79rwDczJ4Mh+T/oc09zS3u/c90KfA/repd8+Le/98Aewn+M18P/Ie7b2tbV3OKe/+XAReZWZGZjQNOAka1sa+5xKX/7YmvPeLe/46WmP6b2SkEI36plq4DQQbanViOMm9hncbKcyXFTdVvahu3ALe7+x6zXFXaLe79fxsY7e5pMzsJ+F8zm+zuu3LUb4u4978IOA04GagBnjazl9396Rz12yLu/Q8Wmp0K1Lj7GznqtUfc+38KUEdwuGUAsMDM5rp7VY76bRH3/v8UOAZYQnA7xkqgNkfdtopL/xvTlnVaI+7972iJ6L+ZHQE8CFzl7vUtWadBdxuh28j7/+MbCWxurI6ZFQH9gW1NrNtYeTVQFraRva3GtnEq8H0zWwvcAPyTmV3Xtq7mFOv+u/sBd08DuPvLBP+dVLSxr7nEuv9h+Xx3r3b3GuCPwAfa2Ndc4t7/BpeR/9G59203R0yH1Ymg/58E/uTuh9z9XWAhkM/bCcW6/+5e6+43uvsUd78IKANWt7m3h4tL/9sTX3vEvf8dLfb9N7N+wOPAN9x9cYt6lSnfx6nj/CAYAakimNTYMIFxclada3n/pMhfha8n8/5JkVUEEyIbbZPgxIbMSZHXNLWNrDi+Rf7n0MW6/8AQwpMACCaZbiK/c4ji3v8BwCsEE+KLgLnAh7tL/8P3BQS/JMd3t59/gukWPyP477438Gfg+G7U/15A7/D1ucBzXfHzz9jW/Rx+UsRbBL8HBoSvu9zvv8b6n1E+i46ZQxfr/ofrPw3c0OY+5nunxf1BMPF8FcHoz9fDsm8DHw1fl4YfxBrgRTL+sABfD9dbCXyoqTbD8vFhG2vCNns0t42Mdb9FnhO6uPcf+ASwPPyheAW4sDv1P1x2RbgP3gC+3w37PwtY3B1//oE+YflygmTu77tZ/8eGba8g+GdmTBft/8kE/7TsBdLA8ox1PhvWXwN8phv2fwGwFdgX1jm/u/Sf4Hf/IWBpxmNKa/qnO0WIiIiIJFx3m0MnIiIi0uUooRMRERFJOCV0IiIiIgmnhE5EREQk4ZTQiYiIiCScEjoRaZSZfd7M3MxGRh1LZzKzCWb2RzPbHvY/nxf47hBm9ryZzY06jtYIb/PlZvaNqGMRSToldCIJYGa/M7NDZjakiTrXh38cL+zM2LqonxPcpeFbwKeApyKNRkSkGUroRJLhQYKrkl/WRJ0rCG4586c8bvdnQE9335jHNmPNzEqBGcAD7n6nuz/k7quijktEpClK6ESS4XGCewpekWuhmVUQjCg97O6H2rsxM+sN4O517r6/ve0lzFCC22/taMvKFuiV35BERJqmhE4kAdz9IPAr4BQzm5ijyqfC5wcbCszsYjN7zMw2mtnB8Pmu8AbQZNT7TniodpKZ3W9m1cDacNlhc+jM7Ewze8TM1pnZATPbYmYPmNmIrHYb1j3dzP7dzN4xs31m9qSZjc3ugJmVm9lDYb0DZvaWmd3bkFyGdfqZ2X+Y2dqwzvqw7dKW7McwlmfMbI+Z7TazOWZ2aua+ANaFb/8ljL+2ifYa5oDdY2aXmNkyYD/BLZwws8+a2dyMPlWZ2XfNrCSrnYfMbL+ZDTez34SxbQ/7X5pV18zsZjPbEO7Pysw+ZNUdFMb2drj9FWZ2k5lZI334qJktC9tdamZnhHU+bGavhjGuMLPzWri/LzGzF81sp5ntNbPVZvbjFqw3xsx+aWbV4TaXmtmVWXXKw7i/ZmbXmFkqrPuKmZ2To812fXdE4q4o6gBEpMUeBL4I/A3B3K5MnwRWufuLGWWfAw4CPwK2AycCnweOBc7M0f6vgPXAPxPcHL4xlwJlwH3Au8BRwNXAqWZ2Qo4RvdsJ7s34r8AQ4KthX05vqGBmxxLcx7EgbHc1MJLg/r4DgL1m1hN4luBG2PcR3DtxCnATwc2zP9JEzJjZWcCTwAbgO+G2vgjMN7NZ7r6Y4J6L1WHMvwEeBeqbajd0OnAJcBdwN8H9SAG+THBf3icI7t04E/ha2LerstooCON7Dfh7YBrBfn0X+GZGvW8D3wDmAP9HsP+fIBhRrMrobynB/joG+AnwJvBh4FZgNHBD1vanAxeG8e8D/gF43Mz+NlznbqAmLP+NmY1y952N7ZAw6XsEmAf8E1BLcH/L8xtbJ1xvKFBJ8B37EbCZYKrBz81soLvfkbXKZQSjqncTfN+/EMY9y90XhW2267sjkgj5vvmvHnro0XEPghs9r84qmwk48I2s8l451v90WPfUjLLvhGW/zVH/8+Gykc20Oyusd2mOdecDBRnlXw3Lj8ooe44gWajI0XbDPae/SZBoTMpafk3Y3lnN7LulBDfDHpJRdiSwG6jMKBuba3820mZRWLeeHDfSbmRffQuoA47IKHsobOc7WXX/ALyd8X4oQdIyJ2ufNuyDuRllN4Rln8vcl8DvwniPyurDwcz9T5DcOcGIY3mO8s83s29+RDBNoLAF++8bGWV3hGWzM8pKgBcIkuKysKw8I+4JGXWHAbuABRll7fru6KFHEh465CqSLA8B5WY2LaPsCoI/Sr/IrOjuNfDeIbr+ZjYYeD5cfFKOtn/SkgAa2g3b7hu2+wZBYpSr3XvdPXOUa374PD5sYxjBCNf9nuPkA3f38OWlwELgXTMb3PAAGi7VcXZjMVtwyPgE4OfuvjWj7U3AL4HpZjaoiW43p9Ldl+aIveEzKDCzsjDe5whG407M0c7dWe/nA8PtL3PyzgeKgR9m7dP/Itj/mT5CMNp4f0Y8DvyAILH7cFb9Z7L2/6KGGNx9TY7y8Tniz7QD6Au06PBsho8Ar7r70xlxHyRI9Hpx+Of8mLunMupuAf4HmGlmZWFxm787IkmhhE4kWR4Kn68ACOdi/TXwvLu/lVnRgjlxjwF7CP64biU4lAnB4axsqRxlhzGzkWb2CzPbSTASsjV89G2k3XVZ77eHzwPD5/Lw+fVmNl0BzM7YXsNjZbh8aBPrjg2f38yx7M9Zddoi574zsxlm9gzB6ON2gngbEpXsfXXI3TdnlTXsqwHh85jweWVmpTDhWZu17liC0dy6rPKG/o7LKl+f9X5HM+UDaNqPwzj/GM4h/KWZXWpmzU31GUPTn1N23CuzK4Zlxl/2V3u+OyKJoDl0Igni7mvMbBFwqZndAFxAkBg9mFkvHJmYT3CI6psEh2prCA5dPU7uf+b2Nbd9MyskGNUYDPw7wR/ZvQQjhL9upN3shOK95rKevZF6mfXnAf/WyPJNzazfVLst2X5TDtt3ZjaBIHlbTTBXaz3B4cvRwH8l9wWVAAAD00lEQVRz+L5qaq5eS/aV5ShrSnYbjX1OzX1+uRt332JmJxIkUucTjNRdDiwxszPcvdnvWyPby467Jfuio747IrGhhE4keR4kODR3PsFI3QGCZCrTbIKk62J3X9hQaGaT2rntKQST8K9w9/cO8VpwJmr/NrbZcDjv+GbqpYC+7t6WuyGsDZ+PzrGsoSx7JLG9LgZKgQ+Fh3YBMLML2tHm2vD5aP4y2towUjsG2JJV93gzK8wapTsmq60O48EldP4UPjCzLwM/BD5O1hSBDOto+nNa20h5pgqCRK/hM23Pd0ckEXTIVSR5HiGYCH4twXyjx9w9+5ppDaM92T/jf9/ObTfVbmtHiABw93cIznD9tAXX03ufjEtsPAycbGYfzVGn1Mz6NrGNjQQnRVwZzp1qWG8EwVnDle6ebkv8TThsX4UjnDe1o82nCM4W/bKZZX4Gnyc45J3pMYKzit+73Ee4LxtOSnm8HXE0q5E5ia+Gz7kOzTd4DDgxPCu5oa1i4HqCUeZ5WfUvDEdDG+oOIxgJrMz4uWjzd0ckKTRCJ5Iw7r7NzP5IMAIEWYdbQwsIzjD8hZn9iOBQ34UEo3btsZxgRO0OC64lV01whusp/GW+V1tcRxDzS2Z2H7AKGEFw2ZILgI3A9wkm8v/OzB4CXiQ4hHwUwTzCi/jLSR+53ERwWZDFZvafBAnolwhOMvhKO2JvzBPA94Anwj4VEUzOL25rg+FhzB8ANwN/MrNHCUdMgbeyqt9HkOjdFx76XEmwLy8A7nT3XHPP8un+8ND/PILDzUMJLhOzhyBpa8y/EXyej4Xf3c0E+20acGOOf17+DCwws7uAQwSXLSkF/jGjTnu/OyKxp4ROJJkeJEjo0gSJw/u4e3V4aO8/gP9HkND9kWC05u22btTdD5rZRwiu09aQBD1LcJbggna0+5oFF8e9BfgM0IfgD/mThImiu9eY2SyCP9SXEozC7CG49toPCZLNprbxTHjB2VsIrrXnwGKCS60sbmvsTWzvTTO7mOCyMN8DdhJc6++nBKOFbfV1gnmLXyI4O/hV4EME14rL3P6+cH99lyBpGUiQ9H0VuK0d22+pBwgusHx1uO1qgjNN/8Xds0+0eI+7v2tmMwj22dUE10RcCVzl7g/kWOVhgpNzvkJwfb8/Ax/JnGrQ3u+OSBI0XN9JREQkMcysnGAe4c3u/r2o4xGJmubQiYiIiCScEjoRERGRhFNCJyIiIpJwmkMnIiIiknAaoRMRERFJOCV0IiIiIgmnhE5EREQk4ZTQiYiIiCScEjoRERGRhFNCJyIiIpJw/x/aiL6gvRbYLgAAAABJRU5ErkJggg==\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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}