{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Mixed Effects Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:19.619191Z",
     "iopub.status.busy": "2021-02-02T07:04:19.618414Z",
     "iopub.status.idle": "2021-02-02T07:04:23.271247Z",
     "shell.execute_reply": "2021-02-02T07:04:23.272026Z"
    }
   },
   "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\n",
    "from statsmodels.tools.sm_exceptions import ConvergenceWarning"
   ]
  },
  {
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:23.281974Z",
     "iopub.status.busy": "2021-02-02T07:04:23.281318Z",
     "iopub.status.idle": "2021-02-02T07:04:24.763633Z",
     "shell.execute_reply": "2021-02-02T07:04:24.762938Z"
    }
   },
   "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.395    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(method=[\"lbfgs\"])\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:24.771056Z",
     "iopub.status.busy": "2021-02-02T07:04:24.770363Z",
     "iopub.status.idle": "2021-02-02T07:04:26.101859Z",
     "shell.execute_reply": "2021-02-02T07:04:26.102885Z"
    }
   },
   "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:               6.0372    \n",
      "Min. group size:   11       Log-Likelihood:      -2217.0475\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.739    0.550 28.603 0.000 14.660 16.817\n",
      "Time              6.939    0.080 86.925 0.000  6.783  7.095\n",
      "Group Var        19.503    1.561                           \n",
      "Group x Time Cov  0.294    0.153                           \n",
      "Time Var          0.416    0.033                           \n",
      "===========================================================\n",
      "\n"
     ]
    }
   ],
   "source": [
    "md = smf.mixedlm(\"Weight ~ Time\", data, groups=data[\"Pig\"], re_formula=\"~Time\")\n",
    "mdf = md.fit(method=[\"lbfgs\"])\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:26.107896Z",
     "iopub.status.busy": "2021-02-02T07:04:26.106481Z",
     "iopub.status.idle": "2021-02-02T07:04:26.129692Z",
     "shell.execute_reply": "2021-02-02T07:04:26.130682Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10324316832591753"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    ".294 / (19.493 * .416)**.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:26.135220Z",
     "iopub.status.busy": "2021-02-02T07:04:26.133799Z",
     "iopub.status.idle": "2021-02-02T07:04:27.465923Z",
     "shell.execute_reply": "2021-02-02T07:04:27.466990Z"
    }
   },
   "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:               6.0283    \n",
      "Min. group size:   11       Log-Likelihood:      -2217.3481\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.739    0.554 28.388 0.000 14.652 16.825\n",
      "Time              6.939    0.080 86.248 0.000  6.781  7.097\n",
      "Group Var        19.837    1.571                           \n",
      "Group x Time Cov  0.000    0.000                           \n",
      "Time Var          0.423    0.033                           \n",
      "===========================================================\n",
      "\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, method=[\"lbfgs\"])\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:27.472095Z",
     "iopub.status.busy": "2021-02-02T07:04:27.470634Z",
     "iopub.status.idle": "2021-02-02T07:04:28.089654Z",
     "shell.execute_reply": "2021-02-02T07:04:28.090612Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:28.097529Z",
     "iopub.status.busy": "2021-02-02T07:04:28.096856Z",
     "iopub.status.idle": "2021-02-02T07:04:29.125134Z",
     "shell.execute_reply": "2021-02-02T07:04:29.124247Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:29.133637Z",
     "iopub.status.busy": "2021-02-02T07:04:29.132852Z",
     "iopub.status.idle": "2021-02-02T07:04:34.069469Z",
     "shell.execute_reply": "2021-02-02T07:04:34.070325Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/travis/build/statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py:2237: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n"
     ]
    },
    {
     "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"
     ]
    }
   ],
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:34.073792Z",
     "iopub.status.busy": "2021-02-02T07:04:34.073119Z",
     "iopub.status.idle": "2021-02-02T07:04:51.235302Z",
     "shell.execute_reply": "2021-02-02T07:04:51.236266Z"
    }
   },
   "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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:51.240032Z",
     "iopub.status.busy": "2021-02-02T07:04:51.239250Z",
     "iopub.status.idle": "2021-02-02T07:04:51.250715Z",
     "shell.execute_reply": "2021-02-02T07:04:51.251569Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:51.254977Z",
     "iopub.status.busy": "2021-02-02T07:04:51.254251Z",
     "iopub.status.idle": "2021-02-02T07:04:51.954909Z",
     "shell.execute_reply": "2021-02-02T07:04:51.955843Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '-2 times profile log likelihood')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn4AAAHoCAYAAADXB+KjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABjQklEQVR4nO3dd3iV5f3H8fc3mxEIkLA3hD1EgggIKrgBFVetilp33a2tra22arWtrb+2bkWqVq3WLQJOHIDgCiNhj4SdMMIMK2Tcvz/OwaYxhJzkJM8Zn9d1nSs5z8rnOTmcfLnv535uc84hIiIiIpEvxusAIiIiIlI/VPiJiIiIRAkVfiIiIiJRQoWfiIiISJRQ4SciIiISJVT4iYiIiESJOK8DhIPU1FTXuXNnr2OIiIiIHNW8efMKnHNpla1T4VcNnTt3JjMz0+sYIiIiIkdlZuuOtE5dvSIiIiJRQoWfiIiISJRQ4SciIiISJVT4iYiIiEQJFX4iIiIiUUKFn4iIiEiUUOEnIiIiEiVU+ImIiIhECRV+IiIiIlFChZ+IiIhIlFDhJyIiIhIlVPiJiIiIRAkVfiIiIiJRQoWfiIiISJRQ4SciIiISJVT4iYiIiEQJFX4iIh5yzlFcWuZ1DBGJEnFeBxARiTa79xczJ6eAWSu3MXtVAZv3HGRg+6aM6J7K8G6pHNsphcS4WK9jikgEUuEnIlLHSkrLWLhhF7NWFTB71TayNuyizEFyYhzDu7dg7IA2fLtmB098vprHPltNYlwMQzo3Z3j3Fozolkq/dk2JjTGvT0NEIoAKPxGROrBhx35mrdrGrJXbmLt6O4VFJcQYDGifws2j0xmVnsoxHVKIi/3vFTd7DhbzTe4O5uYUMHf1dv7y4QpgBclJcRzftQUjurVgePdU0ls2xkyFoIgEToWfiEgQ7C0q4euc7cxa5eu+XVOwD4C2TZMYO6ANI9PTGNG9BSkNE454jCZJ8ZzapxWn9mkFwLbCIubmFPBVznbm5BTwydItAKQlJzK8m681cFi3FnRo3rDuT1BEIoI557zOEPIyMjJcZmam1zFEJISUlTmW5O35vlVv/vqdFJc6GsTHcnzX5ozqkcbI9DS6pTUKWuvchh37mZtTwJzV25mbs52CvUUAdGzekBHdWzDcXwimNk4Mys8TkfBkZvOccxmVrlPhd3Qq/EQEYMueg98PyPhydQE79h0CoE+bJozqkcao9FQGd25WLwMznHOs2rqXOat9heA3ub7uZIBerZMZ3i2V4d1aMLRrc5KT4us8j4iEDhV+taTCTyQ6HSwu5ds1O74v9lZsKQQgtXEio9JTGdkjlRO6p5GW7H0LW0lpGYvz9jBndQFzcwrIXLuTopIyYmOMAe2bft81fGynZiTFa8SwSCRT4VdLKvxEooNzjpVb9jJr5TZmrdrGt2t2UFRSRkJsDEO6NGNkehqj0tPo1TqZmBAfZXuwuJT563cyd7Xv+sDsjbspLXMkxsWQ0bnZ9y2C/ds1/Z8BJiIS/lT41ZIKP5HItWPfIWb7B2TMXrWNLXt81811b9mYUelpjOyRyvFdWtAgIbxbyQoPFvPtmh3+6wMLWL7Z13qZnBjH0K4tfC2C3VPp0UojhkXCXVWFn0b1ikhUOVRSxvz1O5m9ahuzVhawOG83zkHTBvGckJ7q68JNT6NtSgOvowZVclI8Y3q3Ykxv34jhgr1FfJXjGyQyN6eAGct8I4ZTGycwrFsqI/yFoEYMi0QWtfhVg1r8RMKXc4612/f7r9Pbxlc529l3qJTYGOPYjim+7tseafSP8pskb9y5n7n+1sA5OdvZVuhr+ezQvAHDu6Yy3D9qOBSuZxSRqqmrt5ZU+ImElz0Hi5m7uoBZq3zTom3ceQDwFTGj/IXesG4taKLRrpVyzrF6617m5mxnzuoCvs7dzp6DvhHDPVo1Zni3VEZ0T2Vo1+Z6DUVCkAq/WlLhJxLaSsscWRt3MXtlAbNWbWPhhl2UljkaJ8YxrFuL77tvO6c28jpqWCotcyzetPv7buHv1u7gYHEZMQb926d83y08WCOGRUKCCr9aUuEnEnrydh34fvTtl6sK2HOwBDMY0K7p9923gzqmEK8Rq0FXVFLKgvW7mLva1y2ctWEXJWWOhLgYLspoz91j+6gAFPGQBneISMRwzjFpVi5//nA5zkGrJomc3rc1o3qkMaJ7Ks0bHXlKNAmOxLhYju/aguO7tuDn+Kar+27NDj5eupmXv17PvHW7eOKSQXRNa+x1VBGpQC1+1aAWP5HQUFxaxu+mLObVbzcwtn8bbjslnfSWuv1IKPl8+VZ+9vpCSkodfz6/P+MGtPU6kkjUqarFT30gIhIW9hws5qoXvuPVbzdw08ndeOzHg+jRKllFX4g5uVdLpt86kh6tGnPzKwu4593FFJWUeh1LRPxU+IlIyNuwYz/nPzmXr3K285cLBvDL03uF/MwZ0axdSgNeu34Y147swktfr+P8p+ayfvt+r2OJCCr8RCTELVi/kwlPzmHLnoO8eNVxXJTRwetIUg3xsTH8dmwfJk0czPrt+xn72Gw+XLzZ61giUU+Fn4iErPcX5XPxpK9pkBDL2zeOYHj3VK8jSYBO69ua6beOpGtqI254eR73TV3CoZIyr2OJRC0VfiIScpxzPD0zhxv/PZ++bZvw7o0j6N5SI0TDVYfmDXn9hmFcObwzz89Zy4XPfMXGner6FfGCCj8RCSnFpWXc9fYi/vzBcsYNaMMr1x5Pi8aaJizcJcbFcu/ZfXny0mPJ3bqXsY9+yYylW7yOJRJ1QrLwM7MLzWyJmZWZWUa55aea2TwzW+T/Orrcui/MbIWZLfQ/Wh7h2HeZ2Wr/tqfXx/mISPXsPlDMT57/jv98t4GbT+7OoxcP0o2AI8xZ/dsw9ZYTaJfSgGtezOSP7y+juFRdvyL1JVRv4LwYOA94psLyAmC8cy7PzPoBHwHtyq2/1Dl3xBvumVkf4GKgL9AWmGFmPZxzuteAiMc27NjPVS98x9rt+/jrBQO4UIM4Ilbn1Ea8feNwHpi+lEmzcpm3bieP/XgQbVMaeB1NJOKFZIufc26Zc25FJcsXOOfy/E+XAElmFkgf0DnAf5xzRc65NcBq4LjaJxaR2vjfkbtDVfRFgaT4WB44tz+P/ngQy/P3MPbR2Xy+YqvXsUQiXkgWftV0PrDAOVdUbtnz/m7ee6zyu7q2AzaUe76R/20xFJF6dnjkbsOEON6+cQTDurXwOpLUo7MHtuW9W06gVZMkfvL8d/zlw+WUqOtXpM54VviZ2QwzW1zJ45xq7NsXeAi4vtziS51z/YGR/sfEynatZFmlc9aZ2XVmlmlmmdu2bTv6CYlIQJxzPPWFb+Ruv3ZNeefG4Rq5G6W6pTXm3ZtGcPGQDjz5RQ6XTP6GLXsOeh1LJCJ5Vvg5505xzvWr5DGlqv3MrD3wDnC5cy6n3PE2+b8WAq9QeRfuRqB8H1J7IK+S7XDOTXLOZTjnMtLS0gI7ORGpUnFpGb9+axEPfbic8QPb8u9rhmrkbpRLio/lz+cP4G8XDWTRxt2c9chsvlxV4HUskYgTVl29ZpYCTAfucs7NKbc8zsxS/d/HA+PwDRCp6D3gYjNLNLMuQDrwbZ0HF5Hv7T5QzJXPf8trmRu4ZXR3HvnRMRq5K98779j2vHfzCJo3SmDic9/w909WUlpWaceMiNRASBZ+ZjbBzDYCw4DpZvaRf9XNQHfgngq3bUkEPjKzbGAhsAl41n+ss83sfgDn3BLgdWAp8CFwk0b0itSfDTv2c/5Tc/l2zQ4evnAgd5zWU3Puyg+kt0pmys0jmDCoHY98uorLn/uGbYVFR99RRI7KnNP/pI4mIyPDZWYe8S4xIlIN89fv5Np/ZVJS5nj6ssEaxCFH5ZzjjcyN3DNlMU0axPPoxYP0vhGpBjOb55zLqGxdSLb4iUhkmZ6dz48nfU3jpDjevnG4/nhLtZgZFw3pwLs3jSA5MY5LJ3/N45+tokxdvyI1psJPROqMc44nPl/NTa/Mp3+7prxz4wi6pWnkrgSmd5smvHfLCYwb0JaHP17JlS98x/a96voVqQkVfiJSJ4pLy/jVW9n89aMVnD2wLS9fM5TmjRK8jiVhqnFiHI9cfAwPTujH17nbGfvol3y3dofXsUTCjgo/EQm63QeKueK5b3k9cyO3ju7OIxdr5K7Unplx6dBOvP3T4STGx3DxpK95emaOun5FAqDCT0SC6vDI3e/W7uD/LhzIz0/rSeUT6YjUTL92TZl6ywmc3rcVf/5gOde8mMnOfYe8jiUSFlT4iUjQzFu3k3OfmMO2wiJeunoo5w9u73UkiVBNkuJ54pJjue/svsxetY2xj85m/vqdXscSCXkq/EQkKKZl5/HjZ/87cvf4rhq5K3XLzLhieGfe+ulwYmKMi57+ismzc9FtykSOTIWfiNTK4ZG7N7+ygIHtNXJX6t+A9ilMv2UkJ/dqyQPTl3H9S/PYvb/Y61giIUmFn4jU2KGSMu580zdy95xjNHJXvNO0YTyTJg7m7rG9+Wz5VsY+Npvsjbu8jiUSclT4iUiN7N7vm3P3jXkbuXVMOv/40TEkxmnkrnjHzLhmZFdev2EYZWWOC576in/NXauuX5FyVPiJSMDWb9/PeU/N+e/I3VN7aOSuhIxjOzZj+q0jOSE9ld+/t4SbX1nAnoPq+hUBFX4iEqB563Yy4ck5FOw9pJG7ErKaNUpg8uUZ/PrMXny4ZDNnP/YlS/J2ex1LxHMq/ESk2qZm+UbuJifF8Y5G7kqIi4kxbjixG/+57ngOFJcy4cm5/Pubder6laimwk9EjurwyN1bXvWN3H37xhF01chdCRNDOjfn/VtHMrRLc377zmJuf20h+4pKvI4l4gkVfiJSpUMlZfzSP3L3XI3clTDVonEi//rJcdxxag+mZuUx/vEvWb55j9exROqdCj8ROaLd+31z7r45byO3jUnn7xq5K2EsJsa4ZUw6L18zlMKDJZz7xBxez9zgdSyReqXCT0QqdXjkbua6HfztooH8TCN3JUIM75bK9FtPYFCHZtz5ZjZ3vJ7F/kPq+pXooMJPRH5g3rodnPvkHLbvO8TLVw/lvGM1clciS8vkJF6+Zii3jknn7QUbOefxOazaUuh1LJE6p8JPRP7He1l5/PjZb2iSFMc7N45gqEbuSoSKjTF+fmoPXrzqOHbsO8TZj8/h7fkbvY4lUqdU+IkI4Bu5+/hnq7j11QUc0z6Fd24cQZfURl7HEqlzI9PTeP+2kfRv35Sfv57Fr97M5mBxqdexROqECj8R4VBJGb94I5uHP17JhEHteOma42imkbsSRVo1SeKVa4Zy40ndeC1zAxOenMuu/Ye8jiUSdCr8RKLc7v3FXP7cN7w1fyO3n5LO3y4aqJG7EpXiYmO484xePHdlBjlb93LTK/MpLi3zOpZIUKnwE4li67bvY8JTc5i/bhd//9FAbj9FI3dFRvdqxYMT+jFn9XYemLbU6zgiQRXndQAR8Ubm2h1c99I8ypzj5WuGclyX5l5HEgkZF2Z0YOWWQp6dvYYerZO5dGgnryOJBIUKP5EoNGXhJn75ZjbtUhrw3JVDNIhDpBK/PrM3q7bu5fdTltA1tTHDummEu4Q/dfWKRBHnHI99uorb/rOQYzqk8PZPh6voEzmC2Bjj0R8PolOLhtz473ms377f60gitabCTyRKlJY5fvlmNv/3iX/k7tUauStyNE2S4pl8xRDKHFzz4ncUHiz2OpJIrajwE4kSk2fn8ua8jdw6RiN3RQLRJbURT156LDnb9vGz1xZSWua8jiRSYyr8RKLAqi2F/N8nKzm9byt+dkq6Ru6KBGhE91R+N64PM5Zt5eGPV3gdR6TGNLhDJMKVlJbxizeyaJwYx4MT+qvoE6mhy4d1YsWWQp76IocerRozYZDmsJbwoxY/kQj3zKxcsjbu5g/n9CO1caLXcUTClplx39l9GdqlOb96axEL1u/0OpJIwFT4iUSw5Zv38I8ZKxk7oA1jB7TxOo5I2IuPjeGpywbTqkki1700j/zdB7yOJBIQFX4iEaq4tIw7Xs+iaYN4/nBOP6/jiESM5o0SmHz5EPYXlXDdi/M4cKjU60gi1abCTyRCPfl5Dkvy9vDAuf1prtu2iARVz9bJPHLxIBbn7ebOt7JxTiN9JTyo8BOJQEvydvPYZ6s455i2nNGvtddxRCLSKX1a8cvTezI1K48nPl/tdRyRatGoXpEIc6jE18XbrFEC947v63UckYj20xO7sXJzIQ9/vJLuLZP1Hy0JeWrxE4kwj3+2iuWbC/njhP6amUOkjpkZfz5/AAM7pPDz1xeyLH+P15FEqqTCTySCLNq4mye+yOG8Y9txap9WXscRiQpJ8bE8O3EwyUlxXPOvTAr2FnkdSeSIVPiJRIiiklLueGMhqY0T+P04dfGK1KeWTZJ49vIMCvYW8dOX53GopMzrSCKVUuEnEiEembGKlVv28ufzB9C0YbzXcUSizoD2KTx84UC+W7uTe95drJG+EpI0uEMkAixYv5OnZ+bwo4wOnNyzpddxRKLW+IFtWbmlkMc+W03P1slcdUIXryOJ/A+1+ImEuYPFpfzijSxaN0nit+N6ex1HJOr97JQenN63FQ9MX8rMldu8jiPyP1T4iYS5v32ykpxt+3joggE0SVIXr4jXYmKMv110DD1aJXPzK/PJ2bbX60gi31PhJxLG5q3bwbOzc7lkaEdGpqd5HUdE/BolxjH5igwSYmO49l+Z7N5f7HUkEUCFn0jYOnColF+8kU3bpg34zVnq4hUJNe2bNeTpiYPZsHM/N786n5JSjfQV76nwEwlTf/1oBWsK9vHXCwbQOFHjtERC0ZDOzXng3H7MXlXAg+8v8zqOiEb1ioSjb3K38/zcNVw+rBPDu6d6HUdEqvCjIR1ZsXkvz81ZQ89WyVx8XEevI0kUO2LhZ2ZlQMA3IXLOxdYqkYhUaf+hEn75ZjYdmjXkV2f08jqOiFTDb87qxepte7lnymK6pjXmuC7NvY4kUaqqFr/7+WHhdy7QD/gIWAEY0BM4DVgETAl+RBEp76EPlrN+x35eu+54GqmLVyQsxMXG8NiPBzHhiTnc8PI8ptw0gg7NG3odS6LQEf9qOOfuLf/czK4CWgP9nHMrKqzrDXwOrK+DjCLiN3d1Af/6ah0/GdGZoV1beB1HRALQtEE8k6/I4Nwn5nDti5m8+dPhuj5X6l0ggzvuBB6vWPQBOOeWAU8AvwpWMBH5X3uLfF28XVIbcefp6uIVCUdd0xrz+CXHsnJLIT97bSFlZZrWTepXIIVfJ+BAFev3+7cRkTrwx/eXkbf7AA9fOIAGCbqUViRcjeqRxt1j+/DJ0i387ZOVXseRKBNI4bcSuM7MUiquMLNmwHX4rvsTkSCbtXIbr3yznmtHdmVwJ10ULhLufjKiMxcP6cDjn69mysJNXseRKBLIxQW/Ad4FVpnZS/gKQQf0Ai4DUvAN/hCRINpzsJhfv5VNt7RG/PzUHl7HEZEgMDPuP6cfudv2ceeb2XRu0YiBHVK8jiVRoNotfs656cDp+AZw3A48CTwF3OZfdqZ/GxEJogenLWPznoM8fOFAkuLVxSsSKRLiYnjqsmNJbZzIdS9lsmXPQa8jSRQIaOYO59xnzrnBQBtgGDAcaOOcG+ycm1EXAUWi2ecrtvJa5gauP7Ebgzo28zqOiARZi8aJTL4ig8KDJVz3YiYHi0u9jiQRrkZTtjnntjjnvnHOfe2c2xLsUCICu/f7unh7tGrM7aekex1HROpI7zZN+PuPjiFr425+9VY2zmmkr9SdgAo/M0sxsz+b2WIz22dme/3f/7GyQR8iUnP3T1tKwd5DPHzhQBLj1MUrEslO79uaX5zWgykL83hqZo7XcSSCVbvwM7N2wAJ89/MDmA58gG+Ax6+B+WbWNugJRaLQJ0u38Nb8jdx4UjcGtE/xOo6I1IObTu7O+IFt+etHK/hkqTrTpG4E0uL3J6AVMM451885d5Fz7kLnXH9grH/dn+oipEg02bnvEL95ZxG9Widzy2h18YpECzPjrxcMoH+7ptz+nwUs37zH60gSgQIp/M4AHnHOvV9xhXPuA+Ax4MxgBROJVvdOXcLOfYf4v4sGkhBXo8twRSRMJcXHMmliBo0S47jmX5ls31vkdSSJMIH8VUkGqrrL5Eb/NiJSQx8uzmfKwjxuGZ1O37ZNvY4jIh5o3TSJSZdnsLWwiJ/+ez6HSsq8jiQRJJDCbwVwgZn9YB8ziwUuQDN3iNTY9r1F/PadxfRt24QbT+7mdRwR8dAxHVL46wUD+HbNDn7/3hKN9JWgCWTmjkeBycBnZvZ3/lvk9cJ3E+eRwNXBjScSPX733hL2HCzm39cOJT5WXbwi0e6cY9qxYnMhT36RQ6/WyVwxvLPXkSQCVLvwc849Z2Ytgd8Db5dbZUAR8Bvn3AvBjScSHaZl5zE9O59fnt6TXq2beB1HRELEL07rycote7l/2lK6pTXmhPRUryNJmLNAm4/NrAVwKtDJv2gt8Ilzbkdwo4WOjIwMl5mZ6XUMiVDbCos47e8z6dC8IW//dDhxau0TkXL2FpVw/pNzyd99gCk3n0CX1EZeR5IQZ2bznHMZla0L+C+Mc267c+4/zrmH/I/Xgl30mdmFZrbEzMrMLKPc8lPNbJ6ZLfJ/HV1u3RdmtsLMFvofLSs5bgsz+9x/4+nHg5lZpCacc9z97iL2FZXyfxcOVNEnIj/QODGOyVdkEBtjXP2v79h9oNjrSBLGAv4rY2ZnmNnjZjbdzKb5vz8tyLkWA+cBsyosLwDG++8deAXwUoX1lzrnjvE/tlZy3IPAPcAvgpxXpEbey8rjoyVb+PlpPUhvpUHxIlK5Ds0b8tRlg1m/fT+3vrqA0jIN9pCaCWTmjgQzm4Jvxo4bgSHAUP/3H5jZu2aWEIxQzrllzrkfjBB2zi1wzuX5ny4BkswsMYDj7nPOfYmvABTx1NY9B/ndlCUM6pjCtSO7eh1HRELc8V1bcP85/Zi5cht/en+Z13EkTAXS4vd7YDzwf0BL51xL51wakAY8DJyNrzWtvpwPLHDOlb+75fP+bt57zMzqMYtIQJxz/OadRRwsLuXhCwcSG6O3q4gc3SVDO3Ll8M5M/nINr3+3wes4EoYCKfwuAV52zt3pnCs4vNB/zd+vgJeBy6p7MDObYWaLK3mcU419+wIPAdeXW3ypvwt4pP8xsbpZjvAzrjOzTDPL3LZtW20OJfIDb8/fxIxlW/nl6T3pltbY6zgiEkbuHtubE7qn8tt3F5G5NmLHVUodCaTwawvMrWL9V0Cb6h7MOXeKf87fio8pVe1nZu2Bd4DLnXM55Y63yf+1EHgFOK66WY6Qb5JzLsM5l5GWllabQ4n8j827D3Lv1CUM6dyMn4zo4nUcEQkzcbExPH7JINqlNOCGl+exced+ryNJGAmk8MsDjq9i/VAgv3ZxqmZmKfiuMbzLOTen3PI4M0v1fx8PjMM3QEQkpDjn+PXb2RSXlvHXC9TFKyI1k9IwgclXDKGopIxrX5zHvqISryNJmAik8HsVmGhmD5hZs8MLzayZmf0BX9fqK8EIZWYTzGwjMAyYbmYf+VfdDHQH7qlw25ZE4CMzywYW4ptT+Fn/sc42s/vLHXst8DfgSjPbaGZ9gpFZpDreyNzIFyu28eszetFZ9+ISkVro3rIxj/14ECs27+GO17Mo00hfqYZq38DZP3r2HeAMwAGHL3xLwzd7x4fAhAqDLSKCbuAswbBp1wHO+Pss+rRtwqvXHk+MWvtEJAgmz87lgenLuHVMOj8/tYfXcSQEVHUD50CmbCsCzjKzccBYoLN/1VpgqnPu/VrmFIlYzjl+/VY2pc7x1wsGqugTkaC5+oQurNhcyKOfrqJHq8aMG9DW60gSwqpd+B3mnJsGTKuDLCIR69VvNzB7VQF/OLcfHVs09DqOiEQQM+OBCf3ILdjHL97IonOLRvRr19TrWBKiND+USB3bsGM/D05fyojuLbj0uI5exxGRCJQYF8vTlw2mecMErn0xk617NE+BVC6gws8/V+5rZvadmeWYWW6FR87RjyISPcrKHHe+mQ3AQ+cPUBeviNSZtOREnr0ig137i7nupXkcLC71OpKEoECmbPsZvgEcJ+IbNTsLmFnhUXFuXZGo9vI36/gqdzt3j+tD+2bq4hWRutW3bVP+/qOBLNywi9+8vYjqDuCU6BHINX4/w1fcneGcO1RHeUQixrrt+/jT+8sZmZ7KxUM6eB1HRKLEGf3a8PNTe/C3T1bSs3Uy15/YzetIEkIC6epNBV5T0SdydGVljl++kU1cjPHQ+QPQ1NEiUp9uGd2dsQPa8OcPlzN7laYdlf8KpPCbD3StqyAikeSFuWv5du0Ofje+D21TGngdR0SijJnx8AUD6ZLaiN++s1jX+8n3Ain8bsc3c8epdZRFJCLkbtvLXz5azuheLblgcHuv44hIlGqQEMsD5/Rj/Y79PPn5aq/jSIg44jV+ZvZxJYsLgQ/9056tBSr+F8I5504PWjqRMFNa5vjlm9kkxMbwp/P6q4tXRDw1vHsq5x7Tlqdn5nLuoHZ0TWvsdSTxWFUtfj2A9AqPBGC9f7+ulazXXDES1Z77cg3z1u3kvnP60qpJktdxRET4zdjeJMbHcM+UxRrlK0du8XPOda7HHCJhb/XWQv768QpO7dOKc49p53UcEREAWiYn8cvTe/K7KUt4LyuPc/T5FNU0c4dIEJSUlnHHG9k0TIjlwQn91MUrIiHl0qGdGNC+KQ9MX8aeg8VexxEPqfATCYJJs3PJ2rCL+8/pR8tkdfGKSGiJjTEeOLcfBXuL+NvHK72OIx46YuFnZmVmVmJmCeWelx7lUVJ/0UVCw4rNhfzjk1Wc2a814we08TqOiEilBrRPYeLxnXjxq7Us3rTb6zjikapm7rgfcEBJheci4ldcWsYv3siicVIcfzhXXbwiEtruOK0n7y/azG/fWcTbN44gVvOHR52qBnfcW9VzEYGnv8hh0abdPHnpsaQ2TvQ6johIlZo2iOfusb25/bWFvPLteiYe38nrSFLPdI2fSA0tzdvDo5+tYvzAtpzVX128IhIezjmmLcO7teAvHy5nW2GR13GknlV1A+dRNTmgc25WzeOIhIdDJb4u3qYNErj/7L5exxERqTYz4/5z+nHmI7P40/vL+NuPjvE6ktSjqq7x+4LArukz//axtQkkEg6e+Hw1S/P3MGniYJo1SvA6johIQLq3bMz1o7rx+OeruTCjA8O6tfA6ktSTqgq/k+sthUgYWbxpN098vpoJg9pxWt/WXscREamRm07uzrsLN3HPlMW8f+tIEuJ09Vc0qGpwx8z6DCISDopKSrnj9SyaN0rg9+P7eB1HRKTGGiTEcv85fbnqhUyenZ3LTSd39zqS1IMalfdmlm5mI8ysabADiYSyRz9dxYothfzpvP6kNFQXr4iEt9G9WnF631Y89tkqNuzY73UcqQcBFX5m9iMzWwcsB2YBg/3LU81slZldWAcZRUJC1oZdPPVFDhcMbs+Y3q28jiMiEhS/H9+XGDPum7rE6yhSD6pd+JnZOcCrwHrgHnyDOQBwzhUAy4CJwQ4oEgoOFpdyxxtZtExO4p5x6uIVkcjRNqUBt5+SzoxlW/l4yWav40gdC6TF725glnNuJPBMJeu/AQYGJZVIiPn7jJWs3rqXP5/fn6YN4r2OIyISVD8Z0YWerZK5b+pS9h/S7KuRLJDCry/wehXrNwPq/5KIM2/dTp6dlcvFQzpwUs+WXscREQm6+NgYHpjQj027DvDIp6u8jiN1KJDC7yCQVMX6zsCu2oQRCTUHi0v55RtZtGnagN+O7e11HBGROjOkc3MuHNyef85ew4rNhV7HkToSSOH3JfDjylb4R/deBXwWjFAioeLVb9eTW7CPP53Xn+QkdfGKSGS766zeNE6K4553F+NcIHM4SLgIpPC7F+hrZp8D5/mXZZjZzcBCoAnwh6CmE/FQcWkZk2evIaNTM0b1SPM6johInWveKIFfn9GLb9fu4K35m7yOI3Wg2oWfc24+cDrQmv8O7vgz8ChwCDjdObcs6AlFPPL+onw27TrA9Sd28zqKiEi9uSijA8d2TOGP7y9j1/5DXseRIAvoPn7OudnOud7AIOBH+Lp+hwC9nHNzzaxxHWQUqXfOOZ6emUu3tEaM6aUBHSISPWJijAfO7c/uA8U89OEKr+NIkAVyH7/HD3/vnMtyzr3hnHvNOTfPOefMLAWYURchRerb7FUFLMvfw/WjuhETY0ffQUQkgvRp24Qrh3fm1W/XM3/9Tq/jSBAF0uJ3vZn9sbIVZtYc38COnkFJJeKxSbNyaZmcyDmD2nodRUTEEz87tQetmyRx9zuLKSkt8zqOBEkghd/VwK/M7K7yC82sJb7p2zoBpwUxm4gnFm/azZerC7jqhC4kxsV6HUdExBONE+P43fg+LM3fw4tfrfM6jgRJIIM7XgRuBR40s5sAzKw9MBtoCYx2zn1XJylF6tEzs3JpnBjHJUM7eh1FRMRTZ/ZrzYk90vjbJyvZsueg13EkCAId3PEE8FvgETP7NTATSAZOcs5l1UE+kXq1Ycd+pmfnccnQjjTRfftEJMqZGfed3ZdDpWXcP22p13EkCAIq/ACcc38C/gL8EYgDRjrn9G6QiDB5di6xMcZPRnT2OoqISEjonNqIm07qzvTsfGat3OZ1HKmluCOtMLNJR9m3EFiO77q/w8ucc+76IGUTqVc79h3itcwNnHNMO9o0beB1HBGRkHHDSV15d+EmfjdlMR/ePoqkeF3/HK6OWPgB11Rj/1MrPHeACj8JSy99tY6DxWVcN6qr11FEREJKYlwsfzinH5f98xuenpnD7af08DqS1NARu3qdczE1eOi/ABKWDhwq5V9frWVMr5b0aJXsdRwRkZBzQnoq4we25ckvclhbsM/rOFJDAV/jJxKJ3py3gR37Dqm1T0SkCneP7U1CbAz3TFmMc87rOFIDKvwk6pWUlvHs7DUc0yGF47o09zqOiEjIatUkiTtO68HsVQVMX5TvdRypgaoGd6wByvDNw1vsf3608t455zSjvYSVD5dsZv2O/fzmrF6UG6gkIiKVmHh8J96ct5H7py7lxB5pJOvWV2GlqsEdM/EVemUVnotEDOccz8zMpUtqI07t09rrOCIiIS8uNoYHJ/RnwpNz+Psnq/jd+D5eR5IAHLHwc85dWdVzkUjwVe52Fm3azR8n9Cc2Rq19IiLVcUyHFC45riMvzF3D+YPb0bdtU68jSTXpGj+Jas/MzCW1cQLnHdvO6ygiImHlztN70axhAne/u5iyMnUIhouqrvGr0USlzrn1NY8jUn+W5e9h5spt/OK0HroZqYhIgJo2jOc3Z/Xmjjey+M93GzS/eZio6hq/tdTsmj79BZWwMGlWLg0TYrns+E5eRxERCUvnHduO1zM38NCHyzm9bytaNE70OpIcRVWF31VoMIdEqE27DvBeVh5XDOtMSsMEr+OIiIQlM+OBc/tx5iOz+dMHy3n4woFeR5KjqGpwxwv1mEOkXj335RoArh7ZxeMkIiLhLb1VMteO6spTX+RwUUYH3Q81xGlwh0Sd3fuLefXb9Zw9sC3tUhp4HUdEJOzdMro77VIacPe7iyguLTv6DuIZFX4SdV7+Zh37D5Vy7UhNzyYiEgwNE+K49+y+rNyyl3/6e1QkNKnwk6hysLiU5+esYVSPNPq0beJ1HBGRiHFqn1ac0rsVj8xYxaZdB7yOI0egwk+iytvzN1Gw9xA3jFJrn4hIsN17tm8Wj/veW+JxEjkSFX4SNUrLHM/OzqV/u6YM69bC6zgiIhGnfbOG3DomnY+XbuHTZVu8jiOVUOEnUeOTpVtYU7CP60/sipmmZxMRqQtXn9CF9JaN+f17SzhwqNTrOFJBtQs/M+t4lEcHM0sz/UWVEOSc4+mZOXRs3pAz+rb2Oo6ISMRKiIvhD+f2Y+POAzz22Sqv40gFgbT4rQXWVPFYC2wG9prZ+2Z2fFCTitTCd2t3snDDLq4Z2YW4WDV0i4jUpeO7tuC8Y9vx7OxcVm8t9DqOlBPIX8CrgSxgF/AEcDvwM+Ap/7IFwG3AZGAoMNPMRgYvqkjNPTMzh2YN47lwcAevo4iIRIXfnNWbhglx3P3uYpzTRGChIpDCrzXQEEh3zt3qnHvMOfeoc+5moCeQDCQ6524DegFbgd8HPbFIgFZtKeTT5Vu5YnhnGiRoKmkRkfqQ2jiRO8/oyde5O3h34Sav44hfIIXfDcCzzrkdFVc45wrwtfTd7H++DfgnMCQYIUVqY9KsXJLiY7h8WGevo4iIRJUfD+nIwA4pPDh9Gbv3F3sdRwis8GsJxFexPg5fq+BhG6liLmCR+rB590HeXbiJH2V0oHmjBK/jiIhElZgY48Fz+7Fj3yH++vFyr+MIgRV+2cBNZtax4goz6wTchO8awMPSgfzaxROpnefnrKG0zHGNpmcTEfFEv3ZNuXxYZ/79zXqyNuzyOk7UC6TwuwNoDqwwszfN7E/+x5vAcqAF8AsAM0sCLgU+DXZgkerac7CYf3+znrP6t6FD84ZexxERiVp3nNaDtMaJ3P3uYkrLNNDDS9Uu/JxzXwLDgY+BM4Ff+R9nAh8Bx/u3wTl30DnX3jl3ffAji1TPK9+sZ29RCdeP6uZ1FBGRqJacFM894/qwaNNuXv56nddxolpANzRzzmU5587BN4K3rf+R7Jw71zmXVfXe1WdmF5rZEjMrM7OMcstPNbN5ZrbI/3V0uXVfmNkKM1vof7Ss5LhH3F8iS1FJKc99uYYR3VvQv31Tr+OIiES9cQPaMDI9lYc/WsHWPQe9jhO1anQnW+dcGVAKlPi/D7bFwHnArArLC4Dxzrn+wBXASxXWX+qcO8b/2FrJcY+2v0SIKQvz2FpYpNY+EZEQYWbcd3ZfikrKeGD6Mq/jRK2ACj8z62pm/zGz3fhm6dhiZrvN7BUzC9rV8865Zc65FZUsX+Ccy/M/XQIkmVliAMet1f4SHsrKHJNm5dK7TRNGpqd6HUdERPy6pjXmhpO68V5WHnNWF3gdJyoFMldvL+A74AJgNvB/wN/8318IfOvfpr6cDyxwzhWVW/a8v5v3nmrMGVzZ/hIBPlu+ldVb93LDiV3R1NEiIqHlxpO60alFQ+55dzFFJaVex4k6gbT4/RkoAwY558Y55+50zv3SOTcOGAQ44I/VPZiZzTCzxZU8zqnGvn2Bh4Dyg0cu9XfhjvQ/Jga4f8VtrjOzTDPL3LZtW3VPS0LAM7NyaJfSgLP6t/E6ioiIVJAUH8v95/Qjt2Afk2bmeh0n6gRS+J0IPOacW1RxhXNuMfA4cHJ1D+acO8U516+Sx5Sq9jOz9sA7wOXOuZxyx9vk/1oIvAIcF8j+leSb5JzLcM5lpKWlVfe0xGPz1u3ku7U7ufqELsTH1ugSVhERqWMn9khjbP82PP75atZv3+91nKgSyF/GBGBPFet3+7epM2aWAkwH7nLOzSm3PM7MUv3fxwPj8A0Qqdb+EjkmzcqhaYN4fjSkg9dRRESkCveM60NcjPG79xbjnO7tV18CnbnjCjNrUHGFf9kV/m1qzcwmmNlGYBgw3cw+8q+6GegO3FPhti2JwEdmlg0sBDYBz/qPdbaZ3X+U/SUC5Gzby8dLt3D5sE40StRsgSIioax10yR+dmoPvlixjY+WbPY6TtSw6lbZZjYeeBdYBTwFHB512wu4AV9Bda5zblrwY3orIyPDZWZmeh1DjuKut7N5a/4m5v56NKmNNVhbRCTUlZSWMf7xOezaf4gZPz9R/2kPEjOb55zLqGxdIDN3TAUuA5oAf8fXZTod38jeJsBlkVj0SXjYWniQt+Zt4oLB7VX0iYiEibjYGB44tx/5uw/yjxkrvY4TFQIqrZ1zr5rZG8BgoLN/8Vog0zmnMdnimRfmrKW4rIxrRwbtdpIiIlIPBndqxo+P68Bzc9Zy3rHt6d2mideRIlrAwx6dcyXOuW+cc6/5H9+o6BMv7S0q4aWv13FG39Z0SW3kdRwREQnQnaf3ommDeO5+dzFlZRroUZeO2OJnZh1rckDn3PqaxxEJ3H++XU/hwRKuG6XWPhGRcNSsUQK/PrMXd76ZzZvzNnKR7sxQZ6rq6l2L76bMgYqtWRSRwBWXlvHPL9cwtEtzBnVs5nUcERGpoQuObc8bmRv40wfLOLVPK5o1qtM7xEWtqgq/q6hZ4SdSb6Zm5ZG/+yB/nNDf6ygiIlILMTHGA+f2Z+yjs/nzB8t56IIBXkeKSEcs/JxzL9RjDpGAOed4ZmYuPVo15qSeml1FRCTc9WydzNUndOGZWblcmNGejM7NvY4UcTSnlYStL1ZuY8WWQq4b1Q0z8zqOiIgEwa1j0mnbNIm7311MSWmZ13Eijgo/CVvPzMyhdZMkzh7Y1usoIiISJI0S4/jd+L4s31zIC3PXeh0n4qjwk7CUtWEXX+fu4OoTupAQp7exiEgkOb1vK0b3asnfP1lJ/u4DXseJKPqLKWFp0qxckpPiuPg4DfkXEYk0ZsZ9Z/el1Dnun7rU6zgRRYWfhJ21Bfv4YHE+lx3fieSkeK/jiIhIHejQvCG3jE7ng8Wb+XzFVq/jRAwVfhJ2Jn+ZS1xMDD8Z3tnrKCIiUoeuGdmFrmmN+P2UJRws1iRhwRBQ4WdmCWZ2tZn928w+MbNB/uXNzOxyM2tfNzFFfAr2FvFG5kYmDGpHyyZJXscREZE6lBgXywPn9GP9jv08PTPH6zgRodqFn5k1A74GngXOAkYDh6dK2A38Abg52AFFynvxq3UUlZRxraZnExGJCsO7p3Jmv9b8c/Yadu8v9jpO2Aukxe/PQA/gNCAd+P7Gac65MuBt4IygphMpZ/+hEl78ai2n9mlF95aNvY4jIiL15NYx6RQWlfDcnDVeRwl7gRR+ZwOPOudmUPlUbquBTkFJJVKJ17/bwK79xdxwolr7RESiSe82TTi9byuem7OG3QfU6lcbgRR+zYDcKtbHAppRWepESWkZz85eQ0anZgzupCl8RESiza1j0ik8WMK/dFPnWgmk8FsD9K9i/YnAitrFEanc9EX5bNp1gOt0bZ+ISFTq27Ypp/RuxT+/XEPhQbX61VQghd/LwDVmNrrcMgdgZj8DzgVeCFoyET/nHM/MzKVrWiNO6d3K6zgiIuKR28aks/tAsVr9aiGQwu8h4FPgE3yjex3wuJltAf4PmAI8FvSEEvXmrN7O0vw9XD+qKzExdvQdREQkIvVv35TRvVoy+cs17C0q8TpOWKp24eecK3HOnQ1cBiwGlvv3/xaY6Jw7zzlX2aAPkVp5ZlYOacmJnDuonddRRETEY7eNSWfX/mJe/Gqt11HCUlygOzjnXgVerYMsIj+weNNuZq8q4Fdn9CIxLtbrOCIi4rGBHVI4qWcak2ev4YphnWmUGHApE9U0ZZuEtEmzcmmUEMslQzt6HUVERELErWPS2bHvEC9/vc7rKGHniGWymT1Xg+M559zVtcgj8r0NO/YzfVE+V43oTNMG8V7HERGREHFsx2aMTE9l0qxcJg7rRMMEtfpVV1Wv1Ggqv1FzVXSNnwTNP79cgwFXndDF6ygiIhJibj8lnfOf+opXvlnPNSN1q6/qOmLh55zrXI85RP7Hzn2HeO27DZxzTDvaNG3gdRwREQkxgzs1Z0T3Fjw9M5dLh3aiQYKuA68OXeMnIemlr9dxoLhUN2wWEZEjum1MDwr2FvHqt+u9jhI2VPhJyDlYXMoLc9cyuldLerZO9jqOiIiEqOO6NOf4rs15emYOB4tLvY4TFo5Y+JnZGjPLMbP4cs9zj/LIqb/oEqnemLeRHfsOqbVPRESO6rYxPdhaWMR/1OpXLVUN7piJb7BGWYXnInWmtMwxeXYuAzukMLRLc6/jiIhIiBvWrQXHdWnOUzNzuPi4jiTF61q/qlRV+N0K7HfOlQI4566sl0QS1T5aspl12/fz6zN6Yabp2URE5OhuG5POpZO/4Y3MDUwc1tnrOCGtqmv8dgIXHn5iZs+Z2dC6jyTRyjnHMzNz6NyiIaf1be11HBERCRPDu7Ugo1Mznvwih6ISXetXlaoKv2IgsdzzK4FudZpGotrXuTvI2riba0d1JTZGrX0iIlI9ZsatY9LJ332QN+dt9DpOSKuqq3c1cKmZzQN2+5e1MLMq585yzunqSqmRZ2blkNo4gfOPbe91FBERCTMj01MZ1DGFJz/P4cLBHUiI041LKlPVq3I/cCKwEFiDb2DHP/zfV/UQCdjyzXv4YsU2rhjWWRfmiohIwMyM28aks2nXAd6er1a/I6lq5o7X/a19JwGtgAeAN4Gs+okm0WTSrFwaxMcycVgnr6OIiEiYOrFHGgPbN+Xxz1dz/uD2xMeq1a+iKmc1ds7lADkAZnYt8LJz7r36CCbRI2/XAd5bmMfEYZ1IaZjgdRwREQlTZsZtp6Rz1QuZvDN/ExcN6eB1pJBT7VLYOddFRZ/Uhee+XIMDrj6hi9dRREQkzJ3csyX92/la/UpKy46+Q5QJqA3UzOLM7AYzm2ZmS8xssf/768ysytZDkcrs3l/Mq9+uZ/yANrRv1tDrOCIiEuYOj/Bdv2M/7y7M8zpOyKl24WdmTYGvgCeB4cA+4AAwDHgamOvfRqTaXv5mHfsOlXLdKN0pSEREguOU3i3p06YJT6jV7wcCafF7EBgE3Ay0cs4d55wbgm/gx03+dQ8EP6JEqoPFpTw/Zy0j01Pp07aJ13FERCRCHG71W1Owj6nZavUrL5DCbwLwlHPuSedc8eGFzrkS59xTwDPAecEOKJHr3QWbKNhbxA0nqrVPRESC67Q+rejVOpnHPltNaZnzOk7ICKTwawEsrWL9Ev82IkdVVuaYNCuXfu2aMLyb3jYiIhJcMTG+Vr/cbfuYpla/7wVS+K0FTq9i/Rn+bUSO6pNlW8gt2Mf1o7phpunZREQk+M7o25qerdTqV14ghd9zwNlm9rKZDTSzJP/jGDN7ERgHTK6bmBJJnHM8PTOHDs0bcGa/1l7HERGRCBUTY9wypjurt+7lg8X5XscJCYEUfn8FJgGXAPPxjerdB8wDLgMmOeceDnpCiTiZ63ayYP0urjmhK3G6q7qIiNShM/u1oXvLxjz66SrK1OoX0A2cnXPuBqA/8Bt8ReAk//f9nXM/rZuIEmmemZlDs4bxXJjR3usoIiIS4WJjjFtGd2fllr18uGSz13E8V62bLptZQ2A28Kxz7ml8AzlEArZ6ayEzlm3ltjHpNEzQPb9FRKTujRvQlkc+XcWjn67ijL6tiYmJ3mvLq9Xi55zbD3QFdBdEqZVJs3JJio/h8mGdvI4iIiJR4nCr3/LNhXy8dIvXcTwVyAVWnwMn1VEOiQJb9hzknQWbuCijAy0aJ3odR0REosj4AW3pktqIRz9dhXPRe61fIIXfrcBAM/u7mfXQ3LwSqOfmrKG0zHHNCV29jiIiIlEmLjaGm07uztL8PcxYttXrOJ4JpPBbA/TAVwAuA4rM7FCFR1GdpJSwV3iwmFe+Xs+Z/dvQsUVDr+OIiEgUOveYtnRq0ZBHPl0Zta1+gbTa/RuIzldJau3Vb9dTWFTC9aPU2iciIt443Op355vZfL5iK6N7tfI6Ur2rduHnnLuyDnNIBDtUUsY/v1zD8G4tGNA+xes4IiISxSYMasdjn63ikRmrOLlny6ibPUp3z5U6N2XhJrbsKeL6E7t5HUVERKJcfGwMN53UnayNu/li5Tav49S7gAo/M0sxswfNLMvMdvsfWf5lzeoqpISvsjLHpFm59GqdzKj0VK/jiIiIcN6x7WmX0oBHZkTfCN9qF35m1g3IBu4CYoEZwKf+7+8Csv3biHzv8xVbWbV1Lzec2C3qmtNFRCQ0JcTFcOPJ3Vi4YRezVxV4HadeBdLi9yjQDDjNOdfPOXe+c+4851w/4HQgxb+NyPeemZlL26ZJjB3QxusoIiIi37tgcHvaNk3ikSi7r18ghd+JwD+cczMqrnDOfYKv6DsxWMEk/M1fv5Nv1+7g6pFdiY/V5aQiIhI6EuNi+elJ3Zi3bidzc7Z7HafeBPLXeC9Q1VWQW/zbiAAwaWYuTRvEc/GQDl5HERER+YGLhnSgdZOkqLrWL5DC71XgMjNLqLjCzJKAicArwQom4S13214+WrqZicd3olGiJnkREZHQkxgXyw0nduXbtTv4OneH13HqRSCF31QgHlhoZrea2RlmdrqZ3QbMwzfIY6qZDS//qIvQEvqenb2G+NgYrhje2esoIiIiR3TxcR1pmZzII5+u9DpKvQikKab8tX3/4L+zeNgRtjH/NrE1SiZha1thEW/N38gFg9uTlpzodRwREZEjSoqP5YYTu3H/tKV8k7udoV1beB2pTgVS+P2kzlJIRHn563UUl5Zx7UhNzyYiIqHvkqEdefKLHB79bBX/VuHn45z7V10GkcjgnGPKwk0M79aCLqmNvI4jIiJyVL5Wv648MH0ZmWt3kNG5udeR6kxI3mPDzC40syVmVmZmGeWWn2pm88xskf/r6HLrvjCzFWa20P9oWclxjyu3PsvMJtTXOUWLxZv2sHb7fs4e2NbrKCIiItV2ydCOtGiUwCOfrvI6Sp0KycIPWAycB8yqsLwAGO+c6w9cAbxUYf2lzrlj/I+tRzhuhnPuGOAM4Bkz05DTIJqanUd8rHF639ZeRxEREam2hglxXDeqK7NXFTB//U6v49SZkCz8nHPLnHMrKlm+wDmX53+6BEgys2qPHnDO7XfOlfifJvHfASoSBGVljmlZeYxMTyOl4Q/u+iMiIhLSLju+E80bJfDIjMht9QvJwq+azgcWOOeKyi173t+Ne48dYWJYMxtqZkuARcAN5QpBqaUFG3aSt/sg4wdqejYREQk/jRLjuGZkF2au3MbCDbu8jlMnPCv8zGyGmS2u5HFONfbtCzwEXF9u8aX+LuCR/sfEyvZ1zn3jnOsLDAHu8t98urKfcZ2ZZZpZ5rZtVU1YIodNzconMS6GU3q38jqKiIhIjVw+rDMpDeN5NEKv9fOs8HPOneKc61fJY0pV+5lZe+Ad4HLnXE65423yfy3EN4PIcUf5+cuAfUC/I6yf5JzLcM5lpKWlBXZyUai0zDEtO5/RvVqSnBTvdRwREZEaaZwYxzUndOGz5VtZtHG313GCrtqFn5n1NbPzKiw72cw+9Y+wvSP48X6QIQWYDtzlnJtTbnmcmaX6v48HxuEbyFFx/y6HB3OYWSegJ7C2rnNHg29yt1Owt4hxAzSaV0REwtsVwzvTtEF8RI7wDaTF7y/ANYefmFk74D1gANAA+IuZXR6MUGY2wcw2AsOA6Wb2kX/VzUB34J4Kt21JBD4ys2xgIbAJeNZ/rLPN7H7//icAWWa2EF+r4Y3OuYJgZI52U7PzaZgQy+heP7iLjoiISFhJTorn6hO6MGPZFhZviqxWP3OuegNbzSwf+Idz7iH/8zuBe4F059wmM5sGpDnnhtZVWK9kZGS4zMxMr2OErOLSMoY8OIMTe6TxyMWDvI4jIiJSa7sPFHPCQ58xvFsLnpmYcfQdQoiZzXPOVRo6kBa/ZsCWcs/PBGYevrYOmAr0qFlECWdfri5g1/5ixqubV0REIkTTBvFcNaILHy3ZwrL8PV7HCZpACr8CoD2AmTXC1w07o9z6OAKb+1cixNSsPJKT4hjZI9XrKCIiIkFz1YguJCfG8dhnkXOtXyCF32zgp/4BHv8A4vFd43dYT3zX1kkUOVhcyidLtnBG39YkxsV6HUdERCRomjaM58oRnXl/0WZWbC70Ok5QBFL4/QY4ALwJXA38xTm3CsDMYoELgJlBTyghbebKbRQWlTBec/OKiEgEuvqELjRKiOXRCGn1q3bXrHNujZn1AvoAu51z68qtbgj8FMgKcj4JcVOz8mjeKIHh3Vp4HUVERCToUhomcMXwzjw1M4dVWwpJb5XsdaRaCegGzs65EudcdoWiD+dcoXNuinNubVDTSUjbf6iET5dt5cx+rYmLDefZ/0RERI7smpFdaRAfy2OfrfY6Sq0F9NfazJqb2R/MbI6ZrTKzYf7lLczsd/4WQYkSny7byoHiUnXziohIRGveKIHLh3VmanYeq7fu9TpOrQQyc0cHfDdHvhNIBrriu3EzzrntwI+Bm4IfUULV1Kw8WjVJZEjn5l5HERERqVPXjuxCUlwsT3we3q1+gc7ckQQcA4wGrML6Kf7lEgX2HCzmixXbOKt/G2JjKr4VREREIkuLxolMHNaJKQs3kbstfFv9Ain8TgUedc4tAyqb7mMN0CEoqSTkfbJkC4dKy9TNKyIiUePakV1JiIvhic9zvI5SY4EUfo2ArVWsb1zLLBJGpmbn0S6lAYM6pHgdRUREpF6kJSdy6dBOvLtwE+u27/M6To0EUvitAI6vYv1ZwOLaxZFwsGPfIb5cVcD4gW0xUzeviIhEj+tHdSUuxng8TEf4BlL4PQNcZmZXAYenaHBmlmxmfwdOAp4Mcj4JQR8u3kxJmWPcgDZeRxEREalXLZsk8ePjOvL2gk1s2LHf6zgBq3bh55x7CpgETAZy/YvfBHYCtwGPOedeDnpCCTnTsvPomtqIvm2beB1FRESk3v30pG7ExlhYjvAN9AbONwMj8BV/HwDfAk8BI51ztwc9nYScrXsO8lXudsapm1dERKJUqyZJXDykA2/O28jGneHV6hfwdAvOua+cc7c758Y65850zt3inJtTF+Ek9Ly/KB/nYLy6eUVEJIr99KRuxJjx5BfhNcJX82xJQKZm59OrdXLYz1UoIiJSG22aNuCiIe15I3MDm3Yd8DpOtQU6Zds1ZjbXzDabWZGZHarwKKqroOK9TbsOMG/dTt27T0REBPjpSd0BeDqMWv3iqruhmT0M/AzIA74CdtVRJglR07PzADSaV0REBGiX0oALBnfgte82cOPJ3WjTtIHXkY6q2oUfcBW+AR3nOOdK6yiPhLCpWfkMbN+UTi0aeR1FREQkJNx4UjfeyNzAMzNzuffsvl7HOapAunoNmKqiLzqtLdjHok27GTdA3bwiIiKHdWjekPOPbc8r365ny56DXsc5qkAKv4+AIXUVRELbNH8371h184qIiPyPm07uTmmZ4+mZoX+tXyCF363AEDP7vZm1q6tAEpqmZuUzpHMz2qaE/vULIiIi9alji4ZMGNSOV75Zz9bC0G71C2Tmjq3AS8DvgPVmVqxRvdFhxeZCVmwp1GheERGRI7j55O4Ul5YxaWbu0Tf2UCCjeu8HfotvVG8mGtUbNaZl5xFjcGY/dfOKiIhUpnNqI849ph0vf7OO60/sRlpyoteRKhXIqN7r0ajeqOOcY1p2PsO6tQjZN7GIiEgouHl0d95duInJs3O566zeXsepVCDX+CWhUb1RZ0neHtYU7GO8RvOKiIhUqWtaY84e2JYXv1rH9r2hefVbIIXfx2hUb9SZmpVHXIxxRr/WXkcREREJeTeP7s7BklImf7nG6yiVCqTwuxkYbGb3aVRvdDjczTsyPZWUhglexxEREQl53VsmM25AW16cu5ad+w55HecHAin8NgL9gLvRqN6oMH/9LjbtOqDRvCIiIgG4ZXR39heXMvnL0BvhG8jgjn8Drq6CSOiZmpVHQlwMp/Zp5XUUERGRsNGjVTJn9WvDv+au49qRXUOq16zahZ9z7so6zCEhprTMMX1RPif3TCM5Kd7rOCIiImHlljHdmb4on+e+XMPPT+vpdZzvBdLVK1HkmzXb2VZYpG5eERGRGujVugln9G3N83PWsvtAsddxvnfEFj8zGwXgnJtV/vnRHN5ewtu07HwaJsQyuldLr6OIiIiEpVvHpPPhks08P2cNt5/Sw+s4QNVdvV8AzswaOOcOHX5exfbmXx8btHTiieLSMj5YlM8pvVvRMCGQy0BFRETksD5tm3Ban1Y89+UarjqhC01C4NKpqv6qnwzgL/q+fy6Rb87qAnbuL2bcAE3RJiIiUhu3jknn46Vb+NectdwyJt3rOEcu/JxzM6t6LpFrWnY+yUlxnNgzzesoIiIiYa1fu6ac0rslk79cw5UjOns+YLLagzvM7DMzG1PF+pPN7LPgxBKvFJWU8tHizZzetzWJceq1FxERqa1bx6Sz+0AxL361zusoAY3qPQmo6oZuLYETa5VGPDdzxTYKi0o0mldERCRIBrRP4eSeaUyencu+ohJPswR6O5eqBnd0A/bWIouEgKnZ+TRrGM/wbi28jiIiIhIxbh2Tzs793rf6VTlk08wmAhPLLbrLzH5SyaYpwCDgk+BFk/q2/1AJM5ZuYcKx7YiP1S0eRUREgmVQx2aM6pHG4rzdnuY42r06mgOHh6A4oDWQXGEbB+zDN6Xb3UFNJ/Xqs+VbOVBcyvgB6uYVEREJtqcvO9bz26RV+dOdc48AjwCYWRlwu3PulfoIJvVvalYeLZMTOa5Lc6+jiIiIRByviz4IbK5e9f1FsMKDxXy+YhuXHNeR2BjzOo6IiIjUARVzAsAnS7dwqKRMo3lFREQimAo/AXzdvO1SGnBsxxSvo4iIiEgdUeEn7Nx3iNmrChg3oA1m6uYVERGJVCr8hA+XbKakzKmbV0REJMKp8BOmZefRJbURfds28TqKiIiI1CEVflFua+FBvsrZznh184qIiES8ahV+ZpZiZseZWbcqtuliZpcHL5rUhw8WbabMwTh184qIiES8oxZ+ZvYgsAX4ClhpZt+a2bGVbDoceD7I+aSOTc3Ko2erZHq0qjghi4iIiESaKgs/M7sIuAuYCdwE/BHoCMw1s8vqPp7UpbxdB8hct5PxA9t4HUVERETqwdFm7rgd+Nw5d9rhBWb2N+A14AUzS3XO/aPu4kldmp6dD8A4zc0rIiISFY7W1dsLeKv8AufcTuAM4EXg/8zsD3WUTerY1Ow8+rdrSufURl5HERERkXpwtBa/ssoWOufKgKvMbCfwWzNrBnwT7HBSd9YW7CN7425+c1Yvr6OIiIhIPTla4bcSOAF4srKVzrk7zGwP8HtgXJCzSR2avsjXzTtW3bwiIiJR42hdvR8A55hZiyNt4Jy7D/gZ0CGYwaRuTc3KI6NTM9qlNPA6ioiIiNSTo7X4PQfsAFoC24+0kXPuETNbBwwMYjapIyu3FLJ8cyH3ju/jdRQRERGpR1UWfs65TcAT1TmQc+5d4N3aR5K6Ni0rjxiDswboNi4iIiLRpMZTtplZkpldbmatghlI6pZzjmnZ+RzftQUtk5O8jiMiIiL1qDZz9TbFN1NH3yBlkXqwJG8PuQX7GK8p2kRERKJObQo/AAtKCqk3U7PziIsxzujb2usoIiIiUs9qW/i5oKSQeuGcY1pWPiekp9KsUYLXcURERKSeqcUviizYsItNuw4wXvfuExERiUpHu51LVbYBXYDNQcoidWxqVh4JcTGc2lfjcURERKJRjVv8nHNlzrl1zrmiYAYCMLMLzWyJmZWZWUa55aea2TwzW+T/Orrcui/MbIWZLfQ/WlZx/I5mttfMfhHs7KGqtMwxPTufk3qk0SQp3us4IiIi4oHatPjVpcXAecAzFZYXAOOdc3lm1g/4CGhXbv2lzrnMahz/7/hmJYka363dwdbCIo3mFRERiWIhWfg555YBmFnF5QvKPV0CJJlZYiCtjmZ2LpAL7Kt90vAxNSuPBvGxjOl9xIZQERERiXC1HdzhpfOBBRWKvuf93bz3WMWqETCzRsCvgPvqK2QoKC4t44PFmxnTuyUNE0Ky1hcREZF64FkVYGYzgMpuJvdb59yUo+zbF3gIOK3c4kudc5vMLBl4C5gIvFhh1/uAvzvn9lZSF1b8GdcB1wF07Nixym1D3dyc7ezYd0jdvCIiIlHOs8LPOXdKTfYzs/bAO8Dlzrmccsfb5P9aaGavAMfxw8JvKHCBmf0FSAHKzOygc+7xSvJNAiYBZGRkhPX9Cqdl5ZGcGMeJPdK8jiIiIiIeCqt+PzNLAaYDdznn5pRbHgekOOcKzCweGAfMqLi/c25kuX3uBfZWVvRFkqKSUj5cspnT+rYmKT7W6zgiIiLioZC8xs/MJpjZRmAYMN3MPvKvuhnoDtxT4bYticBHZpYNLAQ2Ac/6j3W2md1f7ycRImatLKDwYAnjBrbxOoqIiIh4zJwL617MepGRkeEyM6tzl5jQc+urC5i1ahvf/fYU4mNDss4XERGRIDKzec65jMrWqRKIYAcOlTJj2RbO7NdGRZ+IiIio8Itkny3fyv5DpYxXN6+IiIigwi+iTc3KIy05kaFdWngdRUREREKACr8IVXiwmM9WbGVs/zbExlR9z0IRERGJDir8ItSMZVs4VFKmbl4RERH5ngq/CDU1K592KQ0Y1KGZ11FEREQkRKjwi0C79h9i1sptjB3Qhhh184qIiIifCr8I9OHizZSUOcYP0Ny8IiIi8l8q/CLQtOx8OrdoSL92TbyOIiIiIiFEhV+E2VZYxNycAsYPbIuZunlFRETkv1T4RZgPFudT5mCcunlFRESkAhV+EWZqVh49WjWmZ+tkr6OIiIhIiFHhF0Hydx/gu7U7NahDREREKqXCL4JMz84HYNxAFX4iIiLyQyr8IsjUrDz6tWtCl9RGXkcRERGREKTCL0Ks276PrI271c0rIiIiR6TCL0JM83fzjh2guXlFRESkcir8IsTUrDyO7ZhC+2YNvY4iIiIiIUqFXwRYtaWQ5ZsLGa9BHSIiIlIFFX4RYGp2PmYwtr+6eUVEROTIVPiFOecc07LzOL5LC1o2SfI6joiIiIQwFX5hbmn+HnK37WPcQLX2iYiISNVU+IW5qVn5xMYYZ/ZT4SciIiJVU+EXxg53857QPZXmjRK8jiMiIiIhToVfGFu4YRcbdx7QaF4RERGpFhV+YWxqVj4JsTGc1reV11FEREQkDKjwC1NlZY7pi/I4sWcaTZLivY4jIiIiYUCFX5j6bu0OtuwpUjeviIiIVJsKvzA1NTuPBvGxnNK7pddRREREJEyo8AtDJaVlvL9oM6N7t6RhQpzXcURERCRMqPALQ3NztrNj3yHGD1A3r4iIiFSfCr8wNC07j8aJcZzUM83rKCIiIhJGVPiFmaKSUj5cvJnT+rYiKT7W6zgiIiISRlT4hZnZKwvYc7BE3bwiIiISMBV+YWZqdh4pDeMZ0T3V6ygiIiISZlT4hZEDh0qZsXQLZ/ZrTUKcfnUiIiISGFUPYeTzFVvZd6iUcermFRERkRpQ4RdGpmblkdo4keO7tvA6ioiIiIQhFX5hYm9RCZ8t38rY/q2JjTGv44iIiEgYUuEXJmYs3UJRSZnm5hUREZEaU+EXJqZm5dGmaRLHdmzmdRQREREJUyr8wsCu/YeYtWob4wa0IUbdvCIiIlJDKvzCwEdLNlNc6tTNKyIiIrWiwi8MTMvOp1OLhvRv19TrKCIiIhLGVPiFuIK9RcxZXcC4AW0wUzeviIiI1JwKvxD3waJ8yhzq5hUREZFaU+EX4qZm5ZPesjE9WyV7HUVERETCnAq/EJa/+wDfrdvB+IFt1c0rIiIitabCL4RNz87HORg3oI3XUURERCQCqPALYVOz8+nbtgld0xp7HUVEREQigAq/ELV++36yNuzSoA4REREJGhV+IWraojwAxvZXN6+IiIgEhwq/EDU1K59BHVPo0Lyh11FEREQkQqjwC0Grt+5lWf4exg9QN6+IiIgEjwq/EDQtOw8zGKvRvCIiIhJEKvxCjHOOqVl5HNe5Oa2aJHkdR0RERCKICr8Qsyy/kJxt+zSaV0RERIJOhV+ImZqdR2yMcWa/1l5HERERkQijwi+EOOeYlp3HiO6ptGic6HUcERERiTAq/EJI1sbdbNhxQFO0iYiISJ1Q4RdCpmblER9rnN5X3bwiIiISfCr8QkRZmWN6dj4n9mhJ0wbxXscRERGRCKTCL0RkrtvJ5j0HGT9Q3bwiIiJSN1T4hYipWXkkxcdwSu9WXkcRERGRCKXCLwSUlJbx/qJ8xvRqRaPEOK/jiIiISIRS4RcCCvYeomtaI920WUREROpUSBZ+ZnahmS0xszIzyyi3/FQzm2dmi/xfR5db94WZrTCzhf5Hy0qO29nMDpTb5un6OqeqtG6axBs3DOcM3bRZRERE6lCo9isuBs4DnqmwvAAY75zLM7N+wEdAu3LrL3XOZR7l2DnOuWOCllREREQkTIRk4eecWwZgZhWXLyj3dAmQZGaJzrmieownIiIiEpZCsqu3ms4HFlQo+p73d+HeYxWrxv/qYmYLzGymmY2sh5wiIiIiIcGzFj8zmwFUdlHbb51zU46yb1/gIeC0cosvdc5tMrNk4C1gIvBihV3zgY7Oue1mNhh418z6Ouf2VPIzrgOuA+jYsWN1T0tEREQkZHlW+DnnTqnJfmbWHngHuNw5l1PueJv8XwvN7BXgOCoUfv7WwSL/9/PMLAfoAfzgukDn3CRgEkBGRoarSVYRERGRUBJWXb1mlgJMB+5yzs0ptzzOzFL938cD4/ANEKm4f5qZxfq/7wqkA7n1EF1ERETEcyFZ+JnZBDPbCAwDppvZR/5VNwPdgXsq3LYlEfjIzLKBhcAm4Fn/sc42s/v9+48Css0sC3gTuME5t6PeTkxERETEQ+acejGPJiMjw2VmHu0uMSIiIiLeM7N5zrmMytaFZIufiIiIiASfCj8RERGRKKHCT0RERCRKqPATERERiRIq/ERERESihAo/ERERkSihwk9EREQkSqjwExEREYkSKvxEREREooQKPxEREZEooSnbqsHMtgHrqtgkFSiopzihRucevaL5/HXu0Smazx2i+/zD7dw7OefSKluhwi8IzCzzSHPiRTqde3SeO0T3+evcde7RKJrPP5LOXV29IiIiIlFChZ+IiIhIlFDhFxyTvA7gIZ179Irm89e5R6doPneI7vOPmHPXNX4iIiIiUUItfiIiIiJRQoVfFczsDDNbYWarzezXlaw/ycx2m9lC/+N35dalmNmbZrbczJaZ2bD6TV97tTz/n5nZEjNbbGavmllS/aavnaOdu3+bk/znvcTMZgaybyir6bmbWQcz+9z/fl9iZrfVb/Laq83v3b8u1swWmNm0+kkcXLV834f1Z14tzz2iP+/M7JflPucXm1mpmTWvzr6hrqbnHtafd845PSp5ALFADtAVSACygD4VtjkJmHaE/f8FXOP/PgFI8fqc6uv8gXbAGqCB//nrwJVen1OQzz0FWAp09D9vWd19Q/lRy3NvAxzr/z4ZWBkt515u/c+BV470uRDKj9qefzh/5tXyfR/xn3cVth8PfFaTfUPtUctzD9vPO7X4HdlxwGrnXK5z7hDwH+Cc6uxoZk2AUcA/AZxzh5xzu+oqaB2p8fn7xQENzCwOaAjk1UHGulKdc78EeNs5tx7AObc1gH1DWY3P3TmX75yb7/++EFiG749iuKjN7x0zaw+MBSbXU95gq/H5R8BnXq1+90T+5115PwZereG+oabG5x7On3cq/I6sHbCh3PONVP5LHWZmWWb2gZn19S/rCmwDnvd3+0w2s0Z1nDfYanz+zrlNwMPAeiAf2O2c+7iuAwdRdc69B9DMzL4ws3lmdnkA+4ay2pz798ysMzAI+KaugtaB2p77P4A7gbI6TVl3anP+4f6ZV+Nzj5LPOwDMrCFwBvBWoPuGqNqce/l1nQmjzzsVfkdmlSyrOAR6Pr5pUQYCjwHv+pfHAccCTznnBgH7gHC79qHG529mzfD9r6kL0BZoZGaX1V3UoKvOuccBg/G18JwO3GNmPaq5byirzbn7DmDWGN+H4+3OuT11FbQO1PjczWwcsNU5N6+OM9al2vzuw/0zrza/+2j4vDtsPDDHObejBvuGotqcu+8AYfh5p8LvyDYCHco9b0+F5nvn3B7n3F7/9+8D8WaW6t93o3PucPX/Jr4PxXBSm/M/BVjjnNvmnCsG3gaG10/soDjqufu3+dA5t885VwDMAgZWc99QVptzx8zi8X0I/ts593Y95A2m2pz7COBsM1uLr7totJm9XPeRg6q27/tw/syrzblHw+fdYRfz327eQPcNRbU59/D9vPP6IsNQfeD7310uvv/FHb7os2+FbVrz33shHoevqf/w89lAT//39wJ/9fqc6uv8gaHAEnzXuhi+i75v8fqcgnzuvYFP/ds2BBYD/aqzbyg/annuBrwI/MPr86jvc6+wzUmE5+COWp1/OH/m1fJ9H/Gfd/7tmgI7gEaB7huqj1qee9h+3sUhlXLOlZjZzcBH+Eb+POecW2JmN/jXPw1cAPzUzEqAA8DFzv+OAG4B/m1mCfjeWD+p95OohVqe/zdm9ia+ruASYAFhdNfz6py7c26ZmX0IZOO7pmuyc24xQGX7enIiNVCbczezE4CJwCIzW+g/5G+crzU45NX29x7ugnD+YfuZF4R/8xH9eeffdALwsXNu39H2rd8zqLnanDu+Vv6w/LzTzB0iIiIiUULX+ImIiIhECRV+IiIiIlFChZ+IiIhIlFDhJyIiIhIlVPiJiIiIRAkVfiJSK/4prL7wOkd9M7OLzGyJmRWZmTOzFK8zVcXMOvtzXul1lkCY2ZX+3J29ziISCVT4iUQQM3vbzIrNLK2KbW7z/yEdX5/ZIomZdQf+DWwGforvfl77qtxJRCQE6AbOIpHlJXw3G70Y3/zJlbkMKAA+DNLPPC1IxwknJ+H7/LzDObfQ2ygiItWnFj+RyDId39RClU4Sb2Y9gAzgP843r2iNmVlDAOfcIefcodocKwy19H/dVZOdzayBmenzV0TqnT54RCKIvwB7HTjOzNIr2WSi/+tLZpZgZveZ2bdmtsPMDpjZwsquATOztWY2w8xGmdlcMzsA/NG/7gfX+JnZHWY228y2+a+BW25mvzAzq7DdF2a22sy6m9lHZrbPzLaa2Z8rK4zM7Hwz+9LMCs1sj5llmtnVFbY51szeM7Od/nPKNLNzq/P6mVmMmd1pZiv8ufPM7Iny1++Z2VrgQf/TNf5u8xeqOObha9ROMbO/mVkevm7hJmbW3Mz+YmZZ/vPZ5399x1VyHGdmk83sTDNbYGYH/a/dJZVs28bM3vC/TjvM7J9AkyPkG2Fmn/q33ev/ftgRzmG0mT1kZpv9279hZilmFmdmf/S/Xgf8r/8RLzcod9xG/uPl+M9nu5l9bWYXVGPf8f5t95vZLjObYma9K2xzrz93PzP7l/89scfMXjWzlpUcs8bvHZFwocJPJPK85P96aSXrLgFWOue+xVcI3AB8DdwD3IWvBet5M7u2kn27AFOAr4Bbgc+ryPBzYBnwgP/75cBf/c8rSgZm4Jvf9Q5gDvAr4H8ymNmvgTeBhviKzl8D84Dx5bYZ6d+/Hb7i7JfAfuAdM/txFXkPexJ4CFjpz/0OvtfoU/PNQQtwuz8HwM/wFdPPVOPY/wCG+Y//W+AQ0BW4CPgYuBP4PdAAeM/MKutCHwK8ALwH/ALYi6+I/77gMbMk4FPgHGCy/5hd8U0o/z/MbBTwmX/9H/G9Zt2Az81sRCU//2HgOHy/x38B5wP/BJ4CjvcfYzIwDnikGq/Jk/hew6nAzf6fvxIYWtVO/t/lFHyv1d34XtsTgLnmu/6yopeB9vje58/jm2f843K/02C8d0TCg3NODz30iLAHsBpYVWHZCMABd/ufxwKJlez7SSX7rvXve14l238BfFFhWcNKtpuMr1BJrLCvA26osO1C4Ltyz7sAJfgKpPgK2x6ec9zwFZuzgNjy64EvgQ2Htz3Ca9bPn+XVCstv8i+/sdyyu/3LOlfjd3Glf9t5lWRPLJ/VvywBWAJ8UmG5878G/cotawUUAX8tt+wW/7Y/Kbcs1v+6OODKcsszgZ1Aq3LL2gC7gW8rOYdZQEy55a8BZfj+E1BxeTHQ+CivzU7giWq+fp39z+OBfGBV+eMDA4BS4PVyy+717/tZhffEtf7l1wfjvaOHHuH0UIufSGR6GehuZseXW3YZvj92/wZwzpU654oAzCze3+2Yiu+PZHcza1rhmJvxtYAdlXNuv/+4cWbWzH/cL4BGQM8KmxfjKwrLm4mvFeqw8/AVL/e6CtcmOuec/9uBQC//+TUzs1T/z20BvI+vxadHFbEPd6/+tcLyZ/G1hP6g+zVAz1aSvcg5VwpgZolm1gJfS+wsYHAlx5jlnFtcbv8t+FpTy79W44Dt/LflF//P+J/BPmbW2v8zXvIf5/C2+fjeP0PMrFUl51BW7vlX+Iqj5ypZHgd0rOQcytsFDDWzDkfZrrzBQGvgKefc3nK5s/ENWDrTfniZwGOHX2e/F/jf32lt3zsiYUOFn0hketn/9TIAf5fWRcCXzrk1hzcysyvMLBs4iK9Y2Ib/2j2gYuGXW67IqpKZnWVmXwMH8A022cZ/C5GUCptvcs6VVFi2E2he7vnh7rtFVfzYw3+Yn/b/vPKPw9fk/eC6rnI6+78uL7/Q+a6bXI2v1bE2ciouMJ87zGwlvt9BgT/vDfzwdQJYV8myiq9VJ3y/q4qv6YoKzzv7vy7nh5ZW2Oaw9RWe7zrK8maVHLu8O4DewDrzXV/6VzOrrOAt73CmI+VuDFS8vvB/zt1fgK8pd6zavndEwoZu5yISgZxzq83sK+BHZnY7cBa+4uD7ViAz+xG+lo/pwN+ALfha387Cd91Vxf8YHqjOzzaz4fiu2foKuBHYhO96tmPxXd9W8bilHJ0dfZPvj/sb4LsjbLP4CMur8/OrVfRWobLX707gz/h+L/fhK/xKgZ/gux6zoiO9Vlbh+8qyVuc1rLhtxeMc6edXJ9cPOOfeNrMv8V2neQpwFXCHmf3WOfen6oat5OdVzH2016Mu3zsiIUWFn0jkegnfxfOn42v5KwLeKLf+YnytHuPLt+SZ2eha/twL8RV6pzjnDpY7btcj73JUq/xf+wNzj7DNav/Xfc65GTX4GWv9X3sB8w8vNLN4fF2pX9fgmEdzMb7rIy8vv9DMrqrFMdcCg80srkKrX8WuyrX+r70qOcbhZZW1MAaVc24rvgEi/zTfLYKmA/eZ2cMVu8b91pbL+H6Fdb3wXUdaUMnyZYef+H+nnfFdvwe1f++IhA119YpErtfwFWA34buWaapzble59Yevyfr+c8B/jVltio7Dxy3Dd03e4eMm4Rt0UFNv42tVus//R/t7Zt/fImY+vgLxDqtk+rRq3F5kmv/rzyssvwZfl+XUADNXRxkVPofNdxueCbU45nR816YdvnUPZhZLhdffObcZ3+COieVvbeK/9m8ivsEdW6gjZhZb8TpS/7WhK/AN4Gh0hF0z8V1veoOZfb+NmfUDzgDer3C9IcAt/tfgsCvxdaVP9z+v7XtHJGyoxU8kQjnndpjZ+8C5/kUvVdhkCr5BE9PM7F181zBdB+ThGy1aU+/h6yqeYWYv4btdyxX4rmGrEefcGjP7Pb7biHxjZq/jG3naD98o1POcc2Vm9hN8I3+Xmtlz+FqHWuO7PUgffLcqOdLPWGxmzwDXm1kTfAMF+uC73m4+vlapYJuCr5j9N77BL53wdY8vB46p4TGf9R/jGTPrj+82OedT+X387sA3ivtr/7kbcD2QxA8L4GBLBjaZ2TtAFr5rQQfhK7Q/qPCflO8550rM7Of4BmLMMbN/4Tu3W4BCfLfKqag5vtu3vIPvetGb8HXdPu8/Zq3eOyLhRIWfSGR7CV/htx34oPwK59yL/ha+m/Ddc209vvu07cb/B7EmnHMzzWwivuul/gZsxXct4Wx8f1hretwHzSwH3330fofv1iYr8HVnH95mjpkdh+9+bdfha9XZgq+wqKwgqOhGfIXSNfhaj7YDk4DfurqZneRP+G7pMhFfcbYS3/3selDDws85d8DMxuD7nV6Hr9X3Hf/zrArbzvJvez++1wzgW+BS59yRutSDZT/wOL5r+8biex3W4xtc9JeqdnTOvWpm+/D9Tv+I7xy/AO5yzq2uZJfL8F1P+Qd8f/feBm47PKrdf8zavndEwsLh+1+JiIhEFDO7F98NrDs45zZ6HEckJOgaPxEREZEoocJPREREJEqo8BMRERGJErrGT0RERCRKqMVPREREJEqo8BMRERGJEir8RERERKKECj8RERGRKKHCT0RERCRKqPATERERiRL/Dzvldao4HXGyAAAAAElFTkSuQmCC\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": {
    "execution": {
     "iopub.execute_input": "2021-02-02T07:04:51.959585Z",
     "iopub.status.busy": "2021-02-02T07:04:51.958856Z",
     "iopub.status.idle": "2021-02-02T07:05:14.522759Z",
     "shell.execute_reply": "2021-02-02T07:05:14.523634Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnQAAAHoCAYAAADAJXJqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABj/ElEQVR4nO3dd3gc5dXG4d9Rd5eLXOXeG+5gm96rKQaCE3oJJZDQ8kEICR0CSSBACC1ACJAAoYViqmkGZGNccS9yt3GR3Kva+f7YFRFCkrXSSrMrPfd17WXvzOzqGUm2jt553zPm7oiIiIhI/EoIOoCIiIiIVI8KOhEREZE4p4JOREREJM6poBMRERGJcyroREREROKcCjoRERGROJcUdICgtWrVyrt06RJ0DBEREZF9mjZtWo67Z5TeXu8Lui5dujB16tSgY4iIiIjsk5mtKGu7LrmKiIiIxDkVdCIiIiJxTgWdiIiISJxTQSciIiIS51TQiYiIiMQ5FXQiIiIicU4FnYiIiEicU0EnIiIiEudU0ImIiIjEORV0IiIiInFOBZ2IiIhInFNBJyIiIhLnVNCJiIiIxDkVdCIiIiJxTgWdiIiISJxTQSciIiIS51TQiUhMcXcKCouCjiEiEleSgg5QFjM7E7gN6Avs7+5Tw9uPBu4FUoA84P/c/ZPwvs+AdsDu8Nsc4+4baje5iFTFtj35ZC3JZeLijUxctJH12/ZwZJ82nD4sk8N6Z5CcqN89RUQqEpMFHTAHGAs8UWp7DjDG3dea2QDgA6BDif1nFxd/IhK7CoucOWu2MnHRRiYu3sj0lVsoLHIapSQyqnsrDuudwXuz1/H+3HW0aJTCyYPac/rQTAZ0aIqZBR1fRCTmxGRB5+7zgR/9x+3uM0o8nQukmVmqu++txXgiUgXrt+3h80Ub+WJxDl8u3sjmXfkADOjQlMsO6cYhvTIY2qk5KUmh0bhbx/Rn4qKNvD59Df/+eiXPZi2nV5vGnD40k1OHdKBN07QgT0dEJKaYuwedoVzhy6i/LmvUzczOAC5396NKHNsSKAReA+7yck7OzC4FLgXo1KnTsBUrVtRIfpH6bE9+Id8s3xQahVuUw8L12wFo1TiVQ3q14pCeGRzUsxWtGqfu87227srn7W/X8tr01cxYuYUEg4N6ZnD60A4c068tDVISa/p0RERigplNc/fhP9oeVEFnZhOAtmXsutnd3wwf8xllFHRm1h94i9A8uezwtg7uvsbMmhAq6F5w9+f2lWP48OE+daqu0opUl7uTvXEHny/KYeKijXy9LJc9+UWkJCYwvEtzDumVwSE9M+jbrkm1Lpsu3biD16ev4Y0Za1izZTeNU5M4cWA7xg7twP5dW+iSrIjUaTFX0FVGWQWdmWUCnwAXuvtX5bzuAmC4u1+1r4+hgk6k6rbuyufLJaEC7ovFG1m7dQ8A3Vo1ChVwvVoxsltLGqZEf3ZHUZEzeVkur09fw7uzv2NXXiEdWzRg7JBMxg7tQOeWjaL+MUVEglYnCjozSwc+B+5w99dKHJcEpLt7jpklAy8CE9z98X19DBV0IpVXUFjErNX/W8wwa9UWihyapCZxYI9WHNIrg4N7tqJji4a1mmtXXgHvz1nH69PX8FV2Du4woktzxg7N5MT92tE0LblW84iI1JS4KujM7DTgr0AGsAWY6e7HmtnvgJuAxSUOPwbYCUwEkoFEYAJwnbsX7utjqaATqdjaLbu/L+C+XJzDtj0FmMF+mekc2jNUxA3umE5SjLQWWbtlN2/MWMNr01ezdONOUpMSOKZ/W8YO7cDBPVrFTE4RkaqIq4KuNqmgE/mh3XmFTF6WG17MsJHsjTsBaNM0lUN6ZnBIrwwO6tGK5o1SAk5aMXdn1uqtvD59NW/NWsuWXflkNEnltCEdGDu0A33aNg06oohIxFTQlUMFndR37s7C9du/X406Zfkm8gqKSE1KYP+uLTi0V6iI69m6cdwuONhbUMinCzbw2vQ1fLpgAwVFTv/2TTl9aCYnD25fqZW2IiKxQAVdOVTQSX20dXc+n4dH4L5YvJH120KtHHu2bhxezJDBAV1bkJZc99qB5O7Yy9uz1vLa9DXMXrOVpATjsN4ZjB2ayZF9W5OaVPfOWUTqDhV05VBBJ/XN/O+2cf4zU9iwfS/NGiRzUI9WHNKrFQf3zKB9eoOg49WqReu389r01fx3xhrWbwt9PsYMasfYoZkM6ZgetyOSIlJ3qaArhwo6qU++XprLJc9NpVFKEg+OG8yILi1ITFDRUljkfLkkh9enr+aDuevYk19Et1aNGDu0A6cNzaRDPSt0RSR2qaArhwo6qS8+nLuOq16cQWbzBjx/8QEqUsqxfU8+781ex6vTVzNl2SbMYFS3lowdmsnxA9rSKDUm75goIvWECrpyqKCT+uA/36ziN69/y8DMdP5xwQhaxPgK1VixMncXb8xYw+szVrMidxcNUxI5bkBbTh+ayahuLUnQ6KaI1DIVdOVQQSd1mbvz+OdLue/9BRzcsxWPnzNMI0xV4O5MW7GZ16av5p1Z37F9bwHtm6Vx2tAOjB2aSfeMxkFHFJF6QgVdOVTQSV1VVOTc/e58nv5yGScPas+fzxxESpKa6lbXnvxCPpq3ntemr2bioo0UOQzumM7pQzswZlB70htq9FNEao4KunKooJO6KL+wiBte/ZY3ZqzhgtFduOWkfro8WAM2bNvDmzPX8tr01SxYt52UxASO6NOa04dlcljvDJJ1VwoRiTIVdOVQQSd1za68An7xr+l8tnAjvz6mF1ce3kPtN2qYuzPvu228Nm0Nb85cQ+7OPFo2SuHkwe05fWgm/ds31ddARKJCBV05VNBJXbJlVx4XPvsNs1Zt4e7TBvLT/TsFHaneyS8sYuKijbw2fTUT5m0gr7CI3m2ahFqgDOlA66ZpQUcUkTimgq4cKuikrvhu627Oe3oKKzbt4uFxgzluQLugI9V7W3fl8/a3oUuyM1ZuIcHg4J4ZjB3agWP7t62Td+IQkZqlgq4cKuikLliyYQfnPf012/YU8PfzhjOqe8ugI0kpSzfu4PXpa3hjxhrWbNlNk9QkThjYjtOHZTKiS3NdkhWRSlFBVw4VdBLvZqzczEXPfkNiQgLPXjiCAR2aBR1JKlBU5Exelstr09bw3pzv2JVXSMcWDRg7JJMzhmXSsUXDoCOKSAxTQVcOFXQSzz5ftJHLn59GRpNUnr94fzq3bBR0JInArrwC3p+zjtenr+Gr7BySExJ49OyhHNWvTdDRRCRGlVfQaU29SJx6c+YaLn72G7q0asSrV4xSMReHGqYkMXZoJi9ccgBf3ngEfdo14fIXpvHe7O+CjiYicUYFnUgc+sdXy7j6pZkM7dycly8bSesmWjkZ7zqkN+CFSw5gv8xmXPXiDN6cuSboSCISR1TQicQRd+fPHyzk9rfncUy/Njx30f40TUsOOpZESdO0ZJ67+ACGd27OtS/P5NVpq4OOJCJxQgWdSJwoLHJ++8ZsHvl0CeNGdOTRs4eq7UUd1Dg1iWcv3J/R3Vvxf6/O4t9frww6kojEARV0InFgT34hv/jXNF6csoqrDu/BH8YOJEm3laqzGqQk8tT5wzmsVwa/fWM2/8xaHnQkEYlx+okgEuO27cnn/Gem8MHc9dw6ph+/Pra3epbVA2nJiTx+7jCO6deGW9+ay98nLg06kojEMBV0IjFsw/Y9jHtiMtNWbOahcYO58MCuQUeSWpSalMjfzh7Kifu14+535/PIJ4uDjiQiMSop6AAiUrYVuTs59+kpbNy+N3T5rXfroCNJAJITE3jorMGkJibw5w8XsbegiOuO7qVRWhH5ARV0IjFo7tqtnP/MNxQUFfHvnx/AkE7Ng44kAUpKTOBPZw4iOTGBv36yhLyCIn5zfB8VdSLyPRV0IjFm8tJcfv7PqTRJS+KlS0fRo3WToCNJDEhMMP4wdiApSQk8MXEpewuKuHVMPxV1IgKooBOJKe/PWcevXppBpxYNee6i/Wmf3iDoSBJDEhKMO07pT0pSAk9/uYy8wiLuOmUACQkq6kTqOxV0IjHipSkr+e0bsxnUMZ1nzh9B80YpQUeSGGRm/O7EvqQkJfDYZ9nkFRRx3+n7kaiiTqReU0EnEjB359HPsvnTBws5tFcGj50zlIYp+qcp5TMzbji2N6lJCTw4YTH5hUXcf+Yg9SYUqcf0U0MkQEVFzh3vzOPZrOWcOrj99xPfRfbFzLjmqF6kJCXwx/cXkldQxEPjhpCSpO8fkfpIBZ1IQPIKivj1K7N4a9ZaLj6oKzef0FdzoSRivzisBymJCdw1fj75/5rG384eSmqSbgknUt/oVzmRAOzcW8Alz03lrVlrufG4PvzuRBVzUnWXHNyNO0/pz4T5G/j5c9PYk18YdCQRqWUq6ERq2aadefzsqa/5cvFG7jt9IFcc1l2tJ6Tazh3VhftOH8gXizdy0bPfsCuvIOhIIlKLVNCJ1KI1W3Zz5uNZzP9uG4+fM4yzRnQKOpLUIWeN6MQDPxnE5KW5nP/MFLbvyQ86kojUEhV0IrVk8frtnPFYFhu27+X5i/bnmP5tg44kddBpQzJ5+KdDmL5yC+c+PYWtu1XUidQHKuhEasG0FZs584lJFBQ5L186igO6tQw6ktRhJ+3XnkfPHsrctVs5+6nJbN6ZF3QkEalhKuhEatinCzdwzlNfk94gmdcuH02/9k2DjiT1wLH92/LkucNZtH4HP/37ZHJ27A06kojUIBV0IjXovzPW8PN/TqVbRiNeuXw0nVo2DDqS1COH92nN0+cPZ3nuTsY9OZkN2/YEHUlEaogKOpEa8vSXy7jm5ZmM6NKCly4dSUaT1KAjST10cM8Mnr1wf9Zu2c1ZT07mu627g44kIjVABZ1IDXhwwiLufGcexw9oyz8uHEGTtOSgI0k9NrJbS56/eH9ytu/lJ09MYtWmXUFHEpEoU0EnEmWzVm3hoY8XM3ZIBx752VDSktW1X4I3rHMLXrjkALbuyuesJyaxPGdn0JFEJIpU0IlEUVGRc+tbc2nVOJXbT+lPou7+IDFkUMd0Xrx0JLvzC/nJE5NYsmFH0JFEJEpU0IlE0WvTVzNz1RZ+c1wfXWaVmNS/fTNeunQURQ7jnpzEwnXbg44kIlGggk4kSrbtyee+9xcwpFM6pw3pEHQckXL1btuEly8bSWKCMe7JScxZszXoSCJSTSroRKLkoQmLyd2Zxx0nDyBBl1olxnXPaMx/LhtFw5Qkfvb3ycxctSXoSCJSDSroRKJg8frt/DNrOeNGdGRgZrOg44hUSueWjXj5spGkN0zhnKe+ZuryTUFHEpEqUkEnUk3uzm1vz6VhSiK/PqZ30HFEIpLZvCEvXzaS1k1SOe+ZKUzKzg06kohUgQo6kWr6YO46vlqSy/XH9KZlYzUPlvjTrlkDXrp0JB3SG3DBP6YwcdHGoCOJSIRU0IlUw+68Qu58Zz592jbh7AM6BR1HpMpaN03jpUtH0i2jMZf8cyqfLFgfdCQRiYAKOpFqePzzbNZs2c1tJ/cnKVH/nCS+tWycyos/P4DebZtw2fPTeH/OuqAjiUgl6SeQSBWt2rSLxz/P5qT92jGyW8ug44hERXrDFP718wMY0KEZV/57Om/PWht0JBGpBBV0IlV01/h5JJhx84l9g44iElVN05J5/uIDGNapOVe/NIPXpq0OOpKI7IMKOpEq+GLxRj6Yu56rjuhBu2YNgo4jEnWNU5N49qIRjOrekl+/OouXpqwMOpKIVEAFnUiE8guLuP3teXRu2ZBLDu4adByRGtMwJYmnzx/BIT0z+M3rs3lu0vKgI4lIOVTQiUTon1nLWbJhB7ec1I/UpMSg44jUqLTkRJ48bxhH9W3DLW/O5akvlgYdSUTKoIJOJAIbtu/hwQmLOax3Bkf0aR10HJFakZqUyGPnDOXEge24a/x8/vbpkqAjiUgpSUEHEIkn9723kL0FhdxyUj/MdL9WqT+SExN4aNxgkhONP32wkL0FRVx7VE/9OxCJESroRCpp2orNvDZ9NZcf2p1uGY2DjiNS65ISE7j/J4NJTkzg4Y8Xk1dQxI3H9VZRJxIDVNCJVEJRkXPbW3Np0zSVXx7RI+g4IoFJTDDuO30/UpISePzzbI1Yi8QIFXQilfCfqauYvWYrD40bTKNU/bOR+i0hwbjr1AGkJCXwj6+Wk1dQxJ2nDCAhQUWdSFD0k0lkH7buyuePHyxkRJfmnDyofdBxRGKCmX2/0vvxz7PJLyziD2P3I1FFnUggVNCJ7MNfJixiy648bjt5f11WEinBzLjxuN6kJP1vTt2fzxyk+xqLBEAFnUgFFqzbxvOTV/CzAzrRv32zoOOIxBwz47qje5GalMCfPlhIXmERD40bQrKKOpFapYJOpBzuzq1vzqVJWhLXH9076DgiMe3Kw3uQmpTAXePnk1cwnb+dPUSNt0VqkX6FEinH+Nnf8fWyTfz6mN40b5QSdByRmHfJwd2445T+TJi/nsuen8ae/MKgI4nUGyroRMqwK6+Au8fPp3/7pvx0/05BxxGJG+eN6sIfxg7k80Ubufif37ArryDoSCL1ggo6kTI8+mk2323dw+0n99eqPZEI/XT/Tvz5jEFMys7lgme+YcdeFXUiNU0FnUgpK3J38uTEpZw6uD3Du7QIOo5IXDp9WCYPjRvCtJWbOffpr9m6Oz/oSCJ1mgo6kVLufGceyYnGTSf0DTqKSFwbM6g9f/vZUOas2co5T33Nll15QUcSqbNU0ImU8OnCDUyYv4FfHtmTNk3Tgo4jEveOG9CWJ84dxsL12xn35GRyd+wNOpJInVRuQWdmRWZWGOkjGqHM7EwzmxvOMLzE9qPNbJqZzQ7/eUSJfSlm9qSZLTKzBWZ2ejSySP2RV1DEnW/Po1urRlx0YNeg44jUGUf0acNT5w1nWc5Oxj05mQ3b9gQdSaTOqagP3R2Al9p2KjAA+ABYCBjQGzgGmA28GaVcc4CxwBOltucAY9x9rZkV5+gQ3nczsMHde5lZAqDJTxKRZ75axtKcnTx74QhSkjR4LRJNh/TK4NkL9+fif37DuCcn86+fH0C7Zg2CjiVSZ5Rb0Ln7bSWfm9lFQFtggLsvLLWvL/ApsDIaodx9fvh9S2+fUeLpXCDNzFLdfS9wEdAnfFwRoeJPpFLWb9vDXz9ezFF923BY79ZBxxGpk0Z1b8lzF+3PBf/4hp88MYl/XzKSji0aBh1LpE6IZBjiBuCR0sUcfF+A/Q24MVrBKuF0YIa77zWz9PC2O81supm9YmZtajGLxLk/vDuf/CLn9ydpIYRITRrepQUvXHIAW3flM+7JyazI3Rl0JJE6IZKCrjOwu4L9u8LHVIqZTTCzOWU8TqnEa/sD9wGXhTclAZnAV+4+FJgE/LmC119qZlPNbOrGjRsrG1nqqG+Wb+K/M9dy6cHd6NyyUdBxROq8wR3T+ffPR7Irr4CzntBCCZFoiKSgWwRcWmI07Htm1hy4lNC8ukpx96PcfUAZjwrn4ZlZJvAGcJ67Z4c35xIqKN8IP38FGFrBx37S3Ye7+/CMjIzKRpY6qLAodL/W9s3S+MXh3YOOI1JvDOjQjOcvPoBNO/O48bXZuJeesi0ikYikoPst0B1YbGYPmNnlZnaZmf2FULHXjdDChBoTLibHAze5+1fF2z30P8HbwGHhTUcC82oyi9QNL05ZybzvtvHbE/vSMKWiNUIiEm0DOjTjxuP7MGH+ev49JSpTsEXqrUoXdO4+HjiW0MKHa4BHgceAq8Pbjg8fU21mdpqZrQZGAePN7IPwrquAHsDvzWxm+FE8g/1G4DYz+xY4F7g+Glmk7tq8M48/f7iQUd1acuLAdkHHEamXLhzdhYN7tuLOd+axZMOOoOOIxC2ryjB3eMFBF0JtS5a5+/oo56o1w4cP96lTpwYdQwLwu//O5sUpqxj/q4Po07Zp0HFE6q0N2/Zw7IMTaZ/egDd+caDaBolUwMymufvw0tur9K/G3de7+9fuPjmeizmpv+au3cq/v17JuSM7q5gTCVjrpmncd/p+zF27jfs/qvRUbBEpIaKCzszSzeze8GrUnWa2I/z3e8paLCESi9yd296aS3rDFK49qlfQcUQEOKZ/W352QCeenLiUrCVqIyoSqUoXdGbWAZhBqB8dhBYnvEfobhK/AaabWfuoJxSJsrdmreWb5Zu54djeNGuYHHQcEQn73Yl96dqqEdf9ZxZbduUFHUckrkQyQvcHoA1wUri9yE/c/Ux3HwicGN73h5oIKRItO/cWcM+789kvsxk/Gd4x6DgiUkLDlCQeHjeE3J17uel1tTIRiUQkBd1xwEPu/m7pHe7+HvBX4PhoBROpCX/9ZAnrt+3ltpP7k5Bg+36BiNSqAR2acf0xvXlvzjpembY66DgicSOSgq4JsKaC/avDx4jEpKUbd/D0l0s5fWgmQzs1DzqOiJTj0oO7MapbS257ay7Lc3RrMJHKiKSgWwicYWY/eo2ZJQJnEMGdIkRqk7tzxzvzSE1K5MbjewcdR0QqkJBg3P+TQSQnJnD1yzPJLywKOpJIzIukoHsYOAT4xMxOMbM+4cepwATgYODB6EcUqb6P52/gs4UbueaonrRukhZ0HBHZh/bpDbjntIHMWrWFhz9eHHQckZhX6Xsdufsz4bsy3Aq8XmKXAXuB37r7s9GNJ1J9e/ILueOdefRo3ZjzR3cJOo6IVNKJ+7Xjs4WZ/O3TJRzcM4P9u7YIOpJIzIqoD5273wtkAmcTurfrb4GfAh3c/b7oxxOpvqe/XMbKTbu4bUx/khPVgV4kntx6cn86tmjItS/PZOvu/KDjiMSsiH+6uXuuu7/k7veFHy+7+6aaCCdSXWu37OaRT5ZwXP+2HNSzVdBxRCRCjVOTePCswazbtodb3pwTdByRmBVxQWdmx5nZI2Y23szeCf/9mJoIJ1Jd97w7nyJ3bj6xb9BRRKSKhnRqzjVH9uTNmWv574yKmi2I1F+R3CkixczeJHSHiF8AI4ADwn9/z8z+a2YpNRNTJHKTsnN559vvuPzQ7nRs0TDoOCJSDb84vAcjujTn9/+dw6pNu4KOIxJzIhmhuxUYA9wPtHb31u6eAWQAfwZOBn4f/YgikSsoLOL2t+fSIb0BVxzWPeg4IlJNiQnGAz8ZDMC1L8+kQK1MRH4gkoLuZ8AL7n6Du39/5+TwnLobgReAc6IdUKQq/vX1Shas287vT+pLWnJi0HFEJAo6tmjIXacNYOqKzTz2WXbQcURiSiQFXXsgq4L9k4B21YsjUn25O/Zy/4cLOahHK47t3zboOCISRacM7sCpg9vz4MeLmbFyc9BxRGJGJAXdWmBkBfsPAL6rXhyR6vvzhwvZlVfIbSf3w0z3axWpa+44dQBtm6Zxzcsz2bG3IOg4IjEhkoLuReBcM7vLzL6/EaaZNTezO4FzgX9HO6BIJL5dvYWXvlnF+aO70KO1bi0sUhc1TUvmwXGDWbVpF7e/NTfoOCIxIZKC7nbgA0LNhHPMbJ2ZrQNygJvD++6IfkSRyikqcm59ay4tG6Vy9VE9g44jIjVoRJcWXHV4D16Ztprx3+rikEgkt/7aC5xgZicBJwJdwruWA2+7+7tRTycSgddnrGHGyi386Yz9aJqWHHQcEalhvzyyJxMX53DT698ypFM67dMbBB1JJDBVuVPEO+5+hbsfH35coWJOgrZ9Tz73vreAIZ3SOX1oZtBxRKQWJCcm8NC4wRQWOdf9ZyaFRR50JJHA6MaWUic8/PFicnfu5faT+5OQoIUQIvVF55aNuO3k/kxeuom/f7E06Dgigan0JVcAMzsauAToBrQASv/kdHdXF1epVUs2bOcfXy3nrOEd2S8zPeg4IlLLzhiWyWcLN3L/hws5sHsrBmY2CzqSSK2L5NZf1wLvA4cCa4CJwOelHhNrIKNIudyd29+eR4OURH59bO+g44hIAMyMu08bQKvGqVz98gx25amVidQ/kYzQXUuoaDvO3fNqKI9IRD6Yu54vFudw65h+tGqcGnQcEQlIesMU7v/JIM5+6mvuGj+fe04bGHQkkVoVyRy6VsDLKuYkVuzJL+Su8fPo3aYJ547sHHQcEQnY6O6tuOyQ7vz765V8OHdd0HFEalUkBd10QnPnRGLCE58vZfXm3dx2cn+SErW+R0TguqN7MaBDU2587Vs2bNsTdByRWhPJT8FrCN0p4ugayiJSaas37+LRz5Zw4n7tGNW9ZdBxRCRGpCQl8NC4IezOL+T6V2ZRpFYmUk+UO4fOzD4sY/N24H0zW06ooXBhqf3u7sdGLZ1IOe4ePx8zuPmEvkFHEZEY0z2jMb8/qR83vzGHf2Qt5+KDugYdSaTGVbQoohdQ1q82KwmN7OnyqwTiqyU5vDdnHdcf3Uud4UWkTD/bvxOfLtjIfe8tYHT3lvRt1zToSCI1qtyCzt271GIOkUr7y0eLyGzegJ8fot8pRKRsZsZ9pw/kuIe+4OqXZvDWVQeRlpwYdCyRGqOZ5BJXluXsZOqKzZwzsrP+cxaRCrVsnMqfzxzEovU7uPe9BUHHEalRKugkrrw2bTUJBqcN6RB0FBGJA4f2yuCiA7vybNZyPl2wIeg4IjWm3ILOzIrMrMDMUko8L9zHQ+25pcYUFjmvTV/NIb0yaNM0Leg4IhInbjiuN33aNuH/Xp1Fzo69QccRqREVLYq4g9CiiIJSz0UCMSk7l++27uHmE7WyVUQqLy05kYfGDWHMI19yw6vf8vT5wzErfStykfhW0aKI2yp6LlLbXp22iqZpSRzVt03QUUQkzvRu24TfHt+H296exwuTV3DuqC5BRxKJKs2hk7iwbU8+789dx8mD22sxhIhUyfmju3BorwzuGj+fxeu3Bx1HJKoqaix8SFXe0N0nVj2OSNne/fY79uQXccawjkFHEZE4ZWb86cz9OP7BL/jVSzP575WjSU3SL4hSN1Q0h+4zIpszZ+Hj9a9Dou7Vaavp0boxgzKbBR1FROJY6yZp/PGM/bj4n1P58wcLufnEfkFHEomKigq6w2sthUgFinvP/eb4PprILCLVdmTfNpw7sjN//2IZh/ZqzUE9WwUdSaTaKloU8XltBhEpj3rPiUi0/faEvkxamst1/5nJB9ccQvNGKUFHEqmWKi2KMLOeZnagmen6l9Qo9Z4TkZrQICWRh8YNZvOuPH7z+re4qyuXxLeICjozO8vMVgALgInAsPD2Vma22MzOrIGMUo8V9547fWhm0FFEpI7p374ZNxzbhw/mruflb1YFHUekWipd0JnZKcCLwErg94QWQQDg7jnAfODcaAeU+u3VaatokpbE0f3Ue05Eou/ig7pyYI+W3P72PJZu3BF0HJEqi2SE7nfARHc/GHiijP1fA4OikkqEEr3nBqn3nIjUjIQE4/4zB5OanMDVL80kr6Ao6EgiVRJJQdcf+E8F+9cBGkaRqPlf7zldbhWRmtO2WRr3jt2P2Wu28uCERUHHEamSSAq6PUBFs9K7AFuqE0akpFenraZ7RiMGd0wPOoqI1HHHDWjLuBEdeezzbCYvzQ06jkjEIinovgR+WtaO8GrXi4BPohFKpLj33BnDOqr3nIjUit+f1I8uLRtx3csz2borP+g4IhGJpKC7DehvZp8CY8PbhpvZVcBMoClwZ1TTSb2l3nMiUtsapSbx4FmD2bB9L7/972y1MpG4UumCzt2nA8cCbfnfooh7gYeBPOBYd58f9YRS7xT3nju4ZwZtm6n3nIjUnkEd07n26F6M//Y7Xp++Jug4IpUWUR86d//C3fsCQ4CzCF2CHQH0cfcsM2tcAxmlninuPafFECIShMsP7c7+XVtwy5tzWJG7M+g4IpUSSR+6R4r/7u6z3P0Vd3/Z3ae5u5tZOjChJkJK/aLecyISpMQE4y9nDSYhwbj25ZkUFKqVicS+SEboLjOze8raYWYtCC2I6B2VVFJvqfeciMSCDukNuOe0gUxfuYVHPl0SdByRfYqkoLsYuNHMbiq50cxaE7oNWGfgmChmk3pIvedEJFaMGdSesUM78PDHi5m2YlPQcUQqFMmiiOeAXwF3m9mVAGaWCXwBtAaOcPdvaiSl1BvqPSciseT2k/vToXkDrnl5Jtv3qJWJxK5IF0X8DbgZeMjMfgN8DjQBDnP3WTWQT+oR9Z4TkVjTJC2ZB88awtote7j1rblBxxEpV0QFHYC7/wH4I3APkAQc7O7zoh1M6h/1nhORWDSsc3N+eUQPXp++hrdmrQ06jkiZksrbYWZP7uO124EFhObVFW9zd78sStmkHlHvORGJZVcd3oOJizZy8xuzGda5OR3SGwQdSeQHyi3ogEsq8fqjSz13QAWdRKy499xvT+gbdBQRkR9JSkzgwbOGcMLDX3DtyzN58ecjSUzQ1BCJHeVecnX3hCo81GdCqkS950Qk1nVq2ZA7TunPlGWbePzz7KDjiPxAxHPoRKJNvedEJF6cNqQDYwa15y8fLWLWqi1BxxH5ngo6CZx6z4lIvDAz7jp1AK2bpHLNyzPZubcg6EgiQAUFnZktM7NsM0su8XzpPh4ag5aIqfeciMSTZg2SeeCswSzP3cmd76jJg8SGihZFfE5okUNRqeciUVPce+7G4/qo95yIxI2R3VpyxaHdefSzbA7rncFxA9oFHUnquXILOne/oKLnItGg3nMiEq+uOaoXXy7J4Tevz2Zwx+ZquSSB0hw6CYx6z4lIPEtJSuDBswazN7+I61+ZSVGRLmJJcCqaQ9epKo/aDC/xrbj3nBZDiEi86pbRmFvH9OOrJbk8/eWyoONIPVbRHLrlVG3OnPpOSKWo95yI1AVnjejIpws38McPFjC6R0v6t28WdCSphyoq6C5CiyCkhhT3njt9aKZ6z4lIXDMz7h27H8c9NJGrX5rJ21cdRIMU/b8mtauiRRHP1mIOqWfUe05E6pLmjVK4/8zBnPP019zz7nzuPHVA0JGknonJRRFmdqaZzTWzIjMbXmL70WY2zcxmh/88Iry9iZnNLPHIMbMHAzsB2Sf1nhORuuagnq34+cFdeX7yCj6evz7oOFLPxGRBB8wBxgITS23PAca4+0DgfOB5AHff7u6Dix/ACuD1WswrESjuPXfGsI7qPScidcqvj+1N33ZNueHVb9mwfU/QcaQeicmCzt3nu/vCMrbPcPe14adzgTQzSy15jJn1BFoDX9R8UqkK9Z4TkboqNSmRh8cNZsfeAv7vlW9x11R0qR0xWdBV0unADHffW2r7T4GXvYJ/RWZ2qZlNNbOpGzdurNGQ8kPqPScidV3PNk343Yl9+XzRRv6ZtTzoOFJPBFbQmdkEM5tTxuOUSry2P3AfcFkZu8cBL1b0end/0t2Hu/vwjIyMqp2AVIl6z4lIfXDOyM4c0ac197y3gIXrtgcdR+qBwAo6dz/K3QeU8XizoteZWSbwBnCeu2eX2jcISHL3aTUYXapBvedEpD4wM/54xn40TUvi6pdmsCe/MOhIUsdVuqCrxF0iOppZhtXgLHczSwfGAze5+1dlHPJT9jE6J8HZHu49d/Kg9uo9JyJ1XqvGqfzpzEEsWLedP77/o2nhIlEVyQjdcmBZBY/lwDpgh5m9a2YjqxrKzE4zs9XAKGC8mX0Q3nUV0AP4fYkWJa1LvPQnqKCLWe/OVu85EalfDu/dmgtGd+GZr5bx+SLN2ZaaY5VdgWNmFwK/BDoD/wIWAwb0IjQytgx4FugJnAM0Bo5y95hebTp8+HCfOnVq0DHqhTMfz2LTzjwmXHeo2pWISL2xJ7+Qkx/5ks278nn/6oNp2Th13y8SKYeZTXP34aW3RzJC1xZoCPR091+5+1/d/WF3vwroDTQBUt39aqAPsAG4NQrZpQ5YnrOTb5ar95yI1D9pyYk8NG4IW3flc+Nrs9XKRGpEJAXd5cDf3X1T6R3ungM8ReiSKO6+EXgaGBGNkBL/Xpuu3nMiUn/1bdeUG4/vw4T563lz5tp9v0AkQpEUdK2B5Ar2JxEaxSu2mgruFSv1R1GR89o09Z4TkfrtwtFd2C+zGX94bz479xYEHUfqmEgKum+BK82sU+kdZtYZuBKYVWJzT+C76sWTumDS0lzWqveciNRzCQnGrWP6s37bXh77LHvfLxCJQCQjaNcDHwALzWw8oUURECrcTgz/fRyAmaUBZxNqMSL13KvTVqv3nIgIMKxzc04b0oEnv1jKWSM60rFFw6AjSR1R6RE6d/8SGA18CBwP3Bh+HE+o0BsZPgZ33+Pume5e1p0cpB7Zvief9+Z8p95zIiJhNx7Xh6QE4+7x84OOInVIRHeKcPdZ7n4KoRWt7cOPJu5+qrvPqvjVUh+p95yIyA+1bZbGlYf34P2568hakhN0HKkjqnTrL3cvAgqBgvDfRcr06rTVdM9oxOCO6UFHERGJGRcf1JXM5g24/e15FBTqx6hUX0QFnZl1M7OXzGwrobtCrDezrWb2bzPrVjMRJV6p95yISNnSkhP53Yl9Wbh+Oy9OWRl0HKkDKr0owsz6AF8BzYD3gXmE7hTRFzgTOMbMDnL3BTURVOKPes+JiJTv2P5tGd29Jfd/tIgxg9qT3jAl6EgSxyIZobsXKAKGuPtJ7n6Du/+fu58EDAEcuKcmQkr8Ke49d5B6z4mIlMnMuGVMP7btzucvHy0KOo7EuUgKukOBv7r77NI73H0O8AhweLSCSXxT7zkRkX3r07Yp54zszAtfr2Thuu1Bx5E4FklBlwJsq2D/1vAxIt/3njtGvedERCp07VG9aJyaxB3vzNV9XqXKIr1TxPlm1qD0jvC288PHSD1X3HtujHrPiYjsU/NGKVx/TC++WpLLh/PWBx1H4lQkBd09wH7ADDO72syOCz+uAWYAA4G7ayCjxBn1nhMRiczP9u9E7zZNuGv8PPbkFwYdR+JQJHeKeBs4B2gK/IXQbb3GAw+Et53j7u/UREiJL69OW023jEYMUe85EZFKSUpM4JYx/Vi1aTdPf7ks6DgShyK9U8SLQCdgFPCz8GMU0NHdX4p+PIk3/+s9l6necyIiETiwRyuO7d+Gv326hHVb9wQdR+JMxHeKcPcCd//a3V8OP752d40PC/C/3nNjh+hyq4hIpG4+oR8FRc4f31dLV4lMuY2FzaxTVd7Q3dXyup5S7zkRkerp1LIhPz+4K3/7NJtzRnVmaKfmQUeSOFHRCN1yYFkVHlJPqfeciEj1/eKwHrRpmsrtb82lqEhtTKRyKrr110WE7v4gUinqPSciUn2NUpP4zfF9uPblWbw2fTVnDu8YdCSJA+UWdO7+bC3mkDhX3Htu7NBM9Z4TEammUwZ14LlJK7jv/YUcN6AtTdKSg44kMS7iRREiZVHvORGR6ElIMG4b05+cHXt55NMlQceROKCCTqJCvedERKJrUMd0zhiWyTNfLmNZzs6g40iMU0En1abecyIiNeOG43qTkpjA3ePnBx1FYpwKOqk29Z4TEakZrZuk8csjezJh/nomLtoYdByJYSropFrUe05EpGZdeGAXurRsyB3vzCO/sCjoOBKjVNBJtaj3nIhIzUpNSuR3J/ZjyYYdPD9pRdBxJEZFVNCZWYqZXWxm/zKzj8xsSHh7czM7z8z0U72eUe85EZGad2Tf1hzSK4O/TFhE7o69QceRGFTpgs7MmgOTgb8DJwBHAMX3JNkK3AlcFe2AEruKe8+NGdRevedERGqQmXHLSX3ZnVfI/R8tCjqOxKBIRujuBXoBxwA9ge+XM7p7EfA6cFxU00lMU+85EZHa06N1E84b1YUXp6xk7tqtQceRGBNJQXcy8LC7T6DsW4ItATpHJZXEBfWeExGpXVcf1ZPmDVO4/e15uOvunPI/kRR0zYGlFexPBFKqF0fihXrPiYjUvmYNkvn1Mb2ZsmwT785eF3QciSGRFHTLgIEV7D8UWFi9OBIv1HtORCQYZ43oSN92Tbnn3fnszisMOo7EiEgKuheAS8zsiBLbHMDMrgVOBZ6NWjKJWeo9JyISnMQE47Yx/VizZTdPTqzowpnUJ5EUdPcBHwMfEVrt6sAjZrYeuB94E/hr1BNKzFHvORGRYB3QrSUn7teOxz5fwpotu4OOIzGg0gWduxe4+8nAOcAcYEH49VOAc919rGuGZr2g3nMiIsG76fg+uMO97y0IOorEgKRIX+DuLwIv1kAWiQPFvefGDs1U7zkRkQBlNm/I5Yd256GPF3PuyM7s37VF0JEkQLr1l0REvedERGLH5Yd2p32zNG5/ey6FRbpIVp+VO0JnZs9U4f3c3S+uRh6Jceo9JyISOxqkJHLTCX355YszeGXqKsbt3ynoSBKQii65HkHZDYQrol8P6rDi3nM3HNdbvedERGLESfu14/lJK/jTBws5fmA7mjVIDjqSBKDcgs7du9RiDokD6j0nIhJ7zIxbxvRjzCNf8tePF/O7k/oFHUkCoDl0UinqPSciErsGdGjGuBEdeTZrOUs27Ag6jgRABZ1UinrPiYjEtl8f05sGKYncNX5e0FEkABUtilgGFAF93D0//Hxfc+Tc3btHM6DEBvWeExGJbS0bp3L1kT25a/x8PlmwniP66P/r+qSiRRGfEyrgiko9l3pGvedEROLDeaO68O8pK7nznfkc1CODlCRdiKsvKirofgXscvdCAHe/oFYSScxR7zkRkfiQkpTALSf144J/fMOzWcu49BBdNKsvKirdNwNnFj8xs2fM7ICajySxRr3nRETix2G9W3NEn9Y8/PESNm7fG3QcqSUVFXT5QGqJ5xcAKvXrmeLec2cMy1TvORGROPG7E/uyt6CQP32g+7zWFxUVdEuAs81soJkVt55uaWadKnrUQmapReo9JyISf7plNObCA7vyyrTVfLt6S9BxpBZUVNDdARwKzASKV7g+GP57RQ+pI4qKnNenr1HvORGROPTLI3rQslEKt789D3etaazrKrpTxH/MbBpwGNAGuAt4FZhVO9EkaJOX5rJmy25uPL5P0FFERCRCTdKSueHYPtzw2re8NWstpwzuEHQkqUEVrXLF3bOBbAAz+znwgru/VRvBJHjqPSciEt/OGJbJ85NX8Id3F3B0vzY0TKnwx77EsUo3qHH3rirm6o/te/J5d853jBnUXr3nRETiVEKCcdvJ/Vi3bQ+PfZYddBypQRF1HDSzJDO73MzeMbO5ZjYn/PdLzUxlfx3y3ux16j0nIlIHDOvcglMHt+eJiUtZtWlX0HGkhlS6oDOzZsAk4FFgNLAT2A2MAh4HssLHSB2g3nMiInXHjcf3IdGMe96dH3QUqSGRjNDdDQwBrgLauPv+7j6C0IKJK8P77op+RKlty3N2MmX5JvWeExGpI9o1a8CVh3fnvTnryMrOCTqO1IBICrrTgMfc/VF3zy/e6O4F7v4Y8AQwNtoBpfa9Pn01ZnDaEK2IEhGpKy45uBuZzRtwx9vzKCgs2vcLJK5EUtC1BOZVsH9u+BiJY0VFzmvT13BQj1a0a9Yg6DgiIhIlacmJ3HxCXxas286L36wKOo5EWSQF3XLg2Ar2Hxc+RuJYce85LYYQEal7jhvQllHdWnL/hwvZsisv6DgSRZEUdM8AJ5vZC2Y2yMzSwo/BZvYccBLwVM3ElNry6rTVNElN4tj+bYOOIiIiUWZm3DKmH9t25/PghMVBx5EoiqSg+xPwJPAzYDqhVa47gWnAOcCT7v7nqCeUWlPce+4k9Z4TEamz+rZrytkHdOb5yStYuG570HEkSiJpLOzufjkwEPgtoeLuyfDfB7r7FTUTUWqLes+JiNQP1x3di8apSdzxzlzd57WOqFQzYDNrCHwB/N3dHye0AELqmE8WbKBDegOGdkoPOoqIiNSg5o1SuO7oXtz61lw+nLde02zqgEqN0Ln7LqAboHXOdVRhkTNpaS6ju7dU7zkRkXrg7AM60atNY+4eP589+YVBx5FqimQO3afAYTWUQwI2/7ttbN2dz4E9WgUdRUREakFSYgK3junPyk27ePrLZUHHkWqKpKD7FTDIzP5iZr1079a6pbhz+KjuaiUoIlJfHNijFcf0a8PfPl3C+m17go4j1RBJQbcM6EWosJsP7DWzvFKPvTWSUmpcVnYu3TMa0aZpWtBRRESkFv3uxH4UFDr3vb8g6ChSDZGMsv0L0FKYOii/sIgpyzZx+lCtbhURqW86tWzIJQd35dHPsjl3ZGeGdGoedCSpgkoXdO5+QQ3m+AEzOxO4DegL7O/uU8PbjwbuBVKAPOD/3P2T8L6fEmqh4sBa4Bx31x2IK+Hb1VvYlVfIaF1uFRGpl35xeA9enbaa296exxtXjCYhQYvj4k0kl1xr0xxgLDCx1PYcYIy7DwTOB54HCM/newg43N33A74Frqq9uPHtqyW5AIzspoJORKQ+apyaxG+O78OsVVt4fcaaoONIFURU0JlZupndbWazzGxr+DErvC1qY7TuPt/dF5axfYa7rw0/nQukmVkqYOFHIwv13GhKaJROKiErO4d+7ZrSvFFK0FFERCQgpw7uwOCO6dz3/gJ27C0IOo5EqNIFnZl1JzTydROQCEwAPg7//Sbg2/AxteV0YIa773X3fOAKYDahQq4f8HQtZolbe/ILmb5iCwf20OiciEh9lpBg3HZyfzZu38sjnywJOo5EKJIRuoeB5sAx7j7A3U9397HuPgA4FkgPH1MpZjbBzOaU8TilEq/tD9wHXBZ+nkyooBsCtOd/hWd5r7/UzKaa2dSNGzdWNnKdNG3FZvIKixjdXf3nRETqu8Ed0zl9aCbPfLmM5Tk7g44jEYikoDsUeNDdJ5Te4e4fESrmDq3sm7n7UeHCsPTjzYpeZ2aZwBvAee6eHd48OPye2R66Kd1/gNEVfOwn3X24uw/PyMiobOQ6KSs7h8QEY0TXFkFHERGRGHDjcb1JTjTuGj8/6CgSgUgKuh1ARcNZ68PH1BgzSwfGAze5+1cldq0B+plZcXV2NKFeebIPWdm5DMpsRuNU9YkWERFo3TSNq47oyYT565m4qH5fxYonkRR0LwLnmNmPZs6bWRpwLvDvaIQys9PMbDUwChhvZh+Ed10F9AB+b2Yzw4/W4YUStwMTzexbQiN290QjS122fU8+367eqsutIiLyAxcd1IXOLRtyxzvzyC/UbdzjQSTDMm8TupfrTDN7HFhEqOdbH+BSYC/wtpn94FKnu2dFGsrd3yB0WbX09ruAu8p5zePA45F+rPpsyrJNFBa5+s+JiMgPpCYl8rsT+/Hz56bywuQVXHhg16AjyT5EUtCVnDv3IP+7a4SVc4yFj0msUjKpcVnZuaQkJTC0s7qCi4jIDx3VtzUH92zFXz5axMmD2tOycWrQkaQCkRR0F9ZYCglEVnYuwzs3Jy1ZNbeIiPyQmXHLSf047qEveOCjRdx92sCgI0kFIrn11z9rMojUrk0785j/3TZ+fUyvoKOIiEiM6tmmCeeN6sw/s5Zz9gGd6de+adCRpByxeusvqWGTl4Zu9zVKCyJERKQC1xzZi2YNkrn97bmEOoNJLFJBV09lZefQKCWR/TKbBR1FRERiWLOGyVx/TG++XraJ9+asCzqOlEMFXT2VtSSX/bu2IDlR3wIiIlKxn+7fiT5tm3D3+PnsyS8MOo6UQT/N66Hvtu5mac5ODuyhy60iIrJvieH7vK7ZspsnJy4NOo6UQQVdPTQpu3j+nPrPiYhI5Yzs1pITB7bj0c+WsHbL7qDjSCkq6OqhrOxc0hsm07etViuJiEjl3XRCH9zh3vcWBB1FSql0QWdm/c1sbKlth5vZx2Y2zcyuj348iTZ3Z1J2LqO6tSQhwfb9AhERkbDM5g257NDuvDVrLd8s3xR0HCkhkhG6PwKXFD8xsw7AW8B+QAPgj2Z2XnTjSbSt3LSLNVt263ZfIiJSJZcf2o12zdK47a25FBapjUmsiKSgGwp8XuL52YRu6zXY3fsB7wFXRjGb1ICvlqj/nIiIVF3DlCRuOqEvc9du45Wpq4KOI2GRFHTNgfUlnh8PfO7ua8LP3wZ024EYl5WdQ5umqXTPaBR0FBERiVNj9mvHiC7N+dMHC9m2Jz/oOEJkBV0OkAlgZo2AUcCEEvuTiOzesFLLiufPje7eCjPNnxMRkaoxM24d059Nu/J4eMLioOMIkRV0XwBXhBdGPAgkE5pDV6w3sKaM10mMWLR+B7k789SuREREqm1Ah2acNbwjz2YtJ3vjjqDj1HuRFHS/BXYDrwIXA39098UAZpYInMEP59hJjMnKzgHQgggREYmKXx/bmwbJidz5zrygo9R7lS7o3H0Z0AcYDHR195tK7G4IXAH8IarpJKqysnPp1KIhmc0bBh1FRETqgFaNU7n6qJ58tnAjny7YEHScei2ixsLuXuDu37r7ilLbt7v7m+6+PKrpJGoKCouYvDRXo3MiIhJV543qQreMRtz5zjzyCoqCjlNvRVTQmVkLM7vTzL4ys8VmNiq8vaWZ3WJmfWomplTX3LXb2L6nQPPnREQkqlKSEvj9Sf1YmrOTf2YtDzpOvRXJnSI6AjOBG4AmQDdCDYVx91zgp6gPXczKCt+/dbT6z4mISJQd3rs1h/fO4OGPF7Nx+96g49RLkd4pIo3QHLojgNJ9L94Mb5cYlJWdQ682jclokhp0FBERqYN+f1I/ducX8ucPFgYdpV6KpKA7GnjY3ecDZd3rYxnQMSqpJKryCor4Zvkmjc6JiEiN6ZbRmAsP7MJ/pq1i9uqtQcepdyIp6BoBFS1haVzNLFJDZq7awp78Is2fExGRGvXLI3vSslEKt789F3fd57U2RVLQLQRGVrD/BGBO9eJITcjKzsEMRnZVQSciIjWnaVoy/3dsb6au2Mxbs9YGHadeiaSgewI4x8wuAhLD29zMmpjZX4DDgEejnE+iIGtJLgPaN6NZw+Sgo4iISB135rCODOzQjD+8u4BdeQVBx6k3Imks/BjwJPAUsDS8+VVgM3A18Fd3fyHqCaVaduUVMGPVZkb30OiciIjUvIQE49Yx/Vi3bQ+Pf5YddJx6I9LGwlcBBxIq6t4DpgCPAQe7+zVRTyfVNnX5ZvILXQsiRESk1gzv0oKTB7XniYlL2bBtT9Bx6oWICjoAd5/k7te4+4nufry7/9Ldv6qJcFJ9Wdm5JCUYI7o0DzqKiIjUI9cf04uCIudRjdLViogLOokvk7JzGNIpnYYpSUFHERGReqRzy0acOSyTf3+9krVbdgcdp86L9NZfl5hZlpmtM7O9ZpZX6qH20DFk6+58Zq/ZyihdbhURkQBcdUQPHOeRT5cEHaXOq/SwjZn9GbgWWAtMArbUUCaJkq+X5lLkMFr950REJACZzRsybkQnXpyykisO7U7HFg2DjlRnRXId7iJCCyFOcffCGsojUZSVnUtacgJDOqUHHUVEROqpKw/vwctTV/Hwx4v505mDgo5TZ0VyydWAt1XMxY9J2bmM6NKC1KTEfR8sIiJSA9o2S+OcAzrz+ow1LMvZGXScOiuSgu4DYERNBZHo2rh9LwvXb9ftvkREJHBXHNadlMQEHpqwKOgodVYkBd2vgBFmdquZdaipQBIdk5fmAqj/nIiIBC6jSSrnje7Mm7PWsnj99qDj1EmR3CliA/A8cAuw0szytco1dmVl59IkNYkB7ZsGHUVERITLDulOw+REHpywOOgodVIkq1zvAG4mtMp1KlrlGtOysnM4oFsLkhLValBERILXolEKFx3Ulb9+soQr126jnwYcoiqSVa6XoVWucWH15l2syN3F+aO6BB1FRETke5cc1I1ns5bzlwmL+Pt5w4OOU6dEMnyThla5xoVJ2eH5cz20IEJERGJHs4bJ/Pzgbnw0bz3frt4SdJw6JZKC7kO0yjUuTMrOpWWjFHq1bhJ0FBERkR+48MAupDdM5oGPtOI1miIp6K4ChpnZ7VrlGrvcnazsXEZ2b0lCggUdR0RE5AeapCVz2SHd+WzhRqat2Bx0nDojkoJuNTAA+B1a5RqzluXsZN22Pbrdl4iIxKzzR3emVeMUHvhoYdBR6oxIFkX8C/CaCiLR8VW2+s+JiEhsa5iSxOWHdueu8fOZvDSXkd00CFFdlS7o3P2CGswhUTIpO4d2zdLo0lI3QBYRkdh1zsjO/P2LpTzw4SJevmwkZpomVB1qUlaHFBU5k7JzGd29lf5hiIhITEtLTuTKw3swZfkmvlySE3ScuFfuCJ2ZHQLg7hNLPt+X4uOl9i1Yt53Nu/I1f05EROLCWSM68vhn2dz/4SIO6qHBiOqo6JLrZ4CbWQN3zyt+XsHxFt6fGLV0EpGs7NBvOKNU0ImISBxITUrkl0f25KbXZ/Ppwg0c0adN0JHiVkUF3eEA4WLu++cSuyZl59K1VSPapzcIOoqIiEilnDEsk8c+y+aBjxZxeO/WGqWronILOnf/vKLnElsKCov4etkmTh7cPugoIiIilZacmMCvjuzJr1+ZxQdz13PcgLZBR4pLlV4UYWafmNmRFew/3Mw+iU4sidS3a7ayY2+B5s+JiEjcOXVwe7q1asRfPlpEUZE6pFVFJKtcDwMqurjdGji0Wmmkyorv3zpKvXxERCTOJCUmcPVRPVm4fjvjZ38XdJy4FGnbkorK5u7AjmpkkWrIys6hT9smtGycGnQUERGRiI3Zrz292jTmwQmLKNQoXcQqbCxsZucC55bYdJOZXVjGoenAEOCj6EWTytqTX8jU5Zs5+4DOQUcRERGpkoQE49qjenHFv6bz5sw1jB2aGXSkuLKvO0W0AHqG/+5AW6BJqWMc2Eno1mC/i2o6qZQZK7ewt6BI8+dERCSuHdu/Lf3aNeWhjxczZlB7khN1/4PKqvAz5e4PuXtXd+9KqM/cNcXPSzy6uftAd7/A3VfXTmwpaVJ2DgkG+3drEXQUERGRKktIMK47uhcrcnfx+nSVFJGodOnr7gnu/u+aDCNV81V2LgMz02malhx0FBERkWo5sm9rBnVM5+GPl5BXUBR0nLihscw4t2NvAbNWbeFAXW4VEZE6wCw0Srdmy25enroq6DhxQwVdnPtm+SYKipzR3VsFHUVERCQqDunZiuGdm/O3T5awJ78w6DhxQQVdnJuUnUtKYgLDOjcPOoqIiEhUmBnXHdOLddv28O+vVwYdJy6ooItzWdk5DOmUToOUxKCjiIiIRM3o7q0Y1a0lj36Wze48jdLtiwq6OLZlVx5z127T5VYREamTrj+mFzk79vLcpOVBR4l5Kuji2OSlubjD6B5aECEiInXP8C4tOKRXBo9/ns2OvQVBx4lplSrozCzdzPY3s+4VHNPVzM6LXjTZl6zsXBqmJDIoMz3oKCIiIjXiuqN7sXlXPs9+tSzoKDFtnwWdmd0NrAcmAYvMbIqZDS3j0NHAP6KcTyqQlZ3LiC4tSEnSQKuIiNRNgzumc1Tf1jw5cSlbd+cHHSdmVVgJmNlPgJuAz4ErgXuATkCWmZ1T8/GkPBu27WHJhh263ZeIiNR51x7di217Cnj6S43SlWdfQzvXAJ+6+zHu/ri7/x7oC0wEnjWza2o4n5Rj0tJcAC2IEBGROq9/+2YcP6Atz3y5jM0784KOE5P2VdD1AV4rucHdNwPHAc8B95vZnTWUTSqQtSSXpmlJ9GvfNOgoIiIiNe7ao3uxM6+AJ79YGnSUmLSvgq7Mm6i5e5G7XwQ8CNxsZo9U4r0kirKW5jCyW0sSEyzoKCIiIjWuV5smjNmvPc9+tZycHXuDjhNz9lWELQIOKm+nu18P3A78AtBIXS1ZtWkXqzbt1vw5ERGpV64+qid7Cwp5/LPsoKPEnH0VdO8Bp5hZuZWDu98OXAt0jGYwKV9Wdg4AB/bQ/DkREak/umc05rQhmTw/eQXrt+0JOk5M2VdB9wxwA9C6ooPc/SHgdOCOaIQyszPNbK6ZFZnZ8BLbjzazaWY2O/znESX2nWVm34Zf98do5IhVWdm5tGqcSo/WjYOOIiIiUquuPrInhUXOo58uCTpKTKmwoHP3Ne7+N3efv683cvf/hkfromEOMJbQatqScoAx7j4QOB94HiA8gvgn4Eh37w+0MbMjo5Qlprg7Wdm5jO7eEjPNnxMRkfqlU8uGnDk8kxenrGLNlt1Bx4kZVV7IYGZpZnaembWJZiAAd5/v7gvL2D7D3deGn84F0swsFegGLHL3jeF9EwiNGNY52Rt3sHH7Xs2fExGReuuqI3oC8MgnGqUrVp2Vqc0I3Rmif5SyROp0YIa77wWWAH3MrIuZJQGnUkfn9GVlq/+ciIjUbx3SGzBu/468MnUVK3N3BR0nJlS31UiVr/mZ2QQzm1PG45RKvLY/cB9wGXzfG+8K4GXgC2A5UO5dfM3sUjObamZTN27cWN5hMemrJTl0SG9AxxYNgo4iIiISmCsP70FigvHwJ4uDjhITqlvQeZVf6H6Uuw8o4/FmRa8zs0zgDeA8d/9+3bK7v+3uB7j7KGAhUO5X2N2fdPfh7j48IyOjqqdQ6wqLnMlLN3FgD82fExGR+q1N0zTOGdmZ16evZunGHUHHCVxgI3RV+mBm6cB44CZ3/6rUvtbhP5sT6ov3VG1mqw3zv9vG1t35utwqIiICXHFYd1KTEnnoY43SVaeg2wh0Bb7a14GRMrPTzGw1MAoYb2YfhHddBfQAfm9mM8OP4pYqD5nZvHCee919UbRzBa24/9woLYgQERGhVeNUzh/dhbdmrWXR+u1BxwlUlQu68O2/VoQXJUSVu7/h7pnunurubdz92PD2u9y9kbsPLvHYEN73U3fvF368FO1MsSArO5fuGY1o0zQt6CgiIiIx4bJDutEoJYkHJ9S5cZyI6P6rcSK/sIgpyzbpcquIiEgJzRulcNGBXXh39jrmrt0adJzAqKCLE7NWbWFXXqH6z4mIiJRy8cHdaJqWxF8+qr9z6VTQxYms7FzMYGQ3FXQiIiIlNWuQzM8P7saE+euZtWpL0HECoYIuTmRl59CvXVOaN0oJOoqIiEjMufCgrjRvmMwDH9XPuXQq6OLAnvxCpq/YosutIiIi5WicmsRlh3bn80UbmbZiU9Bxap0KujgwbcVm8gqLtCBCRESkAueN6kyrxinc/2H9G6VTQRcHsrJzSEwwRnRtEXQUERGRmNUwJYkrDutBVnYuk8L3Pq8vVNDFgazsXAZlNqNxalLQUURERGLa2Qd0ok3TVB74aCHuVb5DadxRQRfjtu/J59vVWzmwhy63ioiI7EtaciJXHd6Db5Zv5ovFOUHHqTUq6GLclGWbKCxy3e5LRESkkn4yoiMd0htw/0eL6s0onQq6GJeVnUtKUgJDOzUPOoqIiEhcSE1K5JdH9GDWqi18smBD0HFqhQq6GJeVncvwzs1JS04MOoqIiEjcOH1YJp1aNOSBejJKp4Iuhm3amcf877ap/5yIiEiEkhMTuPrInsxdu40P5q4LOk6NU0EXwyYvDS25HqX+cyIiIhE7dUgHumU04i8fLaaoqG6P0qmgi2FfLcmhUUoi+2U2CzqKiIhI3ElMMK45qhcL12/nndnfBR2nRqmgi2GTsnM5oFtLkhP1ZRIREamKkwa2o3ebJjw4YREFhUVBx6kxqhRi1Hdbd7M0Z6fmz4mIiFRDQoJx7dE9WbpxJ2/OXBt0nBqjgi5GFd+yRP3nREREqufY/m3p374pD328mPw6Okqngi5GZWXnkt4wmb5tmwYdRUREJK6ZGdcd3YuVm3bx2rTVQcepESroYpC7Myk7l1HdWpKQYEHHERERiXtH9GnN4I7p/PWTJewtKAw6TtSpoItBK3J3sWbLbs2fExERiZLiUbo1W3bzn29WBR0n6lTQxaCs8Py50T3Uf05ERCRaDu7ZihFdmvPIp0vYk1+3RulU0MWgrOwc2jRNpVurRkFHERERqTNCo3S9Wb9tL//6emXQcaJKBV2MKZ4/N7p7K8w0f05ERCSaRnVvyejuLXnssyXsyisIOk7UqKCLMYvW7yB3Z57alYiIiNSQ64/pRc6OPJ6btCLoKFGjgi7GZGXnAGhBhIiISA0Z1rkFh/bK4InPs9mxt26M0qmgizFZ2bl0atGQzOYNg44iIiJSZ113dC8278rnH18uCzpKVKigiyEFhUVMXprLgT00OiciIlKTBnVM56i+bfj7F0vZujs/6DjVpoIuhsxdu43tewoY1V3tSkRERGradUf3YtueAp7+YmnQUapNBV0MKe4/N6qbRuhERERqWr/2TTlhYFue+Wo5m3fmBR2nWlTQxZCs7Bx6tWlMRpPUoKOIiIjUC9cc1YudeQU8MTG+R+lU0MWIvIIivlm+idG63CoiIlJrerVpwsmD2vPPrOVs3L436DhVpoIuRsxctYU9+UXqPyciIlLLrj6yJ3sLCnn88+ygo1SZCroY8dWSHBIMRmr+nIiISK3qltGYsUMzeWHyCtZv2xN0nCpRQRcjJmXnMqBDM5o1SA46ioiISL1z9ZE9KSxy/vbpkqCjVIkKuhiwK6+AGas263KriIhIQDq2aMiZwzvy0pRVrNmyO+g4EVNBFwOmLt9MfqFrQYSIiEiAfnlEDwAe+WRxwEkip4IuBmRl55KUYIzo0jzoKCIiIvVW+/QG/HT/jrwydTUrc3cFHSciKuhiwKTsHIZ0SqdhSlLQUUREROq1Kw/vQWKC8dDH8TVKp4IuYFt35zN7zVZdbhUREYkBrZumce7IzrwxYzXZG3cEHafSVNAF7OuluRQ5jNaCCBERkZhw+WHdSUtO5KEJ8TNKp4IuYFnZuaQlJzC4U3rQUURERARo1TiV80d34e1v17Jw3fag41SKCrqATcrOZUSXFqQmJQYdRURERMIuPbgbjVKSeHDCoqCjVIoKugBt3L6Xheu3q/+ciIhIjGneKIWLDurKe3PWMXft1qDj7JMKugBNXpoLoAURIiIiMejig7rSNC2Jv3wU+3PpVNAFKCs7lyapSQxo3zToKCIiIlJKswbJXHpINybMX8+sVVuCjlMhFXQBysrO4YBuLUlK1JdBREQkFl1wYFeaN0zmgY9iey6dKomArN68ixW5u9SuREREJIY1Tk3i8kO78/mijUxbsSnoOOVSQReQSdnh+XM9VNCJiIjEsvNGdaFV41Tu/zB2R+lU0AVkUnYuLRul0Kt1k6CjiIiISAUapCTyi8O6k5Wd+/2ATKxRQRcAdycrO5eR3VuSkGBBxxEREZF9+NkBnWjbNI0HPlqIuwcd50dU0AVgWc5O1m3bo/lzIiIicSItOZErj+jBN8s388XinKDj/IgKugB8FR6uPVD950REROLGWcM70iG9Afd/tCjmRulU0AVgUnYO7Zul0bllw6CjiIiISCWlJCXwqyN7MGvVFj5ZsCHoOD+ggq6WFRU5k7JzGdW9FWaaPyciIhJPxg7NpHPLhjwQY6N0Kuhq2YJ129m8K1/z50REROJQcmICVx/Zk7lrt/HB3HVBx/meCrpalpUdmkg5SgWdiIhIXDplcAe6ZzTiLx8tpqgoNkbpVNDVsknZuXRt1Yj26Q2CjiIiIiJVkJhgXHNULxau3847s78LOg6ggq5WFRQW8fWyTbrcKiIiEudOHNiOPm2b8OCERRQUFgUdRwVdbfp2zVZ27C1gtNqViIiIxLWE8Cjd0o07eXPm2qDjqKCrTcW3CxnZrUXASURERKS6ju3fhgEdmvLQx4vJD3iUTgVdLcrKzqFP2ya0bJwadBQRERGpJjPjuqN7sXLTLl6btjrQLCroasme/EKmLt+sy60iIiJ1yOG9WzO4Yzp//WQJewsKA8uhgq6WzFi5hb0FRVoQISIiUoeYGdcf04uWjVPYuH1vYDmSAvvI9UxWdg6JCcYBmj8nIiJSpxzUoxUH9Qj2DlAq6GpJVnYuAzs0o0lactBRREREJIpi4VaeuuRaC3bsLWDWqi263CoiIiI1QgVdLfhm+SYKilwLIkRERKRGqKCrBZOyc0lJTGBY5+ZBRxEREZE6KCYLOjM708zmmlmRmQ0vsX1/M5sZfswys9NK7BtmZrPNbImZPWyxcEE7LCs7hyGd0mmQkhh0FBEREamDYrKgA+YAY4GJZWwf7u6DgeOAJ8yseGHHY8ClQM/w47jaiVqxLbvymLt2my63ioiISI2JyYLO3ee7+8Iytu9y94Lw0zTAAcysHdDU3Se5uwPPAafWVt6KTF6aizsc2EMLIkRERKRmxGRBVxEzO8DM5gKzgcvDBV4HoOQ9N1aHtwUuKzuXhimJ7JeZHnQUERERqaMC60NnZhOAtmXsutnd3yzvde7+NdDfzPoC/zSz94Cy5st5BR/7UkKXZ+nUqVNEuSPVIDmRY/u3JSUp7mpnERERiROBFXTuflQ1Xz/fzHYCAwiNyGWW2J0JrK3gtU8CTwIMHz683MIvGm46oW9Nvr2IiIhIfF1yNbOuxYsgzKwz0BtY7u7fAdvNbGR4det5QLmjfCIiIiJ1SUwWdGZ2mpmtBkYB483sg/Cug4BZZjYTeAP4hbvnhPddATwFLAGygfdqN7WIiIhIMCy0KLT+Gj58uE+dOjXoGCIiIiL7ZGbT3H146e0xOUInIiIiIpWngk5EREQkzqmgExEREYlzKuhERERE4pwKOhEREZE4p4JOREREJM6poBMRERGJcyroREREROKcCjoRERGROKeCTkRERCTOqaATERERiXMq6ERERETinAo6ERERkTingk5EREQkzqmgExEREYlz5u5BZwiUmW0EVkThrVoBOVF4n3il89f56/zrL52/zl/nX3s6u3tG6Y31vqCLFjOb6u7Dg84RFJ2/zl/nr/MPOkdQdP46/1g4f11yFREREYlzKuhERERE4pwKuuh5MugAAdP51286//pN51+/6fxjgObQiYiIiMQ5jdCJiIiIxDkVdBEys+PMbKGZLTGz31Rw3AgzKzSzM2ozX03b1/mb2WFmttXMZoYftwSRs6ZU5usf/hzMNLO5ZvZ5bWesSZX4+v9fia/9nPC/gRZBZK0JlTj/Zmb2tpnNCn/9LwwiZ02pxPk3N7M3zOxbM5tiZgOCyFkTzOwZM9tgZnPK2W9m9nD4c/OtmQ2t7Yw1qRLn38fMJpnZXjP7dW3nq2mVOP+zw1/3b80sy8wG1XZG3F2PSj6ARCAb6AakALOAfuUc9wnwLnBG0Llr8/yBw4B3gs4a4PmnA/OATuHnrYPOXZvnX+r4McAnQeeu5a//b4H7wn/PADYBKUFnr8Xz/xNwa/jvfYCPg84dxfM/BBgKzCln/wnAe4ABI4Gvg85cy+ffGhgB3A38Oui8AZz/aKB5+O/HB/H11whdZPYHlrj7UnfPA14CTinjuF8CrwEbajNcLajs+ddVlTn/nwGvu/tKAHevS98DkX79fwq8WCvJakdlzt+BJmZmQGNCBV1B7casMZU5/37AxwDuvgDoYmZtajdmzXD3iYS+nuU5BXjOQyYD6WbWrnbS1bx9nb+7b3D3b4D82ktVeypx/lnuvjn8dDKQWSvBSlBBF5kOwKoSz1eHt33PzDoApwGP12Ku2rLP8w8bFb7k9J6Z9a+daLWiMuffC2huZp+Z2TQzO6/W0tW8yn79MbOGwHGEfrGpKypz/o8AfYG1wGzgancvqp14Na4y5z8LGAtgZvsDnQngB1tAKv3vQ+q8iwmN1taqpNr+gHHOythWepnwg8CN7l4Y+iW9TqnM+U8ndFuSHWZ2AvBfoGdNB6sllTn/JGAYcCTQAJhkZpPdfVFNh6sFlTn/YmOAr9y9ohGNeFOZ8z8WmAkcAXQHPjKzL9x9Ww1nqw2VOf97gYfMbCahgnYGdWeEcl8i+fchdZSZHU6ooDuotj+2CrrIrAY6lnieSeg38ZKGAy+Fi7lWwAlmVuDu/62VhDVrn+df8geXu79rZo+aWSt3rwv3+avM1381kOPuO4GdZjYRGATUhYKuMudfbBx163IrVO78LwTu9dBEmiVmtozQXLIptROxRlX23/+FEFokACwLP+qDSP59SB1kZvsBTwHHu3tubX98XXKNzDdATzPramYphH5ovVXyAHfv6u5d3L0L8CrwizpSzEElzt/M2ob/Iy++5JIA1Po3dg3Z5/kDbwIHm1lS+LLjAcD8Ws5ZUypz/phZM+BQQp+LuqQy57+S0Ogs4bljvYGltZqy5lTm3396eB/AJcDEOjI6WRlvAeeFV7uOBLa6+3dBh5LaYWadgNeBc4O6IqMRugi4e4GZXQV8QGjF1zPuPtfMLg/vr4vz5r5XyfM/A7jCzAqA3cC48GhF3KvM+bv7fDN7H/gWKAKecvcyl7nHmwi+/08DPgyPUtYZlTz/O4FnzWw2oUtwN9aR0enKnn9f4DkzKyS02vviwAJHmZm9SGgVfyszWw3cCiTD9+f+LqGVrkuAXYRHKuuKfZ2/mbUFpgJNgSIzu4bQKug6UdBX4ut/C9ASeDQ8plHg7sNrNWMd+VkrIiIiUm/pkquIiIhInFNBJyIiIhLnVNCJiIiIxDkVdCIiIiJxTgWdiIiISDWZ2TNmtsHMotLZwMwKzWxm+PGjFlGlqaATkXKFb2H2WdA5apuZ/cTM5prZXjNzM0sPOlNFzKxLOOcFQWeJhJldEM7dJegsIlHwLKFbHkbLbncfHH6cvK+DVdCJxAkze93M8s0so4Jjrg7/gBxTm9nqEjPrAfwLWAdcAZwL1KmeeiISfe4+EfjB7Q7NrLuZvR++t/cXZtanpj6+CjqR+PE8oWbg4yo45hwgB3g/Sh/zmPCjPjmM0Of5end/xt1fcPf8gDOJSHx6Eviluw8Dfg08GsFr08xsqplNNrNT93Ww7hQhEj/GE/rt7xzgr6V3mlkvQvcSfqS6BYiZNXT3Xe6eV533iVOtw39uqcqLzawBsNfdi6KWSETijpk1BkYDr4TvHgGQGt43FrijjJetcfdjw3/v5O5rzawb8ImZzXb37PI+nkboROJEuLj6D7C/mfUs45Bzw38+b2YpZna7mU0xs01mtjs8sfaC0i8ys+VmNsHMDjGzLDPbDdwT3vejOXRmdn340sHG8ByzBWb2ayvxP1aJ1y4xsx5m9oGZ7QxPGL7XzH70f4+ZnW5mX5rZdjPbFv7N9OJSxww1s7fMbHP4nKZW5jfX8GsTzOwGM1sYzr3WzP5Wcn6cmS0H7g4/XRa+fP1sBe9ZPAfsKDN7wMzWEro829TMWpjZH81sVvh8doY/vyeV8T5uZk+Z2fFmNsPM9oQ/dz8r49h2ZvZK+PO0ycyeJnS7pbLyHWhmH4eP3RH++6hyzuEIM7vPzNaFj3/FQvdmTTKze8Kfr93hz3+5l/1LvG+j8Ptlh88nNzzScEYlXjsmfOwuM9tiZm+aWd9Sx9wWzj3AzP4Z/p7YZmYvmlnrMt6zyt87IlWUAGwpMQ9usLv3BXD31919QBmP4mIOd18b/nMp8BkwZF8fTETix/PhP88uY9/PgEXuPoXQD/jLgcnA74GbCI04/cPMfl7Ga7sCbwKTgF8Bn1aQ4TpgPnBX+O8LgD+Fn5fWBJhA6Ab11wNfATcCP8hgZr8BXgUaEiomfwNMA8aUOObg8Os7ECq6/o/QPTPfMLOfVpC32KPAfcCicO43CH2OPrb/3VD+mnAOgGsJFclPVOK9HwRGhd//ZiAP6Ab8BPgQuIHQvR8bAG+ZWVmXsUcQmlT9FqFLMzsIFeffFzJmlgZ8DJwCPBV+z27Ac6XfzMwOAT4J77+H0OesO/CpmR1Yxsf/M7A/oa/jP4HTgaeBx4CR4fd4CjgJeKgSn5NHCX0O3wauCn/8RcABFb0o/LV8k9Dn6neEPrcHAVkWmt9Y2gtAJqHv838Qup/0hyW+ptH43hGJWPg+tsvM7EwACxlUmdeaWXMzKx7NawUcSOj+yBV+QD300COOHoRu/r241LYDAQd+F36eCKSW8dqPynjt8vBrx5Zx/GfAZ6W2NSzjuKcIFSCppV7rwOWljp0JfFPieVeggFDhk1zq2OL7TRuhInIikFhyP/AlsKr42HI+ZwPCWV4stf3K8PZflNj2u/C2LpX4WlwQPnZaGdlTS2YNb0sB5gIfldru4c/BgBLb2gB7gT+V2PbL8LEXltiWGP68OHBBie1Tgc1AmxLb2gFbgSllnMNEIKHE9peBIkLFfent+UDjfXxuNgN/q+Tnr0v4eTLwHbC45PsD+wGFwH9KbLst/NpPSn1P/Dy8/bJofO/ooUdlH8CL4e/ffGA1cHH4/7f3gVmECrJbKvleo4HZ4dfNBi7e12s0QicSf14AepjZyBLbziH0Q+xfAO5e6O57AcwsOXz5rxWhH349zKxZqfdcR2jEap/cfVf4fZPCv0W2IlS8NQJ6lzo8n1CxV9LnhEaNio0lVJTc5qXm/nn4fzZgENAnfH7NzaxV+OO2BN4lNELTq4LYxZc5/1Rq+98JjVz+6DJohP5eRva97l4IYGapZtaS0MjpRGBYGe8x0d3nlHj9ekKjnyU/VycBufxvpJbwx/jBnEozaxv+GM+H36f42O8Iff+MMLM2ZZxDyXl/kwgVPc+UsT0J6FTGOZS0BTjAzDru47iShgFtgcfcfUeJ3N8S+qF4vP34cv1fiz/PYc/yw69pdb93RCrF3X/q7u3cPdndM939aXdf5u7Hufsgd+/n7mXNmyvrvbLcfWD4dQPd/el9vUYFnUj8eSH85zkA4UtLPwG+dPdlxQeZ2flm9i2wh1ARsJHw3DigdEG3tETxVCEzO8HMJgO7CS3S2Mj/Coz0UoevcfeCUts2Ay1KPC++jDa7gg9b/AP38fDHK/konvP2o3lTJXQJ/7mg5EYPzUtcQui36Or40UTl8OWV681sEaGvQU447+X8+PMEsKKMbaU/V50Jfa1Kf04XlnreJfznAn5sXqljiq0s9XzLPrY3L+O9S7oe6AussND8zT+ZWVmFbEnFmcrL3RgoPX/vB+ceLqyXlXiv6n7viMQFrXIViTPuvsTMJgFnmdk1wAmEfuh/P2pjZmcRGqkYDzwArCc0WnYCoXlNpX+Z212Zj21mownNiZoE/AJYQ2i+2FBC88dKv28h+2b7PuT79/0t8E05x1S1O7sRGt2sjrI+fzcA9xL6utxOqKArBC4kNN+xtPI+V1bq72VlrcznsPSxpd+nvI9fmVw/4u6vm9mXhOZBHgVcBFxvZje7+x8qG7aMj1c6974+HzX5vSMSM1TQicSn5wlNOj+W0EjdXuCVEvvHERqlGFNy5M3Mjqjmxz2TUAF3lLvvKfG+3cp/yT4tDv85EMgq55gl4T93uvuEKnyM5eE/+wDTizeaWTKhS5qTq/Ce+zKO0PzD80puNLOLqvGey4FhZpZUapSu9CXD5eE/y2piWrytrBHBqHL3DYQWVjxtZg0J/YJxu5n9ufQl6rDlJTK+W2pfH0LzNHPK2D6/+En4a9qF0Pw4qP73jkhc0CVXkfj0MqHC6kpCc4XedvctJfYXz3n6/t94eA5XdYqJ4vctIjTnrfh90whN1q+q1wmNAt0e/mH8PbPvW6FMJ1T4XW9l3IarEm003gn/eV2p7ZcQunT4doSZK6OIUv/HWqjdzGnVeM/xhOZ+FbeowcwSKfX5d/d1hBZFnFuyhUd4bt25hBZFrKeGmFli6Xma4bmXCwktfGhUzkunEprPebmZfX+MmQ0gdEuld/3H/f1+Gf4cFLuA0CXt8eHn1f3eEYkLGqETiUPuvsnM3gVODW96vtQhbxJabPCOmf2X0ByhS4G1hFZPVtVbhC7ZTjCz5wm1JTmf0ByxKnH3ZWZ2K6F2GV+b2X8IrcQcQGhV5lh3LzKzCwmthJ1nZs8QGs1pS6gNRj9CLTnK+xhzzOwJ4DIza0pogn0/QvPZphMaRYq2NwkVqf8itGikM6HL1AuAwVV8z7+H3+MJMxtIqB3M6ZTdh+56QquaJ4fP3YDLgDR+XNhGWxNgjZm9QWiV3iZCPbQuAd4r9cvH99y9wMyuI7SA4Ssz+yehc/slsJ1QS5jSWhBqU/IGofmYVxK6hPqP8HtW63tHJF6ooBOJX88TKuhygfdK7nD358IjclcS6hm2klCfsa2Ef9BVhbt/bmbnEpqP9ACwgdBcvS8I/cCs6vvebWbZhPrA3UKohcdCStwmx92/MrP9CfUbu5TQKMx6QgVDWT/oS/sFoQLoEkKjPbmEbstzs9fMHTH+QKh1ybmEiq5FhPqx9aKKBZ277zazIwl9TS8lNEr7Rvj5rFLHTgwfewehzxnAFOBsdy/v0na07AIeITR37kRCn4eVhBbl/LGiF7r7i2a2k9DX9B5C5/gZcJO7LynjJecQmq94J6Gfaa8DVxev8g6/Z3W/d0RiXnGPJxERkbhhZrcRaqzc0d1XBxxHJHCaQyciIiIS51TQiYiIiMQ5FXQiIiIicU5z6ERERETinEboREREROKcCjoRERGROKeCTkRERCTOqaATERERiXMq6ERERETinAo6ERERkTj3/+koE6ly6U5pAAAAAElFTkSuQmCC\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",
    "with warnings.catch_warnings():\n",
    "    # Parameter is often on the boundary\n",
    "    warnings.simplefilter(\"ignore\", ConvergenceWarning)\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",
    "lbl = 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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}