{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Lesson 10: Fixed vs random effects (`fixed_regressors` / `random_regressors`)\n", "\n", "bauer's regression models let you put a patsy formula on any model parameter — for\n", "example, regress the encoding noise `n1_evidence_sd` on subject group, ISI, or\n", "session. The question this lesson answers is: **for each term in that formula, is\n", "the coefficient a single population-level number, or does it vary from subject to\n", "subject?**\n", "\n", "That is the distinction between a **fixed** and a **random** effect, and getting it\n", "right matters both statistically (it changes what the posterior means) and\n", "numerically (the wrong choice produces a poorly-identified, badly-mixing model).\n", "bauer 0.3.0 introduces an explicit API for it:\n", "\n", "* `fixed_regressors={param: formula}` — the **population-mean** design. Each column\n", " gets one coefficient shared by everyone (`group_mu`), with no per-subject offset.\n", "* `random_regressors={param: formula}` — which of those columns **additionally**\n", " carry a **per-subject random effect** (`group_mu + group_sd * offset`, i.e.\n", " partial pooling). Omitted ⇒ defaults to `'1'` (a random intercept only).\n", "\n", "The old `regressors={param: formula}` keyword still works but is **deprecated**: it\n", "silently put a random slope on *every* column — which is wrong for\n", "between-subjects covariates, as we'll see.\n" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "## Setup\n", "\n", "We import bauer from this worktree (the 0.3.0 release branch) and force JAX onto the\n", "CPU so the tiny example fits run anywhere. Everything below is sized to run in well\n", "under a minute on a laptop.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "bauer 0.3.0 from /private/tmp/bauer-rel-smoke/bauer\n" ] } ], "source": [ "import os, sys\n", "os.environ.setdefault('JAX_PLATFORMS', 'cpu') # tiny CPU fits only\n", "# Make sure we import the 0.3.0 worktree, not an older editable install.\n", "# The notebook lives in /docs/tutorial, so the repo root is two levels up.\n", "for _root in (os.path.abspath(os.path.join(os.getcwd(), '..', '..')), os.getcwd()):\n", " if os.path.isdir(os.path.join(_root, 'bauer')):\n", " sys.path.insert(0, _root)\n", " break\n", "\n", "import warnings\n", "import numpy as np\n", "import pandas as pd\n", "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "import bauer\n", "from bauer.models import MagnitudeComparisonRegressionModel\n", "print('bauer', bauer.__version__, 'from', os.path.dirname(bauer.__file__))" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "## 1. Fixed vs random, conceptually\n", "\n", "Say we regress a parameter $\\theta$ (here, log encoding noise) on a single binary\n", "covariate $x$ with a design `Intercept + x`. For subject $s$ on trial with covariate\n", "value $x$, bauer builds the coefficient vector as\n", "\n", "$$\\theta_s = \\underbrace{\\mu}_{\\text{group mean}} \\; + \\; \\underbrace{\\sigma\\,z_s}_{\\text{per-subject deviation}}$$\n", "\n", "applied column by column, where $z_s \\sim \\mathcal{N}(0, 1)$ is a standardised\n", "per-subject offset and $\\sigma$ is a group-level SD (the non-centred parameterisation).\n", "\n", "* A **fixed** effect on a column keeps only the $\\mu$ term: one number for everyone,\n", " no $\\sigma z_s$. Use this when you believe the effect is the same in the population\n", " (e.g. *the* average difference between two groups).\n", "* A **random** effect on a column adds the $\\sigma z_s$ term: the coefficient is\n", " *partially pooled*, each subject shrunk toward the group mean by an amount the data\n", " decide. Use this for quantities that genuinely differ subject-to-subject (e.g. each\n", " person's baseline noise level — the **intercept** — or each person's sensitivity to\n", " a within-subject manipulation).\n", "\n", "In bauer 0.3.0, `fixed_regressors` columns get *only* `group_mu`; the columns named\n", "in `random_regressors` *also* get `group_sd * offset`.\n" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "## 2. A tiny synthetic dataset\n", "\n", "Six subjects, two groups (3 \"control\", 3 \"patient\"). Group is a **between-subjects**\n", "covariate: it is constant within each subject — a person is a patient on every one of\n", "their trials. We make the patients genuinely noisier, plus a per-subject baseline so\n", "there is real between-subject variability in the *intercept*.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " n1 n2 choice group\n", "subject trial \n", "1 0 5 3.5 0 control\n", " 1 5 3.5 1 control\n", " 2 5 3.5 1 control\n", " 3 5 3.5 0 control\n", " 4 5 3.5 0 control\n", "trials per subject: [150, 150, 150, 150, 150, 150]\n", "group is constant within subject: True\n" ] } ], "source": [ "rng = np.random.default_rng(0)\n", "\n", "n_sub = 6\n", "groups = ['control'] * 3 + ['patient'] * 3\n", "base_ns = np.array([5, 7, 10, 14, 20])\n", "fracs = np.array([0.7, 0.85, 1.0, 1.18, 1.4])\n", "\n", "rows = []\n", "for s in range(n_sub):\n", " grp = groups[s]\n", " grp_offset = 0.5 if grp == 'patient' else 0.0 # true fixed group effect\n", " subj_intercept = rng.normal(0, 0.3) # true random intercept\n", " log_noise = -0.7 + grp_offset + subj_intercept\n", " noise = np.log1p(np.exp(log_noise)) # softplus -> positive SD\n", " for n1 in base_ns:\n", " for f in fracs:\n", " for _ in range(6):\n", " n2 = n1 * f\n", " d = (np.log(n2) - np.log(n1)) / (np.sqrt(2) * noise)\n", " p = 1 / (1 + np.exp(-2.5 * d))\n", " rows.append((s + 1, n1, n2, int(rng.random() < p), grp))\n", "\n", "df = pd.DataFrame(rows, columns=['subject', 'n1', 'n2', 'choice', 'group'])\n", "df = df.set_index(['subject', df.groupby('subject').cumcount().rename('trial')])\n", "print(df.head())\n", "print('trials per subject:', df.groupby('subject').size().tolist())\n", "print('group is constant within subject:',\n", " bool((df.reset_index().groupby('subject')['group'].nunique() == 1).all()))" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "Note the last line: `group` has exactly one value per subject. That is the defining\n", "feature of a *between-subjects* covariate, and it is what makes a random slope on it a\n", "mistake.\n" ] }, { "cell_type": "markdown", "id": "7", "metadata": {}, "source": [ "## 3. The deprecated `regressors=` keyword\n", "\n", "Before 0.3.0 you wrote `regressors={param: formula}`. That single keyword did two\n", "things at once: it added the columns to the **population-mean** design *and* gave\n", "**every** column a per-subject random effect. It still works (bit-for-bit) but now\n", "emits a `DeprecationWarning`.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DeprecationWarning: `regressors` is deprecated: it puts a per-subject random effect on EVERY term (a random slope even on between-subjects contrasts like group)\n" ] } ], "source": [ "with warnings.catch_warnings(record=True) as caught:\n", " warnings.simplefilter('always')\n", " m_legacy = MagnitudeComparisonRegressionModel(\n", " df, regressors={'n1_evidence_sd': 'C(group)'})\n", "\n", "for w in caught:\n", " if issubclass(w.category, DeprecationWarning):\n", " print('DeprecationWarning:', str(w.message).split('.')[0])" ] }, { "cell_type": "markdown", "id": "9", "metadata": {}, "source": [ "The deprecation message tells you exactly what to do: split the formula into a\n", "`fixed_regressors` part (population means) and a `random_regressors` part (per-subject\n", "effects). The legacy call is equivalent to\n", "`fixed_regressors = random_regressors = {'n1_evidence_sd': 'C(group)'}` — a random\n", "slope on the group contrast, which is the pitfall we turn to now.\n" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "## 4. The pitfall: a random slope on a between-subjects covariate\n", "\n", "Consider what a random slope on `C(group)` actually asks for. The design has two\n", "columns, `Intercept` and `C(group)[T.patient]`. A random slope means *every* subject\n", "gets a per-subject offset on **both** columns:\n", "\n", "$$\\theta_s = (\\mu_0 + \\sigma_0 z_{s,0})\\,\\text{Intercept} \\;+\\; (\\mu_1 + \\sigma_1 z_{s,1})\\,\\text{patient}_s$$\n", "\n", "But `patient`$_s$ is **0 for every control subject**. So:\n", "\n", "* For the three control subjects, the `patient` column is identically zero — their\n", " offset $z_{s,1}$ is multiplied by 0 on every trial. It enters the likelihood\n", " nowhere. Those are **non-identified nuisance dimensions**: NUTS just samples them\n", " from the prior, wasting geometry and dragging down the effective sample size.\n", "* For the three patient subjects, the `patient` column is 1, so they each carry an\n", " *extra* random-effect term that controls don't. The model now believes patients are\n", " **heteroscedastically more variable** than controls — an artefact of the\n", " parameterisation, not of the data.\n", "* Worst of all, the per-subject `patient` offset $\\sigma_1 z_{s,1}$ is, for each\n", " patient, perfectly **confounded with that subject's intercept offset**\n", " $\\sigma_0 z_{s,0}$: both are constants added to all of that subject's trials. The\n", " data cannot tell them apart. This degeneracy **inflates the posterior on the group\n", " contrast** $\\mu_1$ and under-identifies it.\n", "\n", "The fix is to make the group difference a **fixed** effect and keep a **random\n", "intercept**:\n", "\n", "```python\n", "fixed_regressors = {'n1_evidence_sd': 'C(group)'} # population mean: intercept + group\n", "random_regressors = {'n1_evidence_sd': '1'} # per-subject: intercept only\n", "```\n", "\n", "Now the group difference is a single population number (as it should be — there is\n", "only *one* control-vs-patient contrast in the world), and the genuine\n", "between-subject variability lives where it belongs: in the random intercept.\n", "bauer warns you if you ask for the wrong thing.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "UserWarning: random_regressors for 'n1_evidence_sd' includes 'C(group)[T.patient]', which is constant within subject (a between-subjects regressor). A per-subject random effect on it is non-identified for off-grou ...\n" ] } ], "source": [ "# WRONG: random slope on a between-subjects covariate\n", "with warnings.catch_warnings(record=True) as caught:\n", " warnings.simplefilter('always')\n", " m_wrong = MagnitudeComparisonRegressionModel(\n", " df,\n", " fixed_regressors={'n1_evidence_sd': 'C(group)'},\n", " random_regressors={'n1_evidence_sd': 'C(group)'})\n", " m_wrong.build_estimation_model(df, hierarchical=True)\n", "\n", "for w in caught:\n", " if issubclass(w.category, UserWarning) and not issubclass(w.category, DeprecationWarning):\n", " print('UserWarning:', str(w.message)[:200], '...')" ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "built cleanly\n" ] } ], "source": [ "# RIGHT: fixed group contrast + random intercept (no warning)\n", "m_right = MagnitudeComparisonRegressionModel(\n", " df,\n", " fixed_regressors={'n1_evidence_sd': 'C(group)'},\n", " random_regressors={'n1_evidence_sd': '1'})\n", "m_right.build_estimation_model(df, hierarchical=True)\n", "print('built cleanly')" ] }, { "cell_type": "markdown", "id": "13", "metadata": {}, "source": [ "## 5. The two parameterisations, in the model graph\n", "\n", "We don't even need to sample to see the difference. The per-subject offset variable\n", "`n1_evidence_sd_offset` carries different dimensions in the two models. In the WRONG\n", "model it spans the full regressor coordinate (Intercept **and** group); in the RIGHT\n", "model it spans a dedicated random-effects coordinate (`_re`) containing only the\n", "intercept.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "14", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "WRONG (random slope on group)\n", " offset dims : ('subject', 'n1_evidence_sd_regressors')\n", " columns with a per-subject random effect: ['Intercept', 'C(group)[T.patient]']\n", "\n", "RIGHT (fixed group + random intercept)\n", " offset dims : ('subject', 'n1_evidence_sd_re')\n", " columns with a per-subject random effect: ['Intercept']\n", "\n" ] } ], "source": [ "def offset_dims(model, name='n1_evidence_sd'):\n", " var = f'{name}_offset'\n", " dims = model.named_vars_to_dims[var]\n", " coords = {d: list(model.coords[d]) for d in dims if d in model.coords}\n", " return dims, coords\n", "\n", "for label, m in [('WRONG (random slope on group)', m_wrong),\n", " ('RIGHT (fixed group + random intercept)', m_right)]:\n", " dims, coords = offset_dims(m.estimation_model)\n", " print(label)\n", " print(' offset dims :', tuple(dims))\n", " last = list(dims)[-1]\n", " print(' columns with a per-subject random effect:', coords[last])\n", " print()" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "The WRONG model gives all 6 subjects a per-subject `C(group)[T.patient]` offset —\n", "including the 3 controls, for whom that column is always zero. The RIGHT model has a\n", "random effect *only* on the Intercept, exactly as intended. The group contrast lives\n", "purely in the population-mean node `n1_evidence_sd_mu`.\n" ] }, { "cell_type": "markdown", "id": "16", "metadata": {}, "source": [ "## 6. A quick fit: the contrast is wider / less identified in the WRONG model\n", "\n", "Now a short real fit of both models (numpyro backend, `draws = tune = 200`, 2 chains).\n", "This runs in a few seconds each on CPU. We compare the posterior on the group contrast\n", "`n1_evidence_sd_mu[C(group)[T.patient]]`.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "17", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "We recommend running at least 4 chains for robust computation of convergence diagnostics\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "We recommend running at least 4 chains for robust computation of convergence diagnostics\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "done sampling\n" ] } ], "source": [ "def quick_fit(m):\n", " m.sample(draws=200, tune=200, chains=2, backend='numpyro',\n", " target_accept=0.9, random_seed=1, progressbar=False)\n", " return m.idata\n", "\n", "with warnings.catch_warnings():\n", " warnings.simplefilter('ignore') # silence the few-chains / low-ESS notices\n", " idata_wrong = quick_fit(m_wrong)\n", " idata_right = quick_fit(m_right)\n", "print('done sampling')" ] }, { "cell_type": "code", "execution_count": null, "id": "18", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
posterior SD (contrast)95% HDI widthESS (bulk)
WRONG: random slope on group0.2691.053194.000
RIGHT: fixed group + random intercept0.2360.922373.975
\n", "
" ], "text/plain": [ " posterior SD (contrast) 95% HDI width \\\n", "WRONG: random slope on group 0.269 1.053 \n", "RIGHT: fixed group + random intercept 0.236 0.922 \n", "\n", " ESS (bulk) \n", "WRONG: random slope on group 194.000 \n", "RIGHT: fixed group + random intercept 373.975 " ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def contrast_samples(idata):\n", " return idata.posterior['n1_evidence_sd_mu'].sel(\n", " n1_evidence_sd_regressors='C(group)[T.patient]').values.ravel()\n", "\n", "c_wrong = contrast_samples(idata_wrong)\n", "c_right = contrast_samples(idata_right)\n", "\n", "summary = pd.DataFrame({\n", " 'posterior SD (contrast)': [c_wrong.std(), c_right.std()],\n", " '95% HDI width': [np.diff(az.hdi(c_wrong, hdi_prob=.95))[0],\n", " np.diff(az.hdi(c_right, hdi_prob=.95))[0]],\n", " 'ESS (bulk)': [\n", " float(az.ess(idata_wrong, var_names=['n1_evidence_sd_mu']\n", " )['n1_evidence_sd_mu'].sel(\n", " n1_evidence_sd_regressors='C(group)[T.patient]')),\n", " float(az.ess(idata_right, var_names=['n1_evidence_sd_mu']\n", " )['n1_evidence_sd_mu'].sel(\n", " n1_evidence_sd_regressors='C(group)[T.patient]')),\n", " ],\n", "}, index=['WRONG: random slope on group', 'RIGHT: fixed group + random intercept'])\n", "summary.round(3)" ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "The WRONG model's group contrast has a **wider** posterior and **lower effective\n", "sample size** — the confounding with the per-subject intercept offsets inflates and\n", "under-identifies it. The exact numbers wiggle with such a tiny fit, but the direction\n", "is the systematic consequence of the bad parameterisation. Let's plot the two\n", "contrast posteriors.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "20", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAGGCAYAAACHemKmAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAApN9JREFUeJztnQd4FNXXxt/0XkmF0HvvvfcuKAKKBRVUrBQRe8OCqIAKigUUAelFQYqANCkivfeWEAjpve/u97zX/+bbJJtkE5JsNjm/55lMdnZ25s69U94599xzrHQ6nQ6CIAiCIAiCYGFYm7sAgiAIgiAIglAURMgKgiAIgiAIFokIWUEQBEEQBMEiESErCIIgCIIgWCQiZAVBEARBEASLRISsIAiCIAiCYJGIkBUEQRAEQRAsEhGygiAIgiAIgkUiQlYQBEEQBEGwSETICqXCokWLYGVlhSNHjhS47hNPPIEaNWqUaHnOnTuH999/Hzdu3DBp/9HR0XjooYfg5+enjmP48OFqOf/ndgoL98vffvHFFybXnWFZjZWRn7lcz+3bt1XZTpw4AXPx119/oU2bNnBxcVHH8NtvvxldL7+y8phcXV1LvKwlsZ9ly5bhyy+/LNZtlmfu9boorXuEYH4OHDig2ic2NtZsZfjkk0/yvKcJpYcIWaFCwofUBx98YPQh9c4772D9+vXZln344Ydq2Zw5c3Dw4EF89tlnajn/Hz9+fImWdfDgwWo/gYGB+a7H8rHshuKQx2guIcvs16NGjYKdnR02bNigjqF79+5G1zV3WUsKEbLmvy5K4h4hlA0hy/YRISvYShUIQnZq166dq0rOnDmjlj/yyCPZlnfo0KHEq8/X11dNBdGyZUuUJShOacm+//770bt3b3MXRyhnmHpdlBcyMjKUBdrWVh7bRSElJQVOTk7F3i6C+RGLrGBW2D1Yv359ODg4oGHDhli8eLHR9dLT0/HRRx+hQYMGal0+wJ588klERETk6l4fMmQItm7dilatWqkbF3/z008/ZdvnyJEj1f89e/ZUDwdOXJ6z217f1bljxw6cP38+a93du3cbdS1geZ5//nk0atRIdVPTFaFXr174+++/jR6XVqvFxx9/jGrVqsHR0VF1w7M7vihdqIauBSxf27Zt1f+sJ325WdYlS5ao/2nNysn06dOVBZUiND/27dunxKmbmxucnZ3RqVMnbNq0Ket77icoKEj9/9prr6n95eUukl9ZDbly5QoGDRqk6rVq1ap45ZVXkJaWVqTzJD/Onj2rjo3uEPz9iy++iOTk5FzW5m+//RYtWrRQ55iXlxcefPBBXLt2LWudHj16qDq5efNm1jFxIjxeWhQNadq0qfr+8OHDWcvWrVunlp0+fTpr2eXLlzFmzBh1bumvm2+++SbXccTHx2Pq1KmoWbMm7O3tUaVKFUyaNAlJSUnZ1uP2eYw8L7gttmfz5s3xxx9/mFRfrBv9fngOe3t7q/N4+fLl2eqCk6luRPdyXfBaZfu5u7urY+ncuXOu35ILFy7g4Ycfhr+/v6pH7uvxxx9X51RB94jjx4+r+4y+DSpXrqza89atW/nWFc8bdkdXr14967i2b9+eq354TXB/bBOe52w77ofXAOH9jG2kr2++LPL+ZIipda6/x7GXqaA6zwtaRVnOWrVqqXKyXnitso718KWW90YeC89HrvvWW2/luoZNOR95b3j11VfV/zzvct6X9c8BXj98wefx0HpLeK1069ZNlZHXOK87HjtfFAwpqI25P15Lv/zyS9b+jdW3UAroBKEU+Pnnn3U83Q4fPpxr2bBhw3QbN27ULV26VFenTh1d1apVddWrV89aT6PR6AYMGKBzcXHRffDBB7rt27frFixYoKtSpYquUaNGuuTk5Kx1+bugoCC1fPHixbo///xTN3LkSLWfPXv2qHXCw8N1n3zyiVr2zTff6A4ePKgmLidjx47N2n9qaqr6rmXLlrpatWplrRsXF6e+5zbee++9rP1fuHBB99xzz+lWrFih2717t+6PP/7QjRs3Tmdtba3btWtX1nrXr19Xv+WxdunSRbd27Vrd6tWrdW3bttXZ2dnpDhw4kKue+Bs9hmU0PHYuJyyf/ndvv/12VrlDQkJ0aWlpuoCAAN0jjzyS7fcZGRm6ypUrq/rKDx4Xy9i6dWvdypUrdb/99puuX79+OisrK3XchPtZt26d2v9LL72k9n3s2DGj28uvrPpjtbe31zVs2FD3xRdf6Hbs2KF799131f54PhTlPDGGfj/VqlXTffzxx7pt27bp3n//fZ2tra1uyJAh2dZ9+umnVR288soruq1bt+qWLVuma9Cggc7f318XFham1jl79qyuc+fOqq71x8SJvP766zpXV1ddenq6+szf8PidnJzUvvXwXOI29XCbHh4euqZNm6rzm2VkGXh+sax6kpKSdC1atND5+PjoZs+erersq6++Ur/t1auXTqvVZq3L/daoUUPXrl073apVq3SbN2/W9ejRQx331atXdQXx7LPP6pydndV+eI7znP/00091c+fOzVqne/fuajJW54bn8b1eF0uWLFHnxfDhw9X5x/sK287GxkbVgZ4TJ06o+udxf/fdd7q//vpL3X9GjRqli4+Pz/cekZiYqKtUqZKuTZs2qr54X+F1MGHCBN25c+fyras33nhDbfOZZ55R582PP/6ozrfAwMBs9cN65Ho8dx988EHdhg0bVL1GRUVllevhhx/Wbdq0SZ0HvDexbS9dulSidW4M1lfjxo3VdTd9+nR1z+U2Jk6cqNu5c6daJyUlRdesWTO1Dq9hnrfvvPOOOscGDRqUbXumnI+8N/C+wnXZzjnvyzw+1inr5aefflL1+e+//6rvJk+erJs/f76qf5Zvzpw56jp58skns8pgShtzf7xeWX79/nl9CqWPCFnBLEKWooOiqVWrVtkeqjdu3FA3T8Mb7fLly9VveXM0hNvi8m+//TZrGX/n6Oiou3nzZtYy3kS9vb3VA1cPb9T8raG4zE8k8oHAm3VOcgrZnGRmZiqB2Lt3b93999+f6+HBOmD5DB8KLGufPn3uScga1g9/nxOWmaLt7t27Wct4ozYU/HnRoUMHnZ+fny4hISHbcTZp0kS9ROjbU3+Mn3/+eb7bK6isPCZ+xweKIXyA1K9fv0jniTH0+6HgM4TCksv37dunPvOBxc+zZs3Kth4frnywTZs2LWvZ4MGDc7UToajiNvbu3as+U0S5ubnpnn/+eV3Pnj2z1qtbt65uzJgxWZ/79++v6lj/wNbz4osvqvM+OjpafZ4xY4YSt4YvjmTNmjVqvxQHeviZYpnnnh4Ka/6e2ykItjuFY34UVlQV5bqgeOc6Q4cOzbYP3muaN2+uhJEeinlPT8+sl1dj5HWPOHLkiFrOF7jCwLZxcHDQjR49Otty/flkTMh269Yt27oxMTFZ4smQ4OBgtW3Dc6Uk6twYFK/8PV8c84IvC8au4ZkzZ6rlFLaFPR95X8l5X9TD4+PLy8WLF/MtO88N3p/5MsD19dePqW1MYW54zxXMg7gWCGbh4sWLqvuaXaT67lbCLjd2UxvCLiVPT08MHToUmZmZWRO7dQMCArK6k/RwObvH9LBbqV69eqqLtzT47rvvlFsD90t/NnbVs4suZ9cfeeCBB9R6ethVz+Pcu3cvNBpNiZXxueeeU/Mff/wxa9m8efNUNxu73fKCXWmHDh1S3eiGI/xtbGzw2GOPqW43tm1xw3OE9WJIs2bNsrVpYc+TvMjpB81zlOzatStrPyzPo48+mm0/3Ae7QE3ZD7u72e7sBif67uUBAwaoQSzsrg8JCVFuBH369FHrpKamqvOI3cjsbjXcN7tx+f0///yTVcYmTZqoYzdcr3///tm6YPWw+5znnh52t7NL1bB+DbfD6T/NAbRr1w5btmzB66+/rrZLX8R7pSjXBeuN3ddjx47NVk66KbBe6bLB85d1u2fPHjUQsSg+tnXq1FGuJHSZ4bXOQWGmwLZhNzr3m9PPPi+3mxEjRmT7THcg1q9hdBJCVxu6MJnqCmCMot6L2Pa8v+rPU2Ps3LlTdePzvmGI/jhyltuU87EgeH9guXJCl4H77rsPlSpVUvct3p/pUsJjvHTp0j21sWAeRMgKZiEqKkrN+fDPSc5ld+/eVT5Y9KviTcdwCgsLQ2RkZLb1eYPKCX2ciuMBWxCzZ89WIrF9+/ZYu3atenjxAcoHqbH953X89PVMTEwssXLywTB69Gh8//336gZ+6tQp5cdL37T8iImJUQLG2Ehx+pAZtm1xQuFm+JDVtynFW1HPE2PwxSPn+aNvI/1xcT+sA9Zhzv2wvU3ZD4+FYlYvZPkg79u3rxKzbA+2BcUt0QsE7p/CbO7cubn2SyFL9PtmGdmmOdejOGDZC3vN0I8y57YoBsnXX3+tHvgMQ0QBQp9NhqejCC8qRbkueMyEYilnWWfOnKmOm0KX5zDrWO/DXVg8PDzUsfMl4c0330Tjxo3Vuf/ee+/l8rM0RH/+8LzJibFlJOd1pt9GXtffvVx7Rb0X0f+8oLpkubgtQ6MFoTjlNZez3MVxDzdWR8HBwejatStCQ0Px1VdfqeuM92e9j7l++0VtY8E8yPBHwSzob1QUGDnJuczHx0etzwFcxjB8czc3S5cuVWJk/vz52ZYnJCQYXT+v46cYK+nYqRMnTlQDKn7//XdVt7Rm5rRG5oRWCmtra9y5cyfXd/oBYmwvc1Ac5wmFIh+qhg9SfRvpl3E/fCDzIciHa06MLTMGByS9++67+Pfff5Ulm0KWZeRAMIpY1ictSrS26eteb/l+4YUXjG6TA1/0ZeQgNMNBjoYUto34EDcchEY4SJPQ0saBNJwoJvXWWVrz9IN9KNzj4uJybTcv0V+U60J/TBT6eUUToWCkiGU9FjQwKz/Yc7FixQoljvnCwEFgHCjJOuexG0N//ugFd85jM2aVzSn89NvI6/ozbNfSqHNCq3ZBdclysyeH9WV4TOHh4eqaK4l7Rs66I3zZolWeg8DY+6fHWNi/orSxYB7EIiuYBT4E+cbMkc36LkrCriN2ERrCkaMUF3wAcSRtzkn/QC0MerFR3FZa3jxzChneBI1FCCC8oRpaFSl4N27cqKwGfNjeCwUdY+vWrZUbB61Vv/76q+rmoyjJD35PazPLbbhddt9SxNMyY6w7717LagrFdZ6wLnLGgiX6EcncD89ZWnWM7YcPQMPjyuuYaGnlQ5yxf1lvjLSgX05LLbtjDbtraZWmxZNdo+w2NbZvvdBhGa9evao+G1uvsAlHKGZybsPYiwGFIs8jRgOgi4k+2gP3x25bwxHqbKuc1/q9XBe0cPNljN3Axo6ZE4+DQoTxjFevXp2v9dyUc5LXO91JGF+a+z527Fie6/K64TZXrlyZbTmt+KZ2mXfs2FGVn9eaIRSSPF8Mw9yVRp2TgQMHqv1w/3nBctGqmzN5gD5KTVHC8xXlnqEXt4b3aF7Lhi5WhWnj0urpE/JHLLKCWaBVj0kGmEyAPn9PP/206hZmWJWcXVzMqEVxwe5TWhHpk8fuQt686bc4bNgwtY3CQP9B8sMPP6gHMq0XtGYZ69IqDBQQPC52QfFhyYc53+K5bYqWnPABQUvclClTlBikqGTYJH2omHuBcW/50GPdMYwNrSq0rOldAAjrky4GvFkzNI4pzJgxQ5WZoophlygOGIqKsXb5YmLMElIcZS2I4jhPeCyzZs1SD11aRvnQZzgvPqy7dOmSJZieeeYZFdaLmeroU0yBTysZw5JRyOp9kPk/BQIt9Hxx4HlPQUX4mVbWbdu2qW3poXjlOaT/3xB2h7IcFBfcB8UKBQfDMlF06MUEw2zRtYVlmzx5shK+PL/Ytcr9MVQShVVxwO3wvOc+eDz0Baeln6KL4pvQikw3FvoV81qnoGLII4bIMkZRrgueM7TG0keWLgR0MWDXNbu+T548qeb6nhK6ALEeWXZa1+gTSUspE3ewnLwn5HWP4Espz3e6TzCEFIUQ25j3L5Y5L+hywePh9cN64rnIc5PHxJd6nhsFQSHFFx92d9Ovky8MrEtug+XjfUdPadS5/lyjOOf1xbrkdUdxx655nhe8T7Cs7L5n29BNhdcFrxWGIuP1mp9/bV7oXxh5TXC7vNb5sppfzwuPj9c4623atGlKuPOcoLuJIfQxN6WNWQb6hfPaYxty30UxrAj3iJkGmQkVDGPhtwjDI3FkNkfQ16tXT4VKMTYinyNLGbaFo485OpuhcxjuiJEILl++nLUef8eR4jkxNoL3yy+/1NWsWVONVjUcMX8vUQsY2mrq1KkqbA7LyagMHPma10hhjtplqCiORGcdMMwXw9cYq7vCRi3Qj+RnPTEShLEICywvRzszbFVh+Pvvv9XIb47a5ShqRjJgqCNDChO1IL+y8pi4n5zw+5y3MFPPE2Po93Pq1CkV7ofHxVHbDIHFcDw54bnavn37rDqoXbu27vHHH1cjnvVwFDTDJ3GEPMNC5SwvI1lw2a+//pq1jCG5uE2O0uYo9ZywXp966il1jrGufH19dZ06ddJ99NFH2dZjmRnOjJEdeG7pw3Yx/JA+RBjh/l944YVc+zF2PhmDocQYpsjLy0udSwx5xH1ERkZmW++XX35RIdTYLgyHxigZxXFdMNKJIYy6wXsA2471w3riZ0YhMIRhlBhqjmGW9GHXnnjiCRVyL797BEPsMfQV25vtznplRIRFixYVWFeM6MF20h8XQ1IxrBbPV8OoJvqoBTnLbHjf5G/17coQhsZCPxV3necFz1OG22Idss4Z1YR1zrrSw9BhDF/FsFgMpcUyMByZYX0X9nzk7xltgdeKYYSJvJ4DhPcp/f2B58arr76q27JlS7bfm9rGDOPGEHsMP5cz8oRQeljxz72KYUEQLBNaEjiCl4H79QOGBMESoCWO1kBapEvan7wkuX79unIroTWVltbShhZSWpo///xz1cMiCJaGuBYIQgWEfoT0y2MXM0fmsutcECwBDmBi9z4H37D735JELF0c6H5D33R28dP1SN/dP27cOHMXTxAsEhGyglABoT/s/v37VbxbfYpFQbAEONiN/qX0yV24cCEsCfpS06+a5aa/JcM8cRAhU8PmFYJLEIT8EdcCQRAEQRAEwSKR8FuCIAiCIAiCRSJCVhAEQRAEQbBIRMgKgiAIgiAIFokIWSMwIhkDQUtkMkEQBEEQhLKLCFkjMC4hR5NyLgiCIAiCIJRNRMgKgiAIgiCYSEZGhkp5zLlgfkTICoIgCIIgmMjp06fh5+en5oL5ESErCIIgCIIgWCQiZAVBEARBEASLRFLUFgGNRiO+MYJQwbCxsYGtra2k8xUEQShDiJAtJImJibh165aE5hKECoizszMCAwNhb29v7qIIgiAIAKx0Eiw1F4why/BbcXFxcHd3z2aJvXz5snqY+fr6imVGECoIvE2mp6erkcq8D9StWxfW1uKZJQgVEd4DkpKS4OLionpqBPMiFtlCwFAbfKBRxDo5OZVcqwiCUObgNW9nZ4ebN28qUevo6GjuIgmCYAYoXg2NXIJ5EZNCEbCysir+lhAEocwjVlhBENgz279/fzUXzI8I2Xsk4/ZtpJw9WyITt50f77//Pp5++umsz7t371Yie//+/VnLxo0bh+nTp6v/a9SogQYNGqBFixaoX78+Pv3002zbCw0NxUMPPYRatWqprtPu3bvjwIEDWd8vWrRIbX/p0qVZy/744w/06NEjmw/xpEmTUKdOHTRt2hTNmzfHo48+iuvXr8Oc8NjPnDlTqvu8ceMGfHx8SnWfgiAIQsnCrJ/btm2T7J9lBHEtuAcoNK8OGgxdaipKAitHR9TevAl2lSsb/b5nz55KqBoK2fbt22PXrl3o3Llz1jIKUD1r1qxBkyZNcPv2bTRq1Ai9evVCu3btlL8PBen48eOxYsUKte5ff/2F++67Dzt37kSzZs2yBOE777yDkSNHwsHBIVt56HYxaNAgNGzYUAWKZlesVqvF4sWL1ZtrzZo1C+2HJP5HZRO2KxELpSAIglBhhezevXvx+eef4+jRo7hz5w7Wr1+P4cOH57n+E088gV9++SXXcgqys2fPqv8p2p588slc66SkpBS7T1tmTIwSsT4vvwy7oKBi3XbGrVuI/PprtY+8hGyHDh2UIGUUhaCgICVa3333XcyePRtvv/02QkJCVL1S3OakcuXKyipLfz8K2eXLl8PLywuvvfZa1jq9e/fGU089pdpoyZIlalnr1q2RmZmJb7/9FpMnT862TQpfWl45py+hXuiw3UyBbUcRzYwp586dw9y5c3Hw4EFVNu6T2+Qy/fHQOkyr8rp16xAeHq6OXd/2f//9N55//nklpnl8hmMajxw5gpdeekmJd54Tc+bMUcKfFtQ2bdpgwoQJ2LRpkzpnaH3+4Ycf8M8//6h1f/vtN1V3OUXdyy+/jB07dihxzxBNhlZxPVu3bsWbb76pjoV1PX/+fHXust0mTpyIVq1a4dSpU+q4Fi5cqKzZhHU/b9485aPt5uaGb775Rr2M5KSg42J98Lg4iPHrr79WLx054Xdsc15PPM4qVarA398fX3zxheoBuHr1qtr+lStXsGXLFvWS89lnn6kyV61aVdUVf8O2pLWeL06E/3MbPNaCjlcQBEEQLMK1gA9EPrz4kDaFr776Sgkz/USh5u3trayDhtAJ23A9TiU5MIMi1qFWrWKdTBHGFE0dO3ZUFti0tDQlaClO9INR9JZZY6GCLly4gMjIyCy3gGPHjqlt5YTL+J0hM2bMwMyZM1V0B0O4HsWJXsQag2LLmIDSs2/fPmXx5Xrc92OPPYbDhw/j+PHjSnwZWqAJ2/XQoUPYvHmzEpMUiawLukhQ9P7777/o1q0bgoOD1fqslwceeECJMoooiv4HH3xQnYskKipK7Zf747769OmjBCDXpRg0dq6ePHlSiXeKb/5PcZezzim06WLBFzFu65lnnsGoUaOyvueysWPHqpe6adOmYcyYMWo5BTHFPV/6WL8fffQRHnnkkVxlMOW4+BLC7fMYcr6E6KEbCkU2z4+1a9eq9jCE59R3332n9hETE4NXX31VCXR+7tSpkzouU8jreAVBEATBYoTswIED1YOZD2BTYEisgICArIlihw/TnBZYWngM1+NUXqF7AS1ctFzqLZVt27ZV4o7L+b0hFDfs+qclkMKPERjyG8RmLDobfz9kyBAlZnNiuA1aRemPS39ZvZ8uxSBFZ1506dJF+efqoaCkry4tkLSUUixStOnRizqWiZbQsLAwXLx4UYVI04t0CkaeO4TfUWTSUV+/P1qAKayIq6srBg8erP6nKKelm8dAKASvXbuWq8z0Kaa1lJZMClX+n7PLne3B7dBvWF9uvnjwJYuwjgzLS39lWtt///13JY7Ztvw9La4MAWVYB6YcF8PEDBs2TP1PoU7LqjEoVPXXEwVtzh4Stju3q1+Xn2mBJRT8FPGmRPTL63iF8snZ23FYeTgY83Zexs/7r+NMaBw02oLPE0Eoi7D3iQYBzgXzY9E+suyOpMWsevXq2ZZzwBGX0ceSD/8PP/wQLVu2zHM7tOBx0pPT0liWoVD96aef1AVFwUc4p8jgZDgYzNBHlt3gQ4cOVT6yFFcUbewWzgm71PldTj744ANlTTese9axvgucVtmuXbvixIkTykpoap1SSOqhWBsxYoQS5BSR+vi+XK63eBpa2ulPS4tsfkKK3xkT7Pplhn6/3J6x7eeEZWJX/J49e1Sdv/HGG8qCSmFt6n6Nwe/4Owpk/YtAUY8r53Hw2ijMdoy1T851Df/nsRvuI9UEP3KJBlL+2HUhHPN2XcHRmzHg2eHuZIfk9ExkaHSo4umIaQMaYGizyrC2lkgwguVAA9ALL7xg7mIIlh61gJYs+uhxcJIhHJVP/7wNGzYo30o+wNm9nl+YDHaVU4zoJ0t6y6L/J7utly1blmXh4vzXX39VljtaQI3BF4DnnntO+dKShx9+WHU/G1pZaV3jy8LUqVNz/Z5WONY9685wm6w7+j8aChd993Zh4TYoivXtQVcBU+A5QP9Wikm9eKfvp/47vrTw2AijMrD+9JbSosB65jH269cPn3zyiRoQR8uxIbSCUtSfP39efaa7AK29+t4C+pwalpf1ywxSfNngYDm60ej9cdkTYeyYi+O4+GKk90OPjY1VFuG8oA81reu0ghO6HHAZBWnt2rWVJZltSPHP89OQvI5XKB9kaLT46I9zeHLRYaRmaPDmwAZY91wnLB3XHiue7ohPH2iKqt7OmLjiBEZ9fxCRif9vSBCEsk50dLQaP8G5YH4s1iJLserp6Zmr65MDoDjpoYilRZEiiD6WxqAFbcqUKVmfafkrjJjlwKzixtRt0vLJY2QXcr169dQyDuLiMbB7OT9/VfqisouXfoq0eNLy+corr6joArSoUWRRyOQ1COf111/Hjz/+mPWZAoYvFxTHjRs3Vt37HJxEUcNBToQijIOy8nMvMPR1piWSYr1atWoqgoIp0KrKlxj9YC8Ke/6e0JJL30+6VegHRa1evVp1vVOQFgWKTFq+KbopNOkrSrcZdpcbvsFz0BZdCmip5Lm7atWqrO/Zc0Bxy/OQlk698KN/L8Ux3QL4O+6Drg85X1CK67j0A+boekJBznNL75aRE7YxX2Qo4Il+sJdeuNPNgdZ/bofWesObfl7HK1g+FK5P/nwY/96IxtNdayqLq6G13d7WGo0re6jp1K1YfLHtIoZ/sx+LnmyLOn5uZi27IJgCB9By/AafnRynI5iXMpOilje6gqIW6GGRKdron8eR2QVBkUF/RIqse0lRS+sSR+VT6FEomDv8llA+4AsErd7GLK2lDYUyBTPPb/3LEAeO0dpeHo+3sOS8BwjZydRo8dyvx7D3UgTeH9oYTaoYfwkyJDw+FdM3nUNscgbWP98JtXz/331FEMoiHHirHzxrzPVOKF0s0iJLX0R2TeYcwZ6X6GWX7r10HecFBSaFJkNklQS2Xl4iYoVShYMnaU2mmKV7Bi3IxSlihfLN+xvPYuf5cLw1uKFJIpb4uTvi0/ubYdq6kxj7079Y/0Jn+Lhmj1EtCIJQJoUsB2VRkOqhpYOik6Z6dgWzy5/ds/QRNIR+mxzFbSyWJgch0bWAI99pUaI7AbfJ2JslAcWsWEyFe4GuD2XFOsmIBLQyVJTjFYqPP07dxtJ/gvFizzpoW6Nw3a2ujrZ4b0hjvLrmJMb/cgSrnu2oXBAEQRAKwqx3Cj7M6DunjyhAfzn+Tz89/YAuffxPPezupy9gXtZYDlBhLEuGY6LvHoUwB5XQz1IQBEEofu7Gp+Kt9WfQpY4P+jXyL9I2/N0d8fbgRio016ztF4u9jIJQXHDsAQ1mnAvmp8z4yJYlTPWRFQShYiH3gNzwEfLEz4dxOjQOcx9qqUJs3Qtrj93CLwdu4Nen26NTbZ972pYgCOUf6bsRBEEQiszWM2HYcylCuRTcq4gl97esgqZBHpi88gTiUzOkZQRByBcRsoIgCEKRQ219vPk82tbwKrRfbJ4PJSsrTOxdF3EpGZiz/ZK0jFAmoxYw0lLO9O2CebDIqAVlidDYFMQkZU8XWlx4udijiqdTiWxbEAThXvl5/w3ciUvFm4MaFmtl+rk54qG21ZSLwcjWVdGo8v+7eAmCIBgiQvYeRWyfWXuQkmE83ee94mRngx2vdM9XzDLYPP11OTFcEoPZM1GBPmgzA+dHRkZmxQhlcH0mC2CaUgbRZ4pZppBlgHpj8T0Nt8F1CFPEXrp0KStqBBMwrFy5Ms8ycmAeU90y9i/L9/fff+Pzzz8vljrKeYyCabDNGTXkiy++KJEqY5Yvno+TJ08usP22bdumBmiWBX777TdUrlxZBoeaQHRSOubtuoxBTQJQ1cu52NvivuaV8deFu3j7tzNYM6GjpLEVBMEoImTvAVpiKWJf6VsPQcV8I78Vk4xZ2y+pfRRklWWKT4rK27dvq4xMvXr1MvogpoikeDl48CC8vLzUso0bN+Ls2bNZIjU/GMbMUDzqP+fH3bt3VYYpRpOwtv7Pk8XUDF2lBVOoMpNZafLEE0+oSZ9WuKyWs6hMmDDBpPV4LjEbWFGEbEnUB4Usz22JclIwP+27Do1Wh9Ft/8uaV9zY2VhjQrfaeOu3M9h0+g6GNpfEMIIg5MYynoplHIrYOn7mz0ZDSxKtozdv3sz1IL58+bLKnMZ0qnoRS4YOHVpi5aF47dmzJ5KTk1X2k4ceekilvf3jjz+U+P71119VZrb9+/cr6zAFLjNJvfbaazh8+LCaM4IE076+9dZbGDFihNouYwLzd4GBgejevXue+2fotccffxxhYWEq0gQZMGAAXnzxRSUiGZGClmXWCcX8Z599hl9++UUJ7mbNmuHbb79V0StyWi/nzZunrNZMk8yJx8FtXb16Va3PVLT6lLj3Cv2wuF++cLRt21YdD1PvMg0tR9AzTSLjLRMeE9MC85gYto4vN0wDy7plBA5axs+dO6dSyTJlLtuC8NiY2vbff/9Vn0eOHIn33ntP/U+hzf0eOnRIRezgevw94zOzfmfOnKnaNSeGdcY6Yi8A40OfOXNGpRBmet5atWopwcuy8kWKdbZhwwZ1rk6aNAnh4eHK+v/ss8+qYzZWH0yHzLTKLB/bjdl2fvrpJ9X7wBTMO3fuVNto0KCBshIzNTDriambr127po6BYXS+//57/PXXX2r/O3bswIIFC9R5Mn78+GJpx/IG/VcXHbiBgU0C4VEMA7zyolmQp/K//fzPixjQJECJW0EQBENEyJYjLly4oLrYjVn5jh8/jjp16hSYF5pCx9A6SxFQEBR1jP27efPmbMspGrjM0HpLUaOHWaMY45dChC4OzCY1bdo0JYApXjZt2qTEKo+JAqVz585K3Hz88cfqePz9/bMEjjEouiikKXb0wo5CVs++ffvU/l1dXVX64p9//llZq1luWgjffPNNkxJpcDs8Pr5EUAxTnOWsi3shLS1NuX2QhIQEJbQoBtl136lTJ/Tt21fVMWE5KMgoXrt166ZiLj/88MOYPn26EttsX9YnXyxGjRqlfvPhhx+qdj516pTaJl8maNmnoCWsO+6fLwS1a9dW7XXgwAElfJlS2piQzQmF5smTJ1U70/WFApjikeLS0J2F58CYMWPUywDFJ1+CKDQ56VNBGtYHexnYftw2hWxERIRaTtcVLteLcx4jxflXX32VVR4eg5OTE+6//361nOceX6ZYlxSxQt4s2n8D6Rot7m9RpcSr6fEONfDyiuNY8W8wHutYQ5pFMDu8P/KFOygoyNxFEUTIlg8efPBBZam6ePGislTS2mYMrqOH1kNaOClcKHh+/PHHrAvUmI9sfvD7ogo3WvZoPaYlTD8SlAKD1jKmSjWMVcnjo2AZPHiwErGEgpPWPWPs2rVLbZ/Q2te7d+9s31PIUewQikMKa4pY8txzz5kk0AiFH0WsvjwUTCyvYX0TWpUpzvXikAJYv3+KOmarM8ZTTz2V9T/bi+KdgpXCjdZk/q9vowceeECJM8J6ZTvr62Lu3Lnqfx8fH7WeHh47hRy3xwDftPpymV7Ics7vaPHnbyleCV8umLSEluGC4iqzjihiSceOHbPKkhO2Ma3jhnVP8U4BrheyhvVB6z4zkendVvTnPl0EaM2n5Z9QqFOE6xk9enRW3XN7tL5TyAoFk5iWiYX7r6nEBxyQWtLU8HFBj/q++Oqvy3igVRBcHMT+IpgX3u9oGBLKBnJHKAfofWQpPugqQB/Zpk2bZluHGdP4BhkTE6NcC/hQpwCihZRiwFzQwsoy0X2AllgKJYpAdu3TWpoTU/xyDckpJg3RCxliTHjqP9MPk5ZCPRRuRYGWZE6F9ZE1LCetxBTxtEizXBSkhuUxFJQc0Ec/Uv3x5UV+x25sm/rP/J/o95EfeZXLWFl4DuTXzob1kRfcDsUpr4V7PU+E7Kw9egtJqRo80LL0rFGPtK+OZ5cexbJDwXi6Wy1pEsGs0M2Krkvs6dG7rQnmQxyOyhF9+vRRlkR2peekbt26GDZsmPKTpGDUQ19Lc0ExQ8sYbwb0e6Tlj93G7C6n6KZ/ox4KG1rV6CpA6y8FMFm4cGGe26dI1Lsy0HJpuL2csHue/qS0/hEOQGJ9Eop+WqkpttnVze56Q+jjS79UQt9KiqeSEkYU/ezOooil9XL79u0m/Y7WaLpOkOjoaOUvbXjstMhT/PF8WLp0adaxlzR0d6D/rh5atunnu3jx4qxlV65cUWU2Bl0B6EbAtiF61wIunz17tmovwjktvXo4AJHHyhcU1ov+eHOWR8iOVqvDz/uvo1PtSvB1cyi16mH62l71/fDD3msqdq0gmBPehzk2gnPB/IiQLaYIA1fCE4t14jaLAt8S2WXN7tacUNTRUssubLoQ0OeUVtxXX331no6fIm/QoEGF/h19JSlcxo4dq7qSWS4O8qHFmIN5KHCbN2+uysp1KVZoqaVVkmKX3dXs7s4LdpdT6HEbU6ZMUcfLwVjGoBsDB06x25t1xG5pvfWULhh+fn6qHLSA5ozwwAFn+hBmLDctgSUFX1IolvUDnUy1OPK84E2Xx0AXCopXw+8ovPXnBkUg3VVKA7YnzwH2KHC/FOisQ7qL8LvGjRurAVd0qTAGXWkoUvl71j/PDcLzhZ95PNwOfWwNrbx0p6GLBLfP8+2ll15Sy3kOLFu2TP2W9SxkZ8/lCNyISjZLBIEHWwchKikNq46ESLMIgpCFlS6/PscKCkUMBQ8tM7TQ5JVnvSzEkRXyhuKHo9MpjujLSfHHgVB6f9biQO+aoffFFMo+dOko6oCunPeAisbjCw/hdlwqZo9sbhZ3jFnbL+JiWAL2vNoT9rZihxHMA8dzcIwADUZ6333BfIiP7D1AgUmhKZm9yiZ0T+DAJb6rMRwTB2EVp4gVhIrE1YhE7L0cicl96pnNp5hZvl5Ydgy/nwjFyDZVzVIGQRDKFiJki0HMisW0bMIu5cIODiss+kFbguVgGAJOMJ2Vh0Pg5miLLnV8zFZt1bydVVzZBX9fV64GMkhPMAcMC0nDCOeC+ZG+GUEQBCFf0jO1WHP0FnrW9zN7l/7wFlVw8W4C9l2RtNSCeaCA5bgIEbJlAxGygiAIQr78df4uopPSVexYc9O0igdq+bjgx73XzF0UoQKPo/nzzz/VXDA/ImQFQRCEfFlxOAQNAtxQvZKL2WuK7gTDW1ZR/rqX7v4XLk8QShOGBGSWSM4F8yNCVhAEQcgTRmfZeykCfcuANVZP1zo+8Haxxy8Hbpi7KIIgmBkRshYMY11yYmxQhpjSf2aSAUsZdKNPJFBcMKUukxnkx+3bt1X6Vn0QfVp4EhMTURYxpWyMytC1a1cVFkoQipvfjofCwdbarIO8cmJrY40BjQOw7lgo4lIyzF0cQRDMiAhZC4Yj8jkx05Wnp2fW55UrV2atY0r60LIqZA3TwhankP3oo49UDFFr6/Jx+lPsTp48GR988IG5iyKUM/iStPbYLXSoVQnO9mUryE3/xgFI12hVylxBECou5eNJLmSjRo0aKisV07kya9bu3btVAHg9Z86cUevoodM6s2QxwDMzIe3du9dojYaGhqqMTwxrxYkZocjdu3dx//33q8xQzLBkKCS5HwosZuJiEHmKSMKsScwI9vLLLysrMsU4hS39jhj7leX9999/VZpRJjJo2bKlsqIeOnQoK9kBLc+0RjNzV79+/dTyCRMm4Ny5c2qbzBRlLKA9hX5ematYJmb34vFxf0w/q2fevHkq1S/LxmP38TFuoeKxsVwsA+tEX+bz58+jf//+WfX33XffqeV5HaOxuLiDBw9W6/KYDTOIDR06VNWhPsWuIBQHp0PjcC0iSUUrKGvQtYCpcpf8c1OlzhWE0sLBwUGlLudcMD9l6xVbKDaCg4Oxc+dOZa2jkM2La9euKaG5detWlcWMzutMuUrLJrNiGfLoo4+qVLT6LFb6vPYUow0aNMD69esRHh6uBDFFHEUZiY2NxYEDB9T6derUwZNPPqnSji5duhRTp07FkCFD1HoUskyve/z4cSUYCddnelnyzz//YNy4cUqIs7xMuUrRSqKjo9Wc4pDbpCA1xuHDh9U2nZ2dc32Xnp6uUtD++OOPSnCyLBS8rJOrV69ixowZqmxMV8tUunnxyiuvKNHK9LlMxJCWlqYs48OGDVNCftSoUWq9yMjIrLSoxo4xp3V6zJgxWLJkiaprpmVl2lVOzCzDtuJLBIU3XwYEoThYfzwUXs52aF7Vs0xW6OCmgXh93WkViqtbPV9zF0eoIDC1tQz0KjuIkL1HmPqUkyHM3U7rI61/eqFliD6l3cWLF5GUlJTtO1owvb29leij+ClqnDqKRVOChVMQ8oJk7nlDQkJCUKtWrazP9NOkGN2+fXvWMl/f/x4cO3bswMmTJ9X/FHkUg0wFqxeyjzzySNb63CZ9OatUqWK0PLQM60UsoXCkdTkqKkr5AbM+KThpkbxw4QKef/55JbwpsE3h1q1bCAgIMPod28Pe3l6JWH1ZeDynTp1SApj74Gd9/VKIG6NXr17Kqkwr6cCBA1GvXj2cPXtWtadexBK9RTevY2RZDMvGbTz00ENZy2h95br684nHxeMThOIgU6PFhhO30a2uL2yszZPJqyAaBbqjRiVnLDt0U4SsIFRQRMjeI99//30u30QKN4ocigpaJ435nRFmhKIFzhBa3Gj5XLVqlRKzDLpcFFxdXbP+pzgy9DelwDYsCy14ixcvxr2QUzQbfjbMSW9jY5Ov365huSnmRowYoSzKrEfG7PPw8FDLKYgp4mh1ppCeNm2aSVm8aImlW4IxWBfGxD+X5fWdMdatW6dycLPcFL+0wtLFwBj5HaOhkOX+KXzzO0a2q5OTk0llFISC+PtKJKKS0tGjDLoV6OE1SV/ZhfuuIzw+FX7u/3+vEYSSgsaN3r17K4MN3cQE8yI+svfIs88+q0SL4fThhx+q74KCgnJ9x0kPu9Jzfqe3LNJyx20XB7QO0wpKi59eLOuhbymtsoZd2fRNNSYwaaGcM2dO1jK9a0GfPn2y/GK5jC4GtEoWBF0Z4uLi8hVm7JqvWvW/nOpz587N+o4vCXyI0Q/2iy++UEKPVuSCtqm35BqDXfZ0A6A4JrRA01WCIrRHjx7KB1XvDvDLL78Y3QZFOt0Q6EdLFwe6JrA+69evr4Tp6tWrs9bltvI7RkP4e4pwwxcOWtL1LhWE7gw8PkEoDjaevI0gLyfU9jV/7Nj8oNCmxXi1DPoSSgne53n/LsuDqSsSYpG9R9j1n1f3Py2R+m7fvMRJXui77YsDduNTVFFc0XXB0I2A3fi0HtNnlZZKWgJZ5l9//TXXdiiAX3rpJeUfRCvv8OHDlTX666+/VoOs+GbKkFZvvfVWlltBfjzzzDPKn/Tzzz/HJ598kut7itLp06erbVWrVi3b4K3Tp0/j9ddfVwKW+6SfKffPGwvrlf6itNpu2LAh2zZ5/KxbdtPzOAyh0Fy7dq3y+aXLB9uPwtPFxUUJRFp96ZPK9qZQp+U0J7R80+2A/rusI+7r559/Vv///vvvKloCj4ki/IUXXlAvK3kdoyH8/caNG1V0Agp37ofb1rcTfZoJj1sQ7pXUDA22nb2LIc0CTe6JMBeuDrYqNNiyQ8F4rnttWJdRNwhBEEoGK52+n1vIQt+9S8sexZQeWs9o2aSF07C7XLAs6LaxZ88efPPNN4X6HX1S3dzc1P90+aBFNC8/2dKGop4vJRwoJpQcFeUesO1sGJ5ZchTfjGmFat65B0aWNS7cicera0/hl6faobsM+hJKmGPHjilXMPai5mesEkoHscgKFQ66bTBkGC25hYklS7HIqAC0WlPIMLpBWYEREmgJFoTi4I9Td9QgKksQsaR+gJsq78rDwSJkBaGCIUJWqJDQRaKwFNaCW5rQHUIQioOUdA22n7+LES2NRxYpi9D9oXdDf5WyNjopXcWYFYSSgpFoOIaCc6GCD/Zi4H2GKKI1iTei3377Ld/1ObKb6+Wccg7eoZ8jA9IzWDHnHHwkCIIgFMyui+FKzHata1lxWfVJG5hSVxBKEg5+ZuIcwyg7gvkwq5DlgBoOomHGpMLAmJr6+K2cDOOOHjx4UGV84uAfxjblnF3JeWVLEgRBEP6fzafvoJaPCyp7WlYoNw8nO7Sv5Y2Vh0OyQhwKQknAqDlMYiNxu8sGZhWyDBbPGJsMoF8YGJSewd/1E2OT6vnyyy/Rt29fvPHGGyqcEueM98blxYXcJAWhYkK/6vIerWDnhXCV+tUS6dPQHxfvJqjUuoJQUjAsI0NRci6YH4v0kWVOeo4eptvA22+/jZ49e2azyDJEkSHM1FQcQpZpQOnKwFipDH1U1sPSCIJQfC+vHOTHa58DBA2TVZQn/r4cieR0DTrV/i/rnKXRsqoXfFztlVW2WVDZTKsrCEIFFrKM38nA+wx7wcD1jGtKayt9Z/WxUcPCwuDv75/td/zM5XnBbXEyDL9lDFp+meSA3Qn6uJ2CIFQcmJSC8X4LE+3Ckthy5g6qejuhqoVEK8gJEyPQV5bJHN4Z0giOdv/fWycIQvnEooQsA90bJhGgszWzOTFAvGGQ/5yW0oLSi86YMSNXmtm8oHM3fXKZjUkQhIoDX2SZmKK89sSkZ2qx49xdDGxiPMGLpdC7gb/K8rXjPBM6VDZ3cQRBKGEsSsgag5mWDIPS02c2p/WVfiw5rbSG0I+WjtuGFll9ytC8HmiGfrmCIAiWzsFrUYhPzbRY/1g9Vbyc0DDQDauP3BIhK5QIPj4+eP7559VcMD8WL2SPHz+eLUUsrbTbt2/P5ie7bds2dOrUKc9tMEwXJ0EQhIrKn2fDEODuiJo+LmYrgzY1FSlnz0Abn6D+h04Hx3p1YV+jJqxsbQtllf129xXcjU+Fv3v5zcAmmAe6F5XluOIVDbMK2cTERJXmUw9TP544cQLe3t7qRKGlNDQ0FIsXL1bfc8BWjRo10LhxYzXwgpZYxozlpGfixInKzWDmzJkYNmyYym+/Y8cO7Nu3zyzHKAiCUNbRanXYfvYuOtepVOquEzrokHzkKBJ370LK0WPQpaf/94VeuGZmwsrREU4tWsBz1Eg4VK9R4Da71PHBD39fw7pjoXiuR+0SPgKhopGcnKzi1zMyEv3mhQosZI8cOZIt4oC+e3/s2LFYtGiRihEbHByc9T3F69SpU5W4dXJyUoJ206ZNGDRoUNY6tLyuWLFCRTN45513ULt2baxcuRLt27cv5aMTBEGwDE7cikVEYho61Cpdt4L0W7cQtWABUk+fhl1QEFx794ZT06aw8fRUFlhdZibSb99GxtWrSDp8GLdfmQqXrl3gPWYMbH3/S4BgDBcHW3SsVQlrj93ChO61yq1fs2AeKGI56Pzo0aNo1aqVNIOZsdJJUNRc0EfWw8MDcXFxcHd3N0e7CIIglBoztpzHyn9DsOjJdmrkf0mj02kRu3oNYtesgY2XFzyGDoVjw4b5/0ajQfLhw0j46y9lpfWdMhnOLVrmuf7RmzF4f+NZbHyxC5oGeZTAUQgVlWPHjomQLUOUzxgygiAIgknQlvHnmTC0reldKiJWm5qC8M+/QOyqVXDt2RO+kycXKGKJlY0NXDp0gO+UKbCrWhV3P/oYsevXK9cEY7So6glvF3tllRUEofwiQlYQBKECczUiETeiktGhZsm7FWRGRuDOm28h5eRJeD8xFu79+sHazq5Q27BxcoL3k08qN4SYpUsR9d33RsUsRXn3er74/USoCi0mCEL5RISsIAhCBebPs3fhaGeN5lU9Sl7Evv0ONAkJ8Hn+eTg2bFTkbdHnlSLYY9QoJGzfjuiffjYqZnvV90NMcgb2XIq4x9ILwv/DhChubm7lNjGKpWHx4bcEQRCEosMkCEzt6mBbcrGxM6Micefd95Sfa6Vnn4Wtl1exbNelTRvoMjIQv349rOzt4PXoo7DC/7tH1PBxQW1fF6w9GoK+jfKOJS4IhaFFixZ5ZgAVSh95nRAEQaighCek4kRILNrX9C6xfWhiY5SI1aanF6uI1ePasSPchw5F3PrfEL9la67ve9T3w18XwhGXLNkYBaE8IkJWEAShgrLrQjgYmapNjZIRskxqEPbJJ9CmpMCnBESsHteuXeHSpQuif/5ZJVQwpHtdX2i1wB+nb5fIvoWKx7lz51T4T84F8yNCVhAEoYKy/dxdNAx0h4dT4QZcmYJOq0HEnDnICL2tBmfZepec1Ze4DxkC+5o1Ef7FLOWPq8fLxR4tq3li7dHQEt2/UHFITU1VIpZzwfyIkBUEQaiApKRrsO9yJNqVkDWWA7CSjx2D15gxsK9cGSWNlbU1vB55RCVSuDtzJrSZ/+9K0LO+H44Fx+BGZFKJl0MQhNJFhKwgCEIFZP+VSKRmatGuBPxj47dtQ/yWLfAYPhyODRqgtLBxcYHXY48hPTgEMb8uy1revpY3XOxtsO64WGUFobwhQlYQBKECsuP8XVTxdEKQV/Hmik+9cAFRCxfCuVMnlcCgtLGvUgXu/fsjfsMGpJw+pZYxIkOn2j5Yd+yWSgAhCEL5QYSsIAhCBUOr1eGv8+FoW8xuBZnRUbj7+eewr1ZNpZ01Fy7dusG+Th1EfD0XmoT/wiT1bOCHWzEpOHIzxmzlEsoHtWrVwu+//67mgvkRISsIglDBOHM7DhGJacXqVqDLzMDdzz5X/ytfVZuSi0trSsIEz9GjoUtNReQPP6hljSu7w9/NQVllBeFe8PT0xH333afmgvkRISsIglDB2HE+HK4OtmgY4FZs24xe9AvSr1+H92OPwsat+LZbVGw9POA+fDiSDxxE0qF/YG1lhe71/fDHqTtIzdCYu3iCBRMWFoYZM2aouWB+RMgKgiBUMP46fxetqnnB1qZ4HgGJf/+tBncxMYF91WooKzg1bw6HRo0Q+eMCaJIS0bO+LxJSM5VbhSAUldu3b+PNN99Uc8H8iJAVBEGoQITFpeLs7Xi0rVE8yQnSb4Ugcv58OLVqZZbBXQW5GHjcfz90KSmIWbJEDWyr7++mUtYKglA+ECErCIJQgfjrwl3YWAGtq9+7kNWmpiD8s89h4+UFjwceUMKxrEEXA7eBA5GwfQdSzpxRg772XIpEZGKauYsmCEIxIEJWEAShAsFudWbzcnO8t2xeOugQOf87ZEZFwfuxx2Btb4+yCi3F9rVqIur779GlpodKy/v7CekWFoTygAhZQRCECgIHOR24ElksYbcS/vwTSfv2wXPECNj6+qIso1wMhg1HRlgYdDv+VNEa1h6V6AVC0WC0ggcffFCiFpQRRMgKgiBUEA5ejVLZvO5VyKZdvapS0Dp37KgGVFkCdoGBcOnSBbGrVqNboCPO3YnHhbD/YswKQmFg/NjVq1dLHNkygghZQRCECuQfG+jhiCAvpyJvQ5OYiPDPP4dtYIBZkx4UBbc+fWDl6Iga29fCw8kO645Jylqh8KSnp+PWrVtqLpgfEbKCIAgVAKZm3Xk+XA3yKuqgLJ1Oi8i5c6FNSoLXo4/BytYWloS1oyPcBw9G+j8H0dETWH8sFJkarbmLJVgYZ86cQdWqVdVcMD8iZAVBECoAF+8m4HZc6j25FcT//juSjxyB5+hRsPUqnvBdpQ1dITjwq+WhzSq72d+XI81dJEEQ7gERsoIgCBWAnRfC4WhnjaZVPIr0+9RzZxG9bDlce/aEY8NGsFRojXYfMhQB186imr0Ga2TQlyBYNJbVLyQIgiAUOexWiyBP2BUhm1dmTDTuzpoN+xo14Na/v8W3gH1QEFxat0brq4fxh8YWcckZ8HC+t3BkgiCYB7HICoIglHNiktJxPDgGbYrgVqDLzET4rNl0kIXnmDGwsi4fjw3XgQPRJuw8NBotNpySmLKCYKmUjzuSIAiCkCd7LkVAqwPaFCGbV/SyZUi7eBFeY8bA1s2t3NQyj6Vylw5oFHUdK/dfMXdxBAuiRYsWSE1NVXPB/IiQFQRBqAD+sXV8XVHJ1aFQv0v65x81wMt98CA41KyJ8oZL167oEHcDZyJScelugrmLI1gI1tbWcHBwUHPB/EgrCIIglGMYXmr3pXC0rlE4a2x66C1EzJ0Lx2bN4NKlK8oj1nZ26NS+AVzTk7Fs81FzF0ewEC5duoQePXqouWB+RMgKgiCUY46HxCI+JRNtq5vuH6tJSUH4Z5/DxsMDniNHFjnurCXg1qol2ieGYP25KKRnasxdHMECSExMxJ49e9RcMD8iZAVBEMq5W4Gnkx3q+ruatL4OOkR++w0yIyPh9fjjsHYonDuCpcHBa32aBSHOxhGbV/9l7uIIglDSQvb69esoLvbu3YuhQ4eicuXK6o3/t99+y3f9devWoW/fvvD19YW7uzs6duyIP//8M9s6ixYtUtvKOdExWxAEoaLBbF6tqnnB2kSravzGP5B84KCyxNr5+qIiULdpHdRIj8XKvZdUlAZBEMqxkK1Tpw569uyJpUuX3rM4TEpKQvPmzTFv3jyThS+F7ObNm3H06FFVDgrh48ePZ1uPIvfOnTvZJkdHx3sqqyAIgqURGpuiMnq1MdE/NuXsGUQvWQKXHj3g1LQpKhK9arrjkFs1XF6Vv0FFEAQLF7InT55Ey5Yt8corryAgIADPPvss/v333yLtfODAgfjoo4/wwAMPmLT+l19+iWnTpqFt27aoW7cuPvnkEzXfuHFjtvVogWXZDCdBEISKxq4L4bCxAlpWK1jIZkZHq3ix9jVrwn3AAFQ0ujcMhK0VsGLjYWjT0sxdHKEMU61aNfz4449qLligkG3SpAlmz56N0NBQ/PzzzwgLC0OXLl3QuHFjtTwiIgKlhVarRUJCAry9sw9ioAN29erVERQUhCFDhuSy2AqCIFQU/9hGld3h6mBbcNKDL76gFUDFiy0vSQ8Kg4udNTr42GKzT2NELVtu7uIIZRgfHx+MHz9ezQXzU+S7la2tLe6//36sWrUKM2fOxNWrVzF16lQlHh9//HHVnV/SzJo1S7knjBo1KmtZgwYNlJ/shg0bsHz5cuVS0LlzZ1y+fDnP7aSlpSE+Pj7bJAiCYMmkZmhw4GokWpsQrSB6yWKkXb0Kr0cegY2raYPCyiN96ngjzKUSdq3cCk1ikrmLI5RRIiMjsWDBAjUXLFjIHjlyBM8//zwCAwOVJZYilmJ2586dylo7bNgwlCQUqe+//z5WrlwJPz+/rOUdOnTAo48+qnxvu3btqoR2vXr1MHfu3Dy3NWPGDHh4eGRNVatWLdGyC4IglDT/XItCaoa2wGxeSfv3If6PTXAfMgQO1atX6IZp6GWHyo5W2OzbFNGLFpm7OEIZJTg4GE8//bSaCxYoZClamzZtik6dOuH27dtYvHgxbt68qXxda9asqayf33//PY4dO1YyJQaUeB03bpwSqX369Ml3XWbeoE9tfhbZN954A3FxcVlTSEhICZRaEAShdP1j/d0dUM3bOc910m/dQsS38+HUsgVcOnas8M3D8RW9q7tif+VmuLF0BTRxcRW+TgSh3AnZ+fPnY8yYMepNhOGy6IOaM00bHaAXLlyIkrLEPvHEE1i2bBkGDx5c4Po6nQ4nTpxQluO8YKo5RjownARBECwV3vfoH8uwW3klM9CmpiL8889h4+kJjxEPluukB4WhRxUn6Kytsd2/KaLEKisIZZ78RwAYYfv27Uqo5hSvvHHSksnv7O3tMXbs2AK3xUFZV65cyRajlqKTg7e4HVpK6aZAq69exNL/9quvvlIuBBxoRpycnJRLAPnggw/Ud4xmQF/Xr7/+Wm3zm2++KeyhCoIgWCRXI5IQEpOCJzvXzDvpwY8/IjM8HD4vvwxre/tSL2NZxcPBBu38HbFV2xP3//IpvB9/HLZehUvvKwhCGbbI1q5d26iDc3R0tHItKKyfLUN5cSJTpkxR/7/77rvqMweMGfqg0GUhMzMTL7zwgrKw6qeJEydmrRMbG4tnnnkGDRs2RL9+/ZQQZvzZdu3aFfZQBUEQLNatwN7WGk2r/PeCn5PEnTuRtHs3PEY8ADuDMQbCf/St5oxgK2ec8aiG6J9+lmoRsuHq6oru3buruWB+rHQ0pRYCWmJpCTUcYEXoJ9uoUSMVRcDSoSWXFl76y4qbgSAIlsbDP/yDdI0W7w9tnOu7tOCbuPPa63Bq2RKeI0aYpXxlHT4WX94bgabJYZjy59eos2M7bCtVMnexBEG4F9cCWksJ/ahoMXV2/v8BBBqNBocOHUKLFi1M3ZwgCIJQAiSkZuDwjWiM75K7h0ybloqI2bNh4+0N9/vuk/rPAz7n+lR1xorLvhhn54yohT/Bf9qrUl/Cf9eRVouMjAzY2dnlcrMUSh+TW4BJBTjxTfX06dNZnzlduHBBhbti/FZBEATBfOy7HIlMrQ6ta+SOHxv98yJk3g1X8WKt7ezMUj5LoWfQf8aavb0fRszy5ciMiTF3kYQyAsfdMEY954IFWWR37dql5k8++aQabCVd7oIgCGUPRiuo5u2EAHfHbMuTDh5Awvbt8BgxAnb+/mYrn6Xgbm+NjgFO2BBbC4N1QPSiX+A3eZK5iyUIQg4KbRNnWloRsYIgCGUPrVaH3RcjcmXzyoyMROT87+DYrCmcZeCryfSr5oyQZC0u9B+FmKVLJa6sIFiqRfaBBx5QbgMUsPw/P9atW1dcZRMEQRAKwdnb8YhITENbg2xeOp0WEcxsaGcn8WILSX1PO9Rws8VG55ZolLkC0YuXwPelF+WcFARLs8hyBL8+WLZhKldjkyAIgmA+twIXexs0DPz/pC7xmzYh9cwZeI4cCRsnJ2maQsDnHq2ye+5mIqnPYEQvXgxNQoLUoSBYcvitioCE3xIEwRIZ9s0+uDrY4vUBDdXn9JBg3H51Gpzat4enRCkoEimZWjy7KxyPVLfDA1+/At+XXoLPM08Xb8MJFkV6ejrCw8NVGFImgBIszEc2JSUFycnJ2eLHfvnll9i2bVtxl00QBEEwkajENJwKiUOb//nH6jIzEfH117CpVAnuAwdKPRYRJ1trlbZ29a1M2PXohehFi1R6X6HiQvEaFBQkItZSheywYcOyUsYyixYzZs2aNUstnz9/fkmUURAEQSgADvJi91rrav/5x8b9/hvSb9xULgUSauve6F/NBbFpWhxqNwia2FjEyliQCs21a9cwcuRINRcsUMgeO3YMXbt2Vf+vWbMGAQEByipLcfv111+XRBkFQRCEAth1MRx1/Vzh5WKvsnfFrFoN1+7dYV+1qtTdPVLF1RbNfeyxMswKLh07InrBQmXxFiomNOJR/3AuWKCQpVuBm5ub+p/uBIxiwMwWHTp0UIJWEARBKF0yNVrsvRSBNtW9oNNkInLePNj6+MCtb19pimJiYDUXnIpKR0if+5Fx+zbiN2+WuhUESxSyderUwW+//YaQkBD8+eef6Nevn1pOx2eJLysIglD6HAuORXxqJtrU8Eb8xo1Iv34DniMfhJWtyTlvhAJo6ecAf2cbrIx3gVOrVoj84QfotFqpN0GwNCH77rvvYurUqahRowbat2+Pjh07ZllnW7ZsWRJlFARBEAoIu+XpbIfq2kTErFwFly5dYF+1mtRZMWJjZYX+1Zyx+WYSMobcj/QrV5G4d6/UsSBYmpB98MEHERwcjCNHjmDr1q1Zy3v37o05c+YUd/kEQRCEAth54S5aVfNC9I8/wtrVFW7/6ykTipdeQc5K0P6OADjUq4eohT9JFVdAKleujE8++UTNBfMjcWSNIHFkBUGwFEJjU9D5052YWF2L2ou+gvdTT8GxQQNzF6vc8sOZOByLSMXGKrcRO+sL1Fi9Ck5Nm5q7WIJQYSm0RTYpKQnvvPMOOnXqpPxla9WqlW0SBEEQSo9dF8JhbQVU+X0pHFs0FxFbwgys7ozIVC3+9msE28BAscpWQBitYMOGDRK1oIxQ6JEA48ePx549e/DYY48hMDAwK3WtIAiCYB7/2LpIglNKEtyH3idNUMJUdbNDs0r2WHIpEd2HDEH0woVIDwmRMGcVCMaPZez8o0ePolWrVuYuToWn0EJ2y5Yt2LRpEzp37lzhK08QBMGcpGZosP9yBPpfOQa3/v1h+7/QiELJMqiGCz49GoOrvTqhkutKRC/6BQHvvC3VLgiW4Frg5eUFb+//UiAKgiAI5uOfq5FI0+jQ3DpRBeoXSodWvg4IdLbBkqsp6gUidu1aaOLipPoFwRKE7IcffqhCcDExgiAIgmA+tmw6CO+UONQf2AtW1oW+nQtFxNrKCgOru2BbcDKSu/WFTqNB7OrVUp+CYAmuBbNmzcLVq1fh7++vYsna2dnlSmErCIIglCyZMTHYHZyIZjYpcKwpUQpKm55BTlhxOUGlrX28SxdEL1kK77FjYZXjmSiUPxwdHdGoUSM1FyxQyA4fPrxkSiIIgiCYzOG5C3HXqSHGNgqUWjMDTrbWKq7sqssJGDdgMBJ37ULC9u1wHzRI2qOcQxF79uxZcxdDKKqQfe+99wr7E0EQBKEYSbtyBTsOX4Nd4/poHuQpdWvGUFybbyRhq6YSujdtiqhFv4iQFYRSxrqoMdQWLFiAN954A9HR0VkuBaGhocVdPkEQBMEAnU6Hu59+isNBTdHUxxEONhIC0Vz4O9uijb8DfrkQB9dBg5F66hSSjx+X87Wcc+LECbi7u6u5YIFC9tSpU6hXrx5mzpyJL774Iisg8Pr165WwFQRBEEqOxD17EP7PUZz1qIpWfuKjZ26G1HDF9fhMHAtsCNvKlRH9y2JzF0koYbRaLRISEtRcsEAhO2XKFDzxxBO4fPlyNkfngQMHYu/evcVdPkEQBOF/6DIyEP7pTJxp0wcaWKkwUIJ5aehlh1rudlhyMQHu/fsrP9mMu3elWQShrArZw4cP49lnn821vEqVKggLCyuucgmCIAg5iFm9Guk3b+JYsx6o6moLP+dCD3MQihlmtxxUwxn77qTibuuusLK3R8yKFVLPglBWhSytsPHx8bmWX7x4Eb6+vsVVLkEQBMEATWIiIufOg1P37tgXZy3W2DJE5wAneDpY49ebGXDt0QOxK1dBm55u7mIJQoWg0EKW+YWnT5+OjIyMrLfR4OBgvP766xgxYkRJlFEQBKHCE7VgAbRJSbgzcBSi07QiZMsQdjZW6F/VGb9fS4Sud39ooqORsGWLuYsllBANGjTA0aNH1VywQCHLAV4RERHw8/NDSkoKunfvjjp16sDNzQ0ff/xxyZRSEAShApMRFobonxfBfcgQ7Eu0g4utFep72Zu7WIIB/ao5Q6PTYX2SG5xaNFcJEoTyibOzM1q1aqXmggUKWYac2LdvH9auXYtPP/0UL774IjZv3ow9e/bAxcWlUNvi4LChQ4eicuXKyrL722+/Ffgb7qd169bKxaFWrVr47rvvcq3DsjFgsYODg5ozooIgCIKlEvH1XFg5OsJj+HDsDk1Bcx8H2FpL2K2yhIeDDboEOmHZxXg49R+I1DNnkCLhmcol7IV+4YUX1FwwP0VOzt2rVy9MnToV06ZNQ58+fYq0jaSkJDRv3hzz5s0zaf3r169j0KBB6Nq1K44fP44333wTL7/8shKueg4ePIjRo0fjsccew8mTJ9V81KhROHToUJHKKAiCYO7kB3G//QbPESMQbeWAs9HpaO0n0QrKIoNquCAsWYP9PvVhGxCAmOXLzV0koQSIjIzEt99+q+aC+bHSMbp2AXz99dcmb5DCskgFsbJSltP8UuC+9tpr2LBhA86fP5+1bMKECUqwUsASilgORtti4J80YMAAeHl5YbmJNxX+3sPDA3FxccoCLQiCYC5CXnxJBdqv8tVXWB+circPRmFBLz9lARTKHu8dioK9jRXmp/6L2BUrUGfvHth6eZm7WEIxwgRQ7BmmnyxdDATzYlLsljlz5mT7TB/Z5ORkeHr+lxqRSRHoK0K/2aIKWVOgWO3Xr1+2Zf3798fChQvV4DM7Ozu1zuTJk3Ot8+WXX+a53bS0NDXpMRaVQRAEobRJOXkSiTt2wOell2BlZ4c9t2JQ19NORGwZZlB1F3x+PAY3unWF54oViFu7FpXGjzd3sQShYrsWsEtfP3FAV4sWLZRVlOlpOfF/vpV8+OGHJVpYxqn19/fPtoyfMzMzs0z8ea2TX4zbGTNmKAusfqpatWoJHYEgCIJpsLMsfNYs2FWvBpcuXZCu0WH/nRSJVlDGYcpaPycb/BqigXOnTohZvgI6jcbcxRKEckuhfWTfeecdzJ07F/Xr189axv9ptX377bdR0tAFwRC9Z4ThcmPr5FxmCFPr0o1AP4WEhBR7uQVBEApD0oEDSP73MLweHgMrGxsci0hFUqYOrXwlLW1ZxsbKCgOqO2PLzSSk9+yPjNBQJO3bZ+5iCcUIe5/Z88u5YIFC9s6dO1kxZA3RaDS4W8Jp+QICAnJZVsPDw2Fra4tKlSrlu05OK60hjG5AX1jDSRAEwVzw5Tviy6/gUK8enFq3Vsv2hKbA28EaNd0lm1dZp1eQsxK0v2l9YV+7NqKXyaCv8kRQUBBmz56t5oIFCtnevXvj6aefxpEjR7KsofyfaWuLGr3AVDp27Ijt27dnW7Zt2za0adNG+cfmt06nTp1KtGyCIAjFReKu3Ug9fRqeDz2U1Zu0+1YyWvo65Nu7JJQNXO2s0a2yE1ZeToRj335I2rsX6bdumbtYQjGRmJioxuNwLligkP3pp59QpUoVtGvXTsVypTWzffv2CAwMxIIFCwq1LZ4EJ06cUBOhDy7/18dmY5f/448/ni1Cwc2bNzFlyhTll8uycKAXw4DpmThxohKuM2fOxIULF9R8x44dmDRpUmEPVRAEodTRabWI+PprODZuDMemTdWy4IQM3EjIRGs/cSuwFAZVd0ZEigb7q7eGtbMzYletNneRhGLi0qVLyjjGuWB+Ct1H5evrqxIgXL58WYlJWmUbNmyIevXqFXrntOT27Nkz6zMFKhk7diwWLVqk3BgMAw7XrFlT7Zu+Kd98841KpMDQYIapcXlyrVixQvnr0p+3du3aWLlypRLbgiAIZZ2E7TuQduECAqZPz7K+0q3AzhpoWkmyeVkKVd3sVHv9ejUFnbp1Q+yaNfB98QVY2UsbCkKpx5GtaEgcWUEQzAFHt1+7bxhs3Fzh//Y7WcvH7QhDUqYW77T9byyAYBn8ezcVnx2LwfJW1vB8dwqqfDkH7gMGmLtYwj0icWTLSWYvQRAEoXiJ37oV6VevwnPU6KxlSRlaHA5PlWgFFkhrXwf4OtlgdZwzHBo2VKG4BEEoXkTICoIglBHf2Mhv58OpVUsVrUDPgTspyND+J4oEy8LG2gr9qjrjj+uJ0PTqj+RDh5B2/bq5iyXcI4yU5OPjo+aC+REhKwiCUAZI+PPP/6yxD47Mtpz+sUEutghwkYemJdK7qjPov7fZuyGs3dxk0Fc5oFmzZirDKeeChQlZZtD64IMPJGGAIAhCcUcq+OYbOLVokc0aq9XpsCc0GS39xBprqbjbW6NzoBNWXkuGU48eiFu/HlqDlOiCIJSikKUZ/fPPP1fJDwRBEITiIWHbdqRfuQqPkdmtseej0xGZqhW3AgtnQDVn3E7S4ESLXtDExqrIFILlcvbsWdSpU0fNBQt0LWDSg927d5dMaQRBECqib+w338CxeXM4GqT+JrtDU+Bsa4UGXhKyyZKp42mPOh52WBVpp+IDx65aae4iCfdAWloarl69quaC+Sm009XAgQNVooIzZ86gdevWcHFxyfb9fffdV5zlEwRBKNck7tqFtMuXEfDhh7m+YzavFj4OsLWWbF6WTv9qzvjmdBxiu/WH4/zZatCXQ82a5i6WIFQ8Ifvcc8+pOfMM54TBu8XtQBAEwTQYxjty/ndwbNQIjg0bZvsuIiUTZ6LT8VIzD6nOckCnQCf8ciEev7vVxcOuripBgv+rr5q7WIJQ8VwLtFptnpOIWEEQBNNJOnAAqWfOwOOBB3J9tzc0BbTD0iIrWD4ONlboGeSMdTdSYNe9J+LWroM2Pd3cxRIEi0fCbwmCIJiJyO++h32dOso/NicMu1XP0w4eDjZmKZtQ/DCmbHy6Fgea9FCDvhL/+kuq2QLhQK+tW7equWChQnbPnj0YOnSoasS6desqv9i///67+EsnCIJQTkk+dgwphw8rayzdsgxJ1+iw/04KWvk5mq18QvET6GKL5j72WB1pB4cGDRCzcpVUswXi7u6O/v37q7lggUJ26dKlKnKBs7MzXn75Zbz44otwcnJC7969sWzZspIppSAIQjm0xtpVqwbnNm1yfceUtMmZOgm7VQ7pW9UZJ6PScbvbICT/8w/SQ0LMXSShkNy5cwfvv/++mgsWKGQ//vhjfPbZZ1i5cqUSshMnTlT/f/rpp/jQyKhbQRAEITupFy8iae9eeAwfDivr3LdhJkHwcbRGdTfJ5lXeaOvnCG8Ha2xwrQMrZ2fErl1r7iIJhYQClsmhRMhaqJC9du2acivICd0LrksOaUEQhAKJWrAANr6+cOnUyWgkg123UtDK1zGXy4Fg+dhYW6m0tRuDU2DVtYca9KXLzDR3sQSh4gjZqlWr4i8jDupcxu8EQRCEvEm/FYr4zVvgMWQIrGxzW1yvxWfgVmImWkta2nJLnyBnpGXqsLdxD2RGRCBxr4wxEYSiUuh+q1deeUW5FJw4cQKdOnVSFoN9+/Zh0aJF+Oqrr4pcEEEQhIpA9KJFsHZxgWvv3ka/330rRYVqalJJwm6VVyo52aCVnwPWRGvRs1YtxK5eDbdePc1dLEGoOAkRAgICMGvWLKxa9d+Iy4YNGyo/2WHDhpVEGQVBEMoFmTExSrS433cfrB2NRyTYHZqMppXslZgVyncork+OxuBmt8GovPgbZNy9Czt/f3MXSzABLy8vPPLII2oumJ8ijSS4//771SQIgiCYTsyv/0V2cR840Oj3sWkaHI9Iw/hGks2rvNPc1wF+Tjb4w70+nrWzQ9z69fCZMMHcxRJMoGbNmiqCk1A2kIQIgiAIpYA2NRUxS5fCtVcv2OQRf3L/7RRodFDdzkL5xsbKCr2DnLA5NA3aTt0Qu3oNdFqtuYslmEBqaiquXLmi5oKFCFlvb29ERkaq/2lK5+e8JkEQBCE3cb/9Dk1cHNwHD86zenaFpqCWux0qOUo2r4pAryBnZGp12NO0FzJCQ5F86JC5iySYwLlz51QyKM4FC3EtmDNnDtzc3NT/X375ZUmXSRAEoVxBS1v0zz/DuX172AUEGF0nQ6vD3tBkDKjuUurlE8yDl6MN2vg5Ym2cHfoEBSF2zVq4dOwozSEIxS1kx44dq+aZ/4t1x9RsHPAlCIIgFEzi7t1Iv3kT3s88k+c6x8NTkZChk7BbFTDT10dHonG92xBUW71QDQi0lUFEglAyPrK2trYqakFaWlphfiYIglChifrpZzjUrw/H+vXzdStgxie6FggVh2Y+9ghwtsEf3o1UMoz4jX+Yu0iCUL4He7Vv3x7Hjx8vmdIIgiCUM1JOn0bKkSNwN5IR0ZBdt5LR0tcB1pLNq0LB9u4d5Iw/72Qgs11HxK5apQStIAglFH7r+eefV0kRbt26hdatW8PFJbs/V7NmzQq7SUEQhHJL9C+/wDYgAM5t2+a5zvX4DNxMyMRDdf8biyBULHoGOWHF5QTsadYXvb97F6mnTsGpeXNzF0vIg1atWsnLhiUL2dGjR6s5s3vpYXYvvkFyrtFoireEgiAIFgqD3Mdv2Qrvxx+HlU3ekQh230qGvTW7mSXsVkXE08EGbf0dsS7ZFn19fRG7Zo0IWUEoKdeC69ev55quXbuWNRcEQRD+PwGClb09XHvmn350561klZJWsnlVXJjp61p8Jq53H4q4PzZBk5hk7iIJeXDx4kV07NhRzQULtMhWr169ZEoiCIJQjtCmpCBm5Uq49eoFa2fnPNeL0WfzaizZvCoyTSr9b9CXR1O8kPozErZugeeDD5q7WIIRkpKS8M8//6i5YKGZvZYsWYLOnTujcuXKuHnzZlZ82d9//724yycIgmCRxG3YCG1CAtwGDcp3vb9D/8vm1cZX3Aoq+qCvPlWdse2uBumt2yFm1WpzF0kQyqeQnT9/PqZMmYJBgwYhNjY2yyfW09NTkiUIgiAwAYJOpwZ5ObdtAzt//wLdCup62qng+ELFpmcVJ2h1wK7m/dWAr9SLl8xdJEEof0J27ty5+PHHH/HWW2/BxmDwQps2bXD69OlCF+Dbb79FzZo14ejoqKIg/P3333mu+8QTT6gBZTmnxo0bZ62zaNEio+tITmRBEEqLpAMHkH7tGtwH5Z2OlqRrdNh3OwWtxRorAPBwsEH7AEesT/OCtacnYteukXoRhJIY7NWyZctcyx0cHArtL7Jy5UpMmjRJiWLGpu3atSsGDhyI4OBgo+t/9dVXuHPnTtYUEhICb29vjBw5Mtt67u7u2dbjRKEsCIJQGkQvWQr7GjXg0KhRvusdDk9FUqZOpSkVBP2gr5uJGlzsfh/ifvsdWklAVOaoUaOGcrHkXLBAIUvr6YkTJ3It37JlCxoVcNPOyezZszFu3DiMHz8eDRs2VK4JVatWVe4LxvDw8FCpcfXTkSNHEBMTgyeffDLberTAGq4n6XQFQSgt0oODkbRnD9wGDlT3ooKSIPg52aC6W6HH3QrllEbe9ghyscUf/s2hjY9Hwrbt5i6SkAMa0B599FE1FyxQyL766qt44YUXlDWVfmD//vsvPv74Y7z55pvqO1NJT0/H0aNH0a9fv2zL+fnAgQMmbWPhwoXo06dPrkgKiYmJallQUBCGDBlSYCYyptyNj4/PNgmCIBQ15Ja1qytcunTJdz3eP/8KSUZrP4cCBa9QceC50LeaM/6K0CG5WWvErl5l7iIJOYiIiMA333yj5oIFCllaP9977z1MmzYNycnJGDNmDL777jvV7f/QQw+ZvJ3IyEg1UMw/x0AIfg4LCyvw93QXoBWY1lxDGjRooPxkN2zYgOXLlyuXAkZYuHz5cp7bmjFjhrL26idahQVBEAqLNikJsWvXwrV3b1g75B+F4HxMOsKSNWgrbgVCDrpXcYKttRX+ajUQyf8eRtr161JHZQi6Nb744otqLlho+K2nn35ahd0KDw9XopONSReBopDTEqHPEFYQFKuMlDB8+PBsyzt06KBM/s2bN1c+t6tWrUK9evXUILW8eOONNxAXF5c1yckpCEJRiNuwAdrkZLj371/gurTGutpZqa5kQTDE1c4anQMdsS7NGzo3d5XpSxCEYhKyvXr1UmG3iI+PD/z8/NT/7I7nd6bC3zLqQU7rK8VxTittTih2f/rpJzz22GOwt8//IWBtbY22bdvma5HlQDUOEDOcBEEQCh1ya+mvcG7bFra+viYJ2RY+DsryJgg56V/NBWEpWpzqPhxx69ZDl54ulSQIxSFkd+/erfxbc8LwVvmFzsoJBSjDbW3fnt2RnZ87deqU72/37NmDK1eumGQF5sOFg9MCAwNNLpsgCEJhST50COlXr6pBXgURmpiBi7EZaOcv0QoE49T2sFPxhTf4NoMmJgYJO3dKVQmCEUweKnvq1Kms/8+dO5fNkkpf161bt6JKlSooDEysQKsqY9Ayb/EPP/ygQm9NmDAhq8s/NDQUixcvzjXIq3379mjSpEmubX7wwQfKvaBu3brKSvz1118rIUvHbEEQhJIi+tdlsKtWDY4Gca3zYuetFNhaQVlkBSEv+ld1xrzTcYhs2hbOK1fBfcAAqawygJubmxqYzrlgQUK2RYsWWckFjLkQODk55euHaozRo0cjKioK06dPV4O3KEw3b96cFYWAy3LGlKUP69q1a9XgMmPQ7eGZZ55RQpsDtxjzdu/evWjXrl2hyiYIgmAqGbdvI/Gvv+A9bpxJPv50K2hcyR7OdkUapiBUEDoGOmHxhQRsaT4Ajy39UIV2s69WzdzFqvDQUPbnn39W+HooK1jp2PduAhzcxVVr1aqlQm75GviA0U2AvrKGmb4sGVpyKYIpmsVfVhCEggif8yWilyxB1e+/h7WTU77rxqRp0HVNCJ5q5K78IAUhP5ZejMeOkGQs3fYRKo98AH5Tp0qFmRn2QjMBlIuLS7nRPZaMyeYAWkmZxUKr1SpXAH7WT/Q/lcYUBKEiok1PR+yqVXDt0aNAEUv2hKZAo4OE3RJMon81ZyRn6LCv6wjErl0ng77KACdPnlTGLs4F81Pofq1ffvkFmzZtyvrMeLIMg8UBWrTaCoIgVCQStm5Vg3FMCblFdgQnoZ6nHbwdxZIjFIyvky3a+DtgnXsDZHLQ144dUm2CcC9C9pNPPlH+sOTgwYOYN28ePvvsMxVOa/LkyYXdnCAIgkXDkFuOzZvDzoTBrimZWuy7kyrRCoRCMbC6C64lA+da90bMihVSe4JwL0KWyQLq1Kmj/v/tt9/w4IMPqsFVzI5VmPBbgiAIlk7KmbNIPXXKZGvs/jspSNPoRMgKhaKJtz2qu9nit7o9/sv0dU0yfQlCkYWsq6urijRAtm3bhj59+qj/mQo2JSWlsJsTBEGwWGKWL4ONry+cWrc2aX0O2qnqaovKLiYHjBEEFQljcA0X7Et1wm3/GohduVJqRRCKKmT79u2L8ePHq+nSpUsYPHiwWn727Fk1GEwQBKEioImNRfwfm+DWpw+sTBi5nKHVYdetFLT1l9ixQuHpEugET3tr/NHhAcSuXw+tGI7MRtOmTVUWUs4FCxSyTCzA5AUREREqnmulSpXU8qNHj+Lhhx8uiTIKgiCUOWKZNlSrVULWFP69m4r4dC06+Bcc2UAQcmJvY4V+1ZyxxaYy4lIzEb95i1SSmbCzs1MhSDkXLCiObEVC4sgKgpAfFLBX+w+AfY0a8J040aTKev9QFHbdSsY33X1NSpogCDmJS9Ngwu5wPB59Ao9GHEfNtWukkszA1atX1eD2OXPmoHbt2tIGZqZIaWWYPWvWrFnKveDpp5/G7NmzVfIAQRCEikDS/v3ICAmBm4mDvDRaHXaEJKF9gKOIWKHIeDjYoGcVZ6yt1Axx5y8ixSB1vFB6UO9s3LhRdI+lCtkjR46oNxC+iURHRyMyMjLrreTYsWMlU0pBEIQyRPSvy2BfsyYc6tc3af0TkWmIStWivb9jiZdNKN/cV9MF8Rpr7GzcCzHLlpu7OIJgeUKW5vT77rsPN27cwLp167B+/Xpcv34dQ4YMwaRJk0qmlIIgCGWE9Fu3kLRnj7LGmuoisC04Gd4O1ioRgiDcCwEutugQ4Ig1dbojZstWlSRBECoyRbLIvvbaa7C1/f/wMfyfGb74nSAIQnkmdsUKWDs7w6VLF5PW5zCE7cFJaOvvCGvxjRWKgeG1XHFb54C9AU0Qt3at1KlQoSm0kHV3d0dwcLDRRAlubm7FVS5BEIQyhzYtDbFr1sK1Rw9YO5rmJnA6Kh13kjXoGCBuBULxUMvDDi18HLCy2WBE/roMusxMqdpSpEqVKmqcEOeCBQrZ0aNHY9y4cVi5cqUSr7du3cKKFSvUwC8JvyUIQnkmfssWFT/WbcAAk3/z580kFf+zobd9iZZNqFiMrOOKGzZu2GPli4Rdu8xdnAqFv78/pkyZouaC+Sl0epkvvvhC+YU9/vjjyPzfWyBjqT333HP49NNPS6KMgiAIZYKYX5fBqUVz2AUGmuxWsDU4SaWktRG3AqEYqe9lj+Y+9ljefAj6LF4M9759pX5LiZiYGOzYsUNlNvXy8pJ6tzSLrL29Pb766ivVkCdOnMDx48dV9AJGLnBwkIw1giCUT1JOn0Hq6dNw62dayC1yJiodt5M0anCOIBQ3I+u44bqDN3beSkXqhQtSwaUEB7iPGjVKzQULErLJycl44YUXlE+In5+fciUIDAxEs2bN4OzsXLKlFARBMDMxv/4KWz8/OLVubfJv/gxOgru9NRqLW4FQAjTwskezSnZY3HQIIhYvkToWKiQmC9n33nsPixYtwuDBg/HQQw9h+/btyp1AEAShvJMZHY34zZvh1q8frGxsTHcruJmMdv4OsLGWTF5CyTCmnjuCnX2w/vhtdZ4KQkXDZB9ZxoxduHChErHk0UcfRefOnaHRaGBj4o1dEATBEmGkAuLau7fJv2G0gtCkTIxr5F6CJRMqOnU87dHBxxZL6vXBiBUrUeV5MTAJFQuTLbKMUNC1a9esz+3atVPxY2/fvl1SZRMEQTA7DG0Us3w5XDp3hk0hQgxuYbQCB2s0riTRCoSS5eGGXohy9MAvuy6pEHFCyeLk5ISWLVuquWBBQpaWVw70MoRCVh+5QBAEoTySuHs3Mu/cgdvAgSb/RqvTYfONJBU7VqIVCCVNFVdb9PKzxq9VO+Pmuo1S4SVMw4YNcezYMTUXLMi1gP5eTzzxRLbIBKmpqZgwYQJcXFyyuSAIgiCUF6KXLIVDg/pwqFXL5N8ci0hDeIoGnQPFYiOUDg839cWBu6GY89dFfDVaCyvrQgclEgSLxOQzfezYsSpagYeHR9ZEP9nKlStnWyYIglBeSLt8GcmHDsFtgOnWWEJrrI+jNep52pVY2QTBEA8HG9wfAGys1AjHN++WyilBGHaURj3OBQuyyP78888lWxJBEIQyRvTiJbDx9oZLhw4m/yZTy2gFSegS6ARrSYIglCKDWwRhxx9X8P5fkfh9sE4lLxKKH/ZQp6enq7lgfqTvQRAEwQiZMTGI27ABbv37w8rW9CSIB8NSEJOmRdfK4lYglC72NtZ4MlCDUw5++PX3f6T6hQqBCFlBEAQjxK5eA2i1cOvTp1D1s/F6EoJcbFHTvdAZwAXhnmndpj7axVzDzIN3EZkoEQyE8o8IWUEQhBzoMjJUJi+Xrl1hUwjf/+RMLXaEJKNLZUfp1hXMAgd5PVbbAdqMTLy77JC0glDuESErCIKQg4Tt25F59y7cBw0qVN3sDElGSqZO3AoEsxLQugUevP0vNl9LwB+nJNZ7ccOwW2fOnJHwW2UEEbKCIAg5iPplMRybNIZ9jRqFdiuo72UHf2dxKxDMB326ezQJQouIy3hr7SmEx6dKcxQjTITQuHFjSYhQRjC7kP32229Rs2ZNODo6onXr1vj777/zXHf37t2quy7ndOHChWzrrV27Fo0aNVLhMThfv359KRyJIAjlgeRjx5F68iTcBw8p1O+iUjXYfycFXSV2rFAGcGnfDg/dOgikpWLq6pPQamWEfXFx8+ZNjB8/Xs2FCi5kV65ciUmTJuGtt95S8diYAnfgwIEIDg7O93cXL17EnTt3sqa6detmfXfw4EGMHj0ajz32GE6ePKnmo0aNwqFD4iskCELBRP38M+wqV4ZT69aFqq4/rieCwY4kCYJQFrC2t4d/h7Z45MwW7L0cifl7rpq7SOWGqKgoLFy4UM2FCi5kZ8+ejXHjxqk3G/qcfPnll6hatSrmz5+f7++YmCEgICBrsrGxyfqO2+jbty/eeOMNNGjQQM179+6tlguCIORH+s2bSNyxA+5DhxY6M9Jv1xLR2s8RbvZm7+gSBIVLp05olByGITaRmLXtIg5dE+EllD/MdsdlMOGjR4+iX79+2Zbz84EDB/L9bcuWLREYGKgE6q5du7J9R4tszm3279+/wG0KgiBE/7IY1u7ucOnWrVCVcTEmHRdiMtCjisSOFcoO1o6OcOnWFX12r0BDXye8sOwYbsemmLtYglA+hGxkZCQ0Gg38/f2zLefnsLAwo7+heP3hhx+UD+y6detQv359JWb37t2btQ5/W5htkrS0NMTHx2ebBEGoeAkQYteuhfuAAbB2cCi0NdbD3hotfQv3O0EoaVw6dYatowOejjmuMs2N/+UIktMzpeKFcoPZ+8ByptBjyre80upRuD799NNo1aoVOnbsqAaKDR48GF988UWRt0lmzJgBDw+PrInuDYIgVCxiV6zgzUJl8ioMTEm74VoiugQ6wtZaUoIKZQu+lLl26wbrndvxescAXItIxJRVMvjrXqBx7PXXX89lNBMqmJD18fFRvq05LaXh4eGFOjk6dOiAy5cvZ32mz2xht0k/2ri4uKwpJCSkUMciCIJlo01OVm4Frj17wsbdvVC/3Ruagug0LXoEOZdY+QThXnDu1Em5GXj+tRGv9KuPbWfDMP2Pc8rIIxSeKlWqKAMY50IFFrL29vYq3Nb27duzLefnTp06mbwdRjugy4EeWmpzbnPbtm35bpNhutzd3bNNgiBUHGLXrIEmIQEew4cX+rdrriSgtocdarrblUjZBKE4Ihi49uqFxJ270Mo+GRO618aiAzfw7W6JZFAUEhISVDhQzgXzY9ao3VOmTFHhsdq0aaMEKP1fGXprwoQJWZbS0NBQLF68WH1m5IEaNWqoQMQcLLZ06VLlL8tJz8SJE9GtWzfMnDkTw4YNw++//44dO3Zg3759ZjtOQRDKLtr0dEQtWKjS0dr6+hbqt3eTM7EnNAXjGsvLr1C2cenQAUn79qnUywNfex2xyRn4/M+LcLG3wROda5q7eBYFe4F79uypBqzT1VGowEKW8V4Zh2369OkqHmyTJk2wefNmVK9eXX3PZYYxZSlep06dqsStPrPGpk2bMMggjSQtrytWrMDbb7+Nd955B7Vr11bxatu3b2+WYxQEoWwT9/vvyIyIgF8RrLEc5GVnY4UukgRBsIBsX/T/jl2+HKkXLuChtvWRkqHB+xvPwdbGGo92+O+5KwiWhpVOnGRywagFHPRFf1lxMxCE8osuMxNXBw6CXZUq8Js6tVC/1ep06P97KOp42OHFZp4lVkZBKC74uI+cOxfWLi4I/ORjQAcs2HcdG07exgf3NcbYToVLyVxROXbsmHKNFIts2cDsUQsEQRDMRfyWrcgICYHH/fcX+reHwlJxKzETvWSQl2AhMHqP28ABSLt4EcmHDqnP47vUxPAWVfDehrP4TrJ/CRaIWV0LBEEQzGmNjfxmHpxatYJD7dqF/v3ySwmo5mqLhl4yyEuwHBzr1oND/foqSgfTMFvb2uGpzjXgYGeNT7dcUL6zrw2on2/IyoqOnZ2diljAuWB+xCIrCEKFJO6PP5B+4yY8R48u9G85yGvnrWT0reYsD3zB4mAK5szISMT/8Yf6TNH6aPvqGNe5prLKvr72NDI1WnMXs8zStGlT3Lp1S80F8yNCVhCECocuIwOR876Bc7t2RbLGrrmSCDtrK3SvLClpBcvDzs8PLh07Inb1GmTGRGctH96yCib3qYs1R2/h2aVHkZKuMWs5BcEURMgKglAhIxVk3LoFz1GjCv1bZvJadTkBXSs7wdlObqGCZeLat6+KZMBwXIb0auCPt4c0xP4rkXj4x38QnZRutjKWVU6fPo2goCA1F8yP3IUFQahQ6NLTEfHtfDh37Aj7GoUfpU2XgvAUDfpVk0xeguVi4+SkwnEl7tqN1IsXsn3Xpro3PhneFDeikvDAt/sRHJVstnKWRTIyMlQYUM4F8yNCVhCECkXM8uXIDAsrkjWWLD4fj0Ze9pLJS7B4nNu2hV3VIER9970a/GhIXX83fDaiGTI0Otz/7X6cuhVrtnIKQn6IkBUEocKgiYtT1lim67SvWrXQvz8blYajEWkYVEOssYLlY2VtDY8HRiD91i3E/7Ex1/eBHk6YOaIZfN0cMPqHf7DrYrhZyikI+SFCVhCECkPk9z9Al5ZWpEgFZMnFePg52aCtv2Oxl00QzIF9lSpw6dwZMStXISM8t1D1cLLDh8OaoGkVD4xfdARrj94ySzkFIS9EyAqCUCGg1SlmyRJ43HcfbL28Cv37yBQNNt9IQv9qzrCRGJtCOcKtXz9YOTsj6ocfoGO6rxw42tngzYEN0buhH15ZfRIL/r6GikzdunWxa9cuNRfMjwhZQRAqBBFzvoS1mxvc77uvSL//9WK8ErC9q4pbgVC+sHZwgOf9w5Fy/DiS9uwxuo6NtRVe7FkHI1sH4aNN5/HFnxdVytuKiJubG3r06KHmgvkRISsIQrkn+fBhxG/apFwKrB0L7xaQlKHFsksJ6FPVGa4Scksohzg2bKSy3EUt/AmZ0f8fW9YQJk54vGMNlQls3q4rmP7HuQopZhmx4I033lBzwfyIkBUEodwnP7jzwQdwqFcXrj17Fmkba68mKjE7uIZLsZdPEMoKqrfCxgaR331n1MVAz/0tg/B8j9r4ef8NvLn+DLTaiiVm7969i08//VTNBfMjQlYQhHJN9JKlSL92Hd7jn1ajtAtLhlaHRefi0DnQCb5ONiVSRkEoC9g4O8NzxANIOXoUiTt35bvuwCaBmNi7Llb8G4y3fqt4YlYoO9iauwCCIAglRcbdu4iYN08FfneoVatI29hyIwl3kjWY2kqssUIFcTFo0wbRP/0Ex8aNYOcfkOe6fRr6q/nXf12GtRXw0fAmyv1AEEoTscgKglBuufvJJ7C2t4fnQw8V6fcarQ7fnYlFa18HVHezK/byCUJZhJE9GMUgfM6XuRIlGBOzL/Wqg18PBeOTzecrpM+sYF5EyAqCUC6J37oVCX9ug9cTT8DGpWjW1D+Dk3E9PhMP1nEt9vIJQlmFAyI9H34Y6VevImbN6gLX79soAM90rYUf/76OuTuvoLxTqVIljBs3Ts0F8yOuBYIglDsyo6IQ9sEHcG7fXgV7LwpanQ7zT8eihY8D6nraF3sZBaEs41CtmoovG7dmLZyaNoVT4yb5rj+0eWUkZ2gwe/sl+Lg6YEz7aiivVK9eHQsWLDB3MYT/IRZZQRDKHWHTp0On0aLSM88U2WdvR0gyrsRliDVWqLC49ugB+9q1lIuBJjamwPVHtQ7CkKaBePu309h2NgzllZSUFJw9e1bNBfMjQlYQhHJF/ObNyqWg0vjxsPHwKLJv7NcnY9Gskj0aeIk1VqiYMMqH58NjAI0G4XPmQKfJ31+WL43ju9ZCp9o+eHn5cZwMiUV55Pz582jSpImaC+ZHhKwgCOWG9JAQ3Hn3XTh37gSXTp2KvJ1NN5JwNS4DD9eTzD1CxcbWzQ1ej4xB6rnziFm+osD1mQFsUp+6qOHjgvG/HEForFgthZJFhKwgCOUCXXo6QidPgbWrG3yeebbI20nX6DD3VCza+olvrCAQh5q14D5gAOLWr0fSv4cKrBQHWxu8NaghGLb5qZ8PIyktf0uuINwLImQFQSgX0I8v9cIF+E6aBOsiRikga68mIDQxEw+JNVYQsnDp3h2OTZog4quvkRZ8s8Ca8XS2xzuDGyE4OhmvrD4pCROEEkOErCAIFk/Cjh2I/vlneD3yCBzq1CnydhLTtZh3MhZdKztJ3FhByOH/6jl6NGy9vRE+41NoEuILrJ/qlVwwpW89bD0TVq7CcrEu7O3tJflDGUGErCAIFk3qpUsInfYanDt0gPuQIfe0rQXn4pCYocUYscYKQi6sHRzgNXYstMnJCP9iVoHJEkiHWpUwpl01fLnjEnZdCC8XtdqyZUukpaWpuWB+RMgKgmCxZMbE4NZzz8POzw8+L7xwTxaSO0mZWHQ+HkNqusDHyaZYyykI5QVbLy94PfqocuOJWvAjdCg4k9fotlXRpoYXJq08gZDo5FIpp1BxECErCIJFosvIQOikydAkJMB32jRYOznd0/ZmH4+Bo40VhteULF6CkB8OtWrBc8QIJGzfgbh16wusLGsrK0zpWx/O9jZ4dslRpGZoLLqCGXarVatWEn6rjCBCVhAEi0On1eL2m28h+ehR+E2dqiyy98Lhu6n440aScilwtpPboiAUhHPr1nDt2wcxv/6KxH37Clzf1cEWrw9ogMvhCfho0zmLrmAmQjh+/LgkRCgjyB1bEASLI2L2bMRv3Ajfl16CY6NG97StTK0OHx2OQl0PO/QMujerriBUJNz69IVT61aInDsXKWdOF7h+LV9XPN21Fpb+E4yNJ2+XShmF8o8IWUEQLIroX35B1IKF8H7ySbh07nzP21t5OQGXYzMwrpG76gIVBKEQkQxGPPhfGttPZyLt+rUCfzOgcQC61fXB6+tO4UZkklS1YPlC9ttvv0XNmjXh6OiI1q1b4++//85z3XXr1qFv377w9fWFu7s7OnbsiD///DPbOosWLVIXV84pNTW1FI5GEISSJGbVKtyd8Snchw2D++DB97y9u8mZmHM8Bn2qOqOOp6SiFYTCYmVrC69HH4ONjw/CPvwIGWF38l/fygov9KwDd0c7vLT8ONIztVLpguUK2ZUrV2LSpEl46623lL9J165dMXDgQAQHBxtdf+/evUrIbt68GUePHkXPnj0xdOhQ9VtDKHLv3LmTbaJQFgTBcold/xvC3nsfbgMHqlHTxQFdCuxsrPBIfUlFKwj3EpaLPSTW9vYIe/8DZEZF5ru+s70tXu1XH+fuxGPWtosWV/E0vq1atUrNBfNjpdPpCo6dUUK0b99ejfybP39+1rKGDRti+PDhmDFjhknbaNy4MUaPHo133303yyJLcRwbG1vkcsXHx8PDwwNxcXFKFAuCYF7iNv6B26+9BtdevVDpmWdgxdyX98j24CS8vDcCU1p4olOg+MYKQnGEw4v87jvYODkh8KMPYePhme/6647dws8HbmDxU+3QrZ6vNIBgWRbZ9PR0ZVXt169ftuX8fODAAZO2odVqkZCQAG9v72zLExMTUb16dQQFBWHIkCG5LLaCIFgOsevW4/a0aXDt3r3YRGxsmgbT/41Ca18HdAyQ3hpBKK4Ysz5PPw1NYiLCpn+o5vkxvGUVtKzmiVdWnURUYprFNMLdu3cxe/ZsNRcqsJCNjIyERqOBv79/tuX8HBYWZtI2Zs2ahaSkJIwaNSprWYMGDZRVdsOGDVi+fLlyKejcuTMuX76c53aYoYNWWMNJEATzE7N8Oe68+SZc+/RBpeeeKxYRq3cpSNHo8GwTD0kzKQjFiK2PDyqNH4/MyEiEffghNMl5D+ji4MpJveshXaPFq2tOwYwdxIUiNDQUr7zyipoL5sfsg71yZuLhiWxKdh6K1Pfff1/52foZxJDs0KEDHn30UTRv3lz53NKPpV69epg7d26e26IbA10J9FPVqlXv8agEQbhXohb+hLAPpsNt0KBis8SSbcFJ2HQjGeMausPbUTJ4CUJxYxcQAO9x45B5+/Z/ltmUvLN5ebvY4+VedbHzQjiW/nNTGkOwHCHr4+MDGxubXNbX8PDwXFbanFC8jhs3TonUPn365LuutbU12rZtm69F9o033lD+sPopJCSkkEcjCEJxwZfZ8FmzEP755/B44AE1iOReUs8aEp6cifcORaGdvwO6Vha/WEEoKeyrVIE3LbOhtwoUs+1qemNw00B8tOk8Lt1NkEYRLEPI2tvbq3Bb27dvz7acnzt16pSvJfaJJ57AsmXLMNiE8Dt8KJ44cQKBgYF5ruPg4KAGdRlOgiCUPjqNBmHvvoeoHxfAa+xYeI0ZU2wiVqvT4fUDkeqmN6GJp7gUCEIJYx8UBO9x45EZEqJ6VzRJefvMPtm5BvzdHfHy8uMWn8JWqECuBVOmTMGCBQvw008/qZzFkydPVqG3JkyYkGUpffzxx7OJWH6mbyxdCGjN5UQrqp4PPvhAxZa9du2aErC03HKu36YgCGUTbWoqbr08EbFr18LnxRfhMXRosW7/53Px+CcsFS8284S7vdm9qgShQmBftSq8n34amaGhKjSXJtG4xdXB1gZT+9XD1YhEfLa1bIfkogsiQ39yLpgfs97NGTbryy+/xPTp09GiRQsVJ5YxYhlxgDD+q2FM2e+//x6ZmZl44YUXlIVVP02cODFrHYbdeuaZZ1QYL0ZAoDM2t9uuXTuzHKMgCAWjiY1F8JNPIWnfPvgxzFaPHsVabcfCUzHnRAzuq+mCZj4O0iSCUMqW2UrPPovM8HAVC1oTZzw8Zk0fV4ztWAM/7b+O3RfDy2wb1a5dWw0o51yo4HFkyyoSR1YQSo+M0FAEP/2MGuXs//rrcKhXr1i3H5WqwQObbqOSozXea1cJttaShlYQzEHG3buI+vFH2Li6IuC991SEA2MuQNM3nsONqCRsndQNvm5l78UzIyNDGc08PT1hZ2dn7uJUeKR/TRAEs5F6/jyuj34I2qQkBH70UbGL2EytDlP3RSBdo8PkFl4iYgXBjNj5+8NnwgTlRnT7rbeMprNlSK6JfepCo9PhlVUnoNWWPVvb6dOnVbQkzgXzI0JWEASzkLhvP2488ihsPD0Q+PHHsKtcudj38cWxGPx7NxWTWnhKqC1BKCtxZhkT2sYGt996G2nXr+Vax8vZHpN718Pey5FYsC/394JgiAhZQRBKndi16xDy7LNwbNAAAe9/ABvP/FNZFoX1VxPwy4V4PNnQHU0qlb3uSUGoqNh6eCifWRs3V4S98y5SzpzJtU6r6l54oGUVNfDrREjRU84L5R8RsoIglBp0yY/4ei7uvPUWXHv1UgO7rB2LP0UsrbCMF9unqhMGVHMu9u0LgnBv0E+20jPPwjYoCHc//BBJBw/mWufRDtVR29cVL/x6FHHJGVLlglFEyAqCUCro0tNx+/U3EPntt/B65JH/snXZFH9mratx6Xhx91008LLHuEaSglYQyirWDg6o9OSTcGzSRCVBiftjY7bv7WysMa1/fcSnZGLqmpMWk8JWKF0kaoERJGqBIBQvmrg43HrpJSQfPwGfF16Aa5cuJVLFd5MzMebPO2pQ14ftK8HFTt7VBaGsQ4Eav2ULknbvhvvQIfAeOxZWVv9/7f57PQofbjqPNwY2wLPdzR/ySqPRICkpCS4uLipDqWBebM28f0EQyjnpISEIeeZZFV4r4N134diwYYnsJyZVg6d23EW6Bni3rbeIWEGwEJi9z2PQIOUrH79hAzLvhsN34sQst6N2NSvhwVZBmLn1AppW8UCnOrnDdpUmFK+SAbTsIOYKQRBKjORjx3Bj1CgVbifwk09KTMQmpGvx9M67Kmbsu+284eMkVhJBsDRcO3VS1tiUU6dw5+23kRkVmc1ftlmQJ15YdgyhsSlmLefly5fRv39/NRfMjwhZQRBKhLiNGxE89gnYBVZWItYuMLBE9hOfrsG4v8JwMyEDb7f1RmUX6WgSBEuFL7s+zz2HzNhYhL72OtIuX1LLbaytMLVffdjbWGP8L4eRlJZptjImJCRg27Ztai6YHxGygiAUKzqtFuFzvsTtV6fBpUsX+L/zDmzc3EqklmPTNBi34y6ux2fg3baVUNNdsuwIgqXDl16fF19UYbruvP0OEnbtUss9nOzw1uBGuBGZjMkry2ayBKH0ESErCEKxoUlMwq0XXkTUDz/A69FHUen552FVQikcObDr0W1hCE7IVCK2loeIWEEoL9i6uanIJk6tWiFy3jxELVwIXWYmavq4YGq/eth+7i4+3XrB3MUUygDSBycIQrGQdv26ikyQcfsO/N54A86tWpVYzV6LS1c+sRzY9WGHSqjiKrcyQShvWNnawmPECNhWqYz4jX8g7epV+L0yBe1q+mB811r4Ye81+Lk5qP+FiotYZAVBuGcS/voLNx4cCV3Kf4O6SlLEHriTgtFb78DGykpErCBUgIgGrh07odKECci4exehU6ci+eQJ3Ne8Mka0CsJHm87jt+OhpVqmqlWrYt68eWoumB+JI2sEiSMrCKahy8hAxNdfI+rHBXDu0AE+zz8Pa2fnEos1ueRiAj47Go1mPg6Y1NxTQmwJQgVCk5SE2BUrkHbpEjyGD4Pn6Ifw9Z7r2HUxHHMfboXBzUpmQKlQtpH+OEEQikT6rVsInfIKUs+ehddjj8H9vvuU9aQkSEzX4u1/IvFncDKG1nDBo/Xd1ChmQRAqDjYuLvB+6ikk7dmDuA0bkXL6DCZMmgSNToeXlx8HbwkDm5a8mI2OjsbmzZsxaNAgeHt7l/j+hPwRi6wRxCIrCAVk4dm4EWHTp8PaxVUFLneoV6/EquxoeCqm7Y9ATJoWzzfxQMdAJ2keQajgpIcEI2b5CmgTE+Hx+ONYgOrYdyUSnz7QDKPalmyX/7Fjx9C6dWscPXoUrUrQjUowDbHICoJgMhl3wxH2/vtI3LULLt26otK48bB2cSmRGkzO1OKbk7FYdD4e9bzs8FYbb/g7yy1LEATAvmo19RIdv2kTYn/4AY+0aA6ntvdj2tpTiEhMw/M9apdYD5FQtpCngiAIBaLTaBC7ahXCZ89RI4n9pk2Dc7t2JWbx3XUrBR8fiUJkigYP13PDfbVc1OAuQRAEPdYODvB84AE4Nm6M2LVrMfjCp3Ae+CQ+//MiroQnYsYDTeFoJ1n+yjsiZAVByJfkI0cQ9uFHSLt4Ea69eil/2JJKcHA2Kg2fH4vBobupaO5jjzdbeyNAMnUJgpAPjvXrw2/yZMRv3Yqev30H76ZdseQkcDUiEfMeboVqlUpmAKpQNhAhKwiCUTiQImLuXCTt3QuHunUROGOGmpcE56LT8O2pWPx1KwVBLrZ4o7UXWvk6SNegIAgmYe3kBM/774dTy5ZovX49vG5cxOLWIzDoq734+IGmKlxXcbkauLi4oEOHDmoumB8Z7GUEGewlVFTYrZ904ACiFy9RI4PtqlSBx8iRcOnUCVbWxRt2WqPV4e/bKcoHlhbYAGcbPFjHFV0DnSQigSAI95QmO+mffxDx126sqNEFR33qok+9SvhoRAsEeDhKzZYzRMgaQYSsUBEHccVv2YzYlauQfv067GvUUOG0XDp3hpVN8fqY3U7KxIZriVh9JQG3kzSo42GH+2q6oL2/owhYQRCKNe4sk7UcuByBNXV7INPBCS/0qotx3euI72w5QoSsEUTIChXB8pp+9SoS9+1D4s6dSD58BLC1hXObNnAfOBAODRsWa7d+XJoG24KTselGEv69mwp7Gyt0CnBEv2rOSsjK6GJBEEqKzOhohG3bifUJzthXpRl87YGJAxrjwQ41YWdT+J4mCb9VthAhawQRskJ5Q5OYhLTLl5B6+jRSTpxE8tGjyLx7F1Z2dmrEr3OnTnBu314FHC8uolM12HkrGX/eTMLBsFTodEDjSvboWtkJHQMc4WQrGbIFQSg9MiMicGn3QaxP9sBxv3oIsMnAuG618VD3BnBztDN5OyJkyxYy2EsQygnatDRkhN5Gxq0QpAeHIP3mTaTfuIG0K1eQeeeOWofC1b52bRU6y7FZMzg2aqRC2BQXoYkZ2BGSjL9CknE0PE0ta+htjycauqNDgCO8HCQUjiAI5sHW1xeNRt6HelFRuLD/KP6IcMCnf1ljzs5ruK+2Gx4d3ApNqnhK81gYImQFwYLcATQxMUqcpt8MRkZIsBKsGbduqXSxmsjI/1/Z1hZ2AQGw9feHc9u2sK9WDXbVqsG+alUlZouzTJdi/xOvO0KScCEmA3bWQJNKDnimiQfa+jnAQ8SrIAhlCNtKldDkvn5omJKC0H+OYVtwErak18SKa/tRxy4DI9pXx7AuDVDZU7IIWgLiWmAEcS0QzI0mMRFp588j9cJFZVFNu3QJaVevQhsfn7WOjbe3Eqq0Mqi5nx/s/PzUnN8V9yAtPVqdDqci07A9JBnbg5MRkpgJZ1srFS6rnb8jWvo6iNuAIAgWA1/Ik69cxeETV7E/1RmnvWog08YWTZ0yMbhtDQxqXzdbLFpxLShbiJA1gghZoTTRpqcj9ezZ//xXT51GyulTyLgZ/N+XdrawDwqCXVBV2HFeubKabAMCitUloCAytTocDU9VA7YoYCNSNPC0t0YbfwcVbYAWWDtrybwlCILlu2jFnDqLQ5fv4rDGHRe8qiLDxg61bNPRt4Ev+nVpiAZ+zrhzOxRBQUFwdJRwXuZGhKwRRMgKJUlmVBSSjx1DyrHjSDl+DKlnz0GXkfGf/2qtWsqH1YHzmjVVHFemhDUH6RodDoalKKvrX7eSEZumhY+jDdpTvAY4or6XvaSNFQSh3KJJSUHcmfM4ej0Cx1MccNazKpLsneGuS0dnb2v0aF4NPTs1hJ+7iFlzIkLWCCJkhWINc3XjhhKtyceOIvnIUWTcvKm+s/HxUakVHTjVq6dit5pLtOqJSdNg3+0U7AxJVskKkjJ1CHS2US4DHKwlobIEQaiI6DQapFy/gQtXbmPPzbvYfmQX7HqNh51nAOpo4tHR2xpdGldB567N4Orlbu7iVijMLmS//fZbfP7557hz5w4aN26ML7/8El27ds1z/T179mDKlCk4e/YsKleujGnTpmHChAnZ1lm7di3eeecdXL16FbVr18bHH3+M+++/3+QyiZAV7sW3VbkInDyFlJMnlcVVExsHWFnBvkZ1ONSrD4cGDeDYoIHybTU3zK51Njod+++kKOF6MiINWgC1PezUQC0K2KquthLnVRAE4X9cuHQejz3zKL58fQYS4YFTiVa4YOuFOAdX2GozUT8lAm2c0tEmyB1tm1RHpSb11f1e4mWXDGY1/6xcuRKTJk1SYrZz5874/vvvMXDgQJw7dw7VqlXLtf7169cxaNAgPP3001i6dCn279+P559/Hr6+vhgxYoRa5+DBgxg9ejQ+/PBDJV7Xr1+PUaNGYd++fWjfvr0ZjlIoj/D9LzMsTA3CSuVArPMXkHL2bJa11crZGQ516sC1b7//rK5168K6DOTlTs3U4nxMOo6Fpymf18PhqUjM0KnBWk0rOeDpJh5o7esAb0cJkyUIgpAflWpVRed6DdGfvrUaDW4Gh+NkSALOwQFrdd745bYDrEJjUX3172iUeBuNHTPRxN8Z9Wv4w6VWDdhXr67GPpTmeIfyiFktshSWrVq1wvz587OWNWzYEMOHD8eMGTNyrf/aa69hw4YNOH/+fNYyWmNPnjypBCyhiKVFdcuWLVnrDBgwAF5eXli+fLlJ5RKL7P/D00OXnAxNfDw08QnQJsRDk5AAbWLif1NSkrJC6lJSoE1OgTY1Fbq0NDVxEBMyM1WXDLS08/0PZoyytYGVja3yC7VysIe1vQOsnBxh7egEa2cnJfqsnZ1h7eIKa1dOLrDh3M1Nfeb/Vk5OJfaGy+NmhAAG0M4MD0d6aCgyOIXcUilcGaOVx64Ox8lJ3ZDo36p8W+vUUQOyrKzNF/CfkQXCkzW4Gp+Ba3EZuBCTjvPR6bgUmw6NDnCwsUI9Tzs08rJHUx8H5TJgK4O1BEEQTLbILvlhKRrUa5jnMyQ0MRPnQmNwMTwJV5OBUK0DdFZWsNVqUC0+DNUTON1FNZsMVPewRzU/V3hUDoBdQCBsA/xh97+oNCoKjRmfJ2Uds1lk09PTcfToUbz++uvZlvfr1w8HDhww+huKVX5vSP/+/bFw4UJkZGTAzs5OrTN58uRc69BlIS/S0tLUpCcuLi5L0FoS6p2EwjE9XYlIznWpqf/9T6HJ//WCMzkZuqQkaJOT/idGk5QwVSI1IeE/cUrRmpgEUIgaw9r6P7HJUZuOjrCypyC1B+ztYUWhamv3XwgoB4dsF6GOolarhS4zE6DwjclUg510GSxzhiqzKvf/RHGecP//E7wUk6osTo7//c992tkrkazKQZFm/T8rI8W5VvPf/jMyoE1NU/Wh/Z9g1/5vYpmysLKCTSVv2Pj4qvisdo0bw4HRA4KCYOvjkyWoKddTOSUno7hdANI0OiRnapGcqUNihhYJ6VrEpWvVICxm0YpM1SA8ORN3kjS4naxRg7UI47pWcbFVLgIP17BDbd4wXW1gkyVc05Gakl6s5RUEQSivpKQkZ80TkxLzXM/TGuhU1VlN+h6x4AQNbiak42aMD27Gu+NgWn0k6/73fMwA3K4ko9LpMPgkX4RXWgI80hLhnpECD3sreDjawt3ZHm5uLnB2c4Srmyuc3Zzh4OoKW3fX/56H/3sWWjnQMOQAKz6b7exh7WCv4otbmnuDm5tbgWU2m5CNjIyERqOBv79/tuX8HBYWZvQ3XG5s/czMTLW9wMDAPNfJa5uE1t8PPvgg1/KqVasW8qiEcs1lWCzXzF0AQRCEcsYzE58xdxHKPTQsuru7l+3MXjmVNq2K+alvY+vnXF7Ybb7xxhtqAJkerVaL6OhoVKpUyeLeXvKC1mUK85CQkAJPCsEykDYtf0iblj+kTcsf0qala5EtCLMJWR8fH9jY2OSylIaHh+eyqOoJCAgwur6tra0Snfmtk9c2iYODg5oM8fQsn/mWKWJFyJYvpE3LH9Km5Q9p0/KHtGnZwGzew/b29mjdujW2b9+ebTk/d+rUyehvOnbsmGv9bdu2oU2bNso/Nr918tqmIAiCIAiCYJmY1bWA3fmPPfaYEqIUoD/88AOCg4Oz4sKyyz80NBSLFy9Wn7l83rx56ncMwcWBXRzoZRiNYOLEiejWrRtmzpyJYcOG4ffff8eOHTtU+C1BEARBEASh/GBWIctQWVFRUZg+fbpKiNCkSRNs3rwZ1atXV99zGYWtnpo1a6rvGZXgm2++UQkRvv7666wYsoSW1xUrVuDtt99WSRGYEIHxait6DFm6Trz33nu5XCgEy0XatPwhbVr+kDYtf0ibli3MntlLEARBEARBEIqCRNgVBEEQBEEQLBIRsoIgCIIgCIJFIkJWEARBEARBsEhEyAqCIAiCIAgWiQjZcszHH3+sojg4OzubnODhiSeeUNnMDKcOHTqUeFmFkmlPjuV8//33VYQPJycn9OjRA2fPnpXqLiPExMSoEIQeHh5q4v+xsbH5/kau0bLHt99+q6LqODo6qvjof//9d77r79mzR63H9WvVqoXvvvuu1MoqFH+b7t69O9dzk9OFCxekuksBEbLlmPT0dIwcORLPPfdcoX43YMAAFfpMPzHkmWCZ7fnZZ59h9uzZKv7y4cOHVea7vn37IiEhoUTLKpjGmDFjcOLECWzdulVN/J9itiDkGi07MLzjpEmT8NZbb+H48ePo2rUrBg4cmC10pCHXr1/HoEGD1Hpc/80338TLL7+MtWvXlnrZheJpUz0XL17M9uysW7euVHFpwPBbQvnm559/1nl4eJi07tixY3XDhg0r8TIJJd+eWq1WFxAQoPv000+zlqWmpqrffvfdd9IEZubcuXMMfaj7559/spYdPHhQLbtw4UKev5NrtGzRrl073YQJE7Ita9Cgge711183uv60adPU94Y8++yzug4dOpRoOYWSa9Ndu3ap6zYmJkaq2QyIRVYw2k3i5+eHevXqqQxq4eHhUksWCC0/YWFh6NevX7ZA3t27d8eBAwfMWjYBKjMh3QkMk7XQjYfLCmofuUbLTi/J0aNHs11jhJ/zakO2e871+/fvjyNHjiAjI6NEyyuUTJvqadmyJQIDA9G7d2/s2rVLqruUECErZIPdJ7/++it27tyJWbNmqe7oXr16IS0tTWrKwqCIJf7+/tmW87P+O8F8sA34wpgTLsuvfeQaLTtERkZCo9EU6hrjcmPrZ2Zmqu0JltemFK8//PCDcg9Zt24d6tevr8Ts3r17S6nUFRsRshYGB+4Ycyo3nPhmfy9pgwcPHqzSBQ8dOhRbtmzBpUuXsGnTpmI9DqF02pNwGzkHgOVcJpinTY21Q0HtI9do2aOw15ix9Y0tFyyjTSlc2XvZqlUrdOzYUQ0U43P0iy++KKXSVmxszV0AoXC8+OKLeOihh/Jdp0aNGsVWrXzTrF69Oi5fvlxs2xRKpz05sIvQisB21ENXkZzWBqH02/TUqVO4e/duru8iIiIK1T5yjZoPHx8f2NjY5LLU5XeN8bo0tr6trS0qVapUouUVSqZNjUE3oaVLl0qVlwIiZC3wIuNUWkRFRSEkJCSbEBIsoz0ZOoYPze3btyvfLb3/F0P/zJw5s0T2KZjeprTcxMXF4d9//0W7du3UskOHDqllDLNmKnKNmg97e3sVmonX2P3335+1nJ+HDRuWZ7tv3Lgx27Jt27ahTZs2sLOzK/EyC8XfpsZgtAN5bpYS5hhhJpQON2/e1B0/flz3wQcf6FxdXdX/nBISErLWqV+/vm7dunXqfy5/5ZVXdAcOHNBdv35djcTs2LGjrkqVKrr4+HhpNgtrT8KIBYxSwGWnT5/WPfzww7rAwEBpzzLCgAEDdM2aNVPRCjg1bdpUN2TIkGzryDVatlmxYoXOzs5Ot3DhQhWJYtKkSToXFxfdjRs31Pcc6f7YY49lrX/t2jWds7OzbvLkyWp9/o6/X7NmjRmPQriXNp0zZ45u/fr1ukuXLunOnDmjvqe8Wrt2rVRsKSBCthzDMD28mHJOFKh6+JnhnEhycrKuX79+Ol9fX3URV6tWTW0jODjYjEchFLU99SG43nvvPRWGy8HBQdetWzclaIWyQVRUlO6RRx7Rubm5qYn/5wzhI9do2eebb77RVa9eXWdvb69r1aqVbs+ePdmu2+7du2dbf/fu3bqWLVuq9WvUqKGbP3++GUotFFebzpw5U1e7dm2do6OjzsvLS9elSxfdpk2bpIJLCSv+KS3rryAIgiAIgiAUFxK1QBAEQRAEQbBIRMgKgiAIgiAIFokIWUEQBEEQBMEiESErCIIgCIIgWCQiZAVBEARBEASLRISsIAiCIAiCYJGIkBUEQRAEQRAsEhGygiAIgiAIgkUiQlYQBEEQBEGwSETICoIgWDA9evTApEmTUBExx7FHRUXBz88PN27cgCXy4IMPYvbs2eYuhiAUGyJkBcGCCQsLw8SJE1GnTh04OjrC398fXbp0wXfffYfk5GRUVEpT4FRkIWlu1q1bhw8//LBU9zljxgwMHToUNWrUgCXy7rvv4uOPP0Z8fLy5iyIIxYJt8WxGEITS5tq1a+jcuTM8PT3xySefoGnTpsjMzMSlS5fw008/oXLlyrjvvvty/S49PR329vbSYFIXFo+3t3ep7i8lJQULFy7E5s2bYak0a9ZMifBff/0Vzz33nLmLIwj3jFhkBcFCef7552Fra4sjR45g1KhRaNiwoRKzI0aMwKZNm5TVSG8xfPHFFzFlyhT4+Pigb9++anlaWhpefvll1U1Kay4tuYcPH87aPh92X375ZbZ9tmjRAu+//37WZ/22OVFQV6pUCW+//TZ0Ol2+ZddqtZg5c6ayJDs4OKBatWrKSmRKufT75TrTpk1TYiYgICCrXE888QT27NmDr776ClZWVmrSdwPnVRdbt25V+9Efw5AhQ3D16tWs/a1Zs0bVrZOTk/q+T58+SEpKyndfObnX4zKVgvaTkJCARx55BC4uLggMDMScOXNMsipfuXJFHR/Prd69e8PZ2Rn169fHoUOHClU+U4/T1PrSlzuvNtLDc/Kzzz5DrVq11DrNmzdXvykMW7ZsUddcx44ds5bt27cPdnZ2qrx6rl+/rurq5s2bKIvwBXf58uXmLoYgFAsiZAXBAqGf3rZt2/DCCy8oQWIMPkj1/PLLL+oBvH//fnz//fdqGUXE2rVr1XfHjh1TorJ///6Ijo4uVFn026ag+frrr5UwWrBgQb6/eeONN5SQfeedd3Du3DksW7ZMuUUUplz8nsfO/VKgTJ8+Hdu3b1eikkLj6aefxp07d9RUtWrVfOuCgofilkLpr7/+grW1Ne6//34luPn7hx9+GE899RTOnz+P3bt344EHHlDCqKB9GXKvx2UqBe2Hx8lj37Bhg9ru33//rdYriJMnT6pzatasWeplhZ/5AvL666+bXLbCHGdhzs/82kgPy/zzzz9j/vz5OHv2LCZPnoxHH31UvYiYyt69e9GmTZtsy06cOKFeIvlCZriML0XVq1dHWaRdu3b4999/s4lvQbBYdIIgWBz//PMPn9C6devWZVteqVIlnYuLi5qmTZumlnXv3l3XokWLbOslJibq7OzsdL/++mvWsvT0dF3lypV1n332mfpcvXp13Zw5c7L9rnnz5rr33nsv6zO33bBhQ51Wq81a9tprr6lleREfH69zcHDQ/fjjj7m+M6Vc+v126dIl22/btm2r9q3/fuLEibm2b6wujBEeHq7q9/Tp07qjR4+q/2/cuGF03bz2VRLHVdD+C9oP657fr169Ouv72NhYnbOzc4HH8M477+g8PT1V3eiZN2+ernHjxur/4cOHq+9HjBiR73ZMOc7C1BfLXVAbcXuOjo66AwcOZFs+btw43cMPP6wzlWHDhumeeuqpbMvGjx+ve/zxx7Mte/fdd1XZisrVq1d1GzZsyPq8fft23ezZs4tlW+TkyZP51pcgWBJikRUEC8bQ6kpoZaE1qHHjxtmsLTmtSOw2z8jIUD62etg9SksNLVqFoUOHDtnKQQvl5cuXodFojK7P7bNs7J7OSWHKRV8/Q9hNHh4eXmB5c9aFfr9jxoxR3c7u7u6oWbOmWh4cHKy6oFlWdluPHDkSP/74I2JiYgrcT2kflyn7oV81v+dnPR4eHspFoCBogaW7iq+vb9Yybo+WUkI3gMWLF8NU8jvOwp6fBbURrf6pqanKlcTV1TVrYnkNXUhM8ZGlm4MhvN7ocmPI8ePHVZmKCl0YLly4kPWZbhK0IBfHtghdK0hFHhAqlB9EyAqCBULxQPGY8wFFIcbv9A8qPTndD/RdrjmFMJfrl7F7PaevK8XFvZKzbIUtl6GwMYTf0xWgIIy5YlCg0V2DAohd3Xq/Tw6Ms7GxUV3eFASNGjXC3LlzlfCjH6SplMZxmbKf/L43Rcga+obqBZtexPXs2RNubm4wlfyOszD1RQpqI/126d9L4amfKHAL4ydLv2pDgcyXNboptGzZMtt6dIUwFLfDhw/H6NGj0bZtW9StW1fVpZ4lS5agffv2SoTTd5WuDnSD4LnI7VI8Dxw4MEvAX7x4EYMGDULr1q2Vj3BkZKRaznXee+899WJJlwYem7FtEb17huFLiSBYKiJkBcEC4WAWWpfmzZuXbUCLqVDsMnIBB6oYilQOHKO/n/4hR99DPQzXY0y8/fPPP7k+82FNcWEMfkcxS1/UopTLFLiNvCzCOaGApUjgA59WPe4np8WV4onWwQ8++ECJN25//fr1Ju+ruI6rIAraT+3atZWApOXesF1pQc+PuLg4NXApp2AzZo0sjeMwRn5tRHFLH1Za2Lltwykvn2Zj8PgpEPVQVFIcMkKInoMHDyI0NDSbRfbUqVOqJ4A+2DzPDAdRUpTyxen06dNKKDPySJMmTdT1wePgtcL24XXDngz6xf/www84evSoigmr90c/c+aM6kng9Uef7Y0bN6J79+65tqVfNygoSO1PECwdCb8lCBbKt99+qx7cfEByxDe7amlF5cOSllpabPKzSjL0zquvvqpGjXPQDgfcsKtx3Lhxap1evXph0aJFylrp5eWlBmYZE6chISFqANGzzz6rLFG0hnFAUF6wa/a1115Tg3koNngMERERyrLFfRdULlNgxAWKA0YQYBcyt8W6MQaPjS8GFAfs3qbYMRzAxO1QCPTr10+NoOdnllcvqEzZlyn1XRwUtB9aTMeOHZv1PY+HVjyW15ilUw8tiGx7Q3FGYUvBXxJCtrD1VVAb8binTp2quudpnWUEBAr4AwcOqDZjnZgCB5txoCKPm+cNhTzhOU/XCkZ24JzoXXsSExPV/7xGCMukF9i0MPO8YzxcWv957vH4bt26lSWw+RLBMnKAIge/UUgzqoZ+H+PHj1fr8AWFUTQIryu6jBDDbenhAD/WlSCUB0TICoKFQusarSyMIcuHKx9YtDrR+sSHNsNz5cenn36qHuqPPfaYCslEQfznn3+qBzThNukDyYcmH4oMPG/MIvv4448rqxT9Fyl2XnrpJTzzzDP57puimA9mBme/ffu2EpATJkwwqVymwOOnOGFdsGwsd14B7CniVqxYoQQIrVfskmb0BXbbEvrMcrQ6rWgUP+y2pVBnV25h9lUcx2UKBe2HWZ1Y12xXHhtfKPgyktP3M6eQbdCgQTa3EJ57HJlfUokBClNfBbUR4flLkcuEBjyvWfZWrVrhzTffNLlM7P5nOVatWqVe3Chk2TPCNue5w3OA5Wb0hG+++Ua5YtAaS591/UsgX/a4HcIXRYpflp11y3JTtFepUiVrn7Se8veEVlseFyM0GMIoFIZ+z/wNr0HeEwy3RegrTCHNuhSEcoG5R5sJgmC5mDJiXyjbcES/h4eHbsGCBfe8rV27dpkUtcCS2bRpk4rKodFodP369dO9/vrr+a4/f/58Xb169VTUBUZ8aNasWVa0gFdeeSWr3r/88ktdUFCQbu/evbqRI0dm/f67777TffLJJ+r/uXPn6saOHZv13alTp7LW+eijj7KWMzJHTExMrm3pI0307du3WOpCEMoC4iMrCIJQgaAllcHwOVqf1kEmRyDDhg27p+2y250RA5j1iv6XOZMXlBfo00prLP1gaanOGX0hJ7TIMiYxXWjorvP5559nxZeltZmWYvqy0lebllpadukTy//pIkSXGy4jTz75JGJjY5V1nG4ejL9MuI7eyksfW7oz0OKcc1uELgh0hRCE8oIV1ay5CyEIgmXC7nf6SObMACaUbSFLv0oOVKIvJX2p6W6gF0KCaYSFhSmXGMOuf2PQH5cvDoUZVCYIgumIkBUEQRCEEoKRBAoTqk0QhMIhQlYQBEEQBEGwSMRHVhAEQRAEQbBIRMgKgiAIgiAIFokIWUEQBEEQBMEiESErCIIgCIIgWCQiZAVBEARBEASLRISsIAiCIAiCYJGIkBUEQRAEQRAsEhGygiAIgiAIgkUiQlYQBEEQBEGwSETICoIgCIIgCBaJCFlBEARBEATBIhEhKwiCIAiCIMAS+T9fGfXl2jQ5fwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(7, 4))\n", "for c, label, color in [(c_wrong, 'WRONG: random slope on group', 'C3'),\n", " (c_right, 'RIGHT: fixed group + random intercept', 'C0')]:\n", " sns.kdeplot(c, ax=ax, label=label, color=color, fill=True, alpha=0.2)\n", "ax.axvline(0.5, color='k', ls='--', lw=1, label='True contrast (log scale)')\n", "ax.set_xlabel('Group contrast on log $n_1$ noise ($\\mu_{patient}$)')\n", "ax.set_ylabel('Posterior density')\n", "ax.set_title('Identifiability of the between-subjects group contrast')\n", "ax.legend(fontsize=8)\n", "sns.despine()\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "id": "21", "metadata": {}, "source": [ "## 7. When *is* a random slope correct?\n", "\n", "Random slopes are not bad — they are bad *only on between-subjects covariates*. For a\n", "**within-subject** manipulation (one that takes different values across a single\n", "subject's own trials — ISI, session, trial type, stimulus condition) a random slope is\n", "exactly right: it lets each subject have their own sensitivity, partially pooled toward\n", "the group. There the per-subject offset is identified, because the covariate genuinely\n", "varies within the subject.\n", "\n", "So the rule of thumb is:\n", "\n", "| Covariate type | Varies within subject? | Put it in… |\n", "|---|---|---|\n", "| Between-subjects (group, sex, patient/control) | No | `fixed_regressors` only; keep `random_regressors={p: '1'}` |\n", "| Within-subject (ISI, session, condition) | Yes | `fixed_regressors` **and** (if you want per-subject slopes) `random_regressors` |\n", "\n", "You can mix them in one formula. For example, a fixed group contrast plus a random\n", "slope on a within-subject ISI:\n", "\n", "```python\n", "fixed_regressors = {'n1_evidence_sd': 'C(group) + isi'}\n", "random_regressors = {'n1_evidence_sd': 'isi'} # per-subject ISI slope + (implicit) intercept\n", "```\n", "\n", "bauer validates that every `random_regressors` term is a subset of the fixed design,\n", "and warns whenever a requested random term turns out to be constant within subject —\n", "catching the footgun before it costs you a night of debugging non-convergence.\n" ] }, { "cell_type": "markdown", "id": "22", "metadata": {}, "source": [ "## Summary\n", "\n", "* **Fixed effect** = one population-level coefficient (`group_mu`), shared by all\n", " subjects. **Random effect** = population mean plus a partially-pooled per-subject\n", " deviation (`group_mu + group_sd * offset`).\n", "* In bauer 0.3.0, declare them explicitly: `fixed_regressors` for the population-mean\n", " design, `random_regressors` for the per-subject part (default: random intercept only).\n", "* The legacy `regressors=` keyword is **deprecated** — it put a random slope on every\n", " column, which is wrong for between-subjects covariates.\n", "* A random slope on a **between-subjects** covariate (constant within subject) creates\n", " non-identified per-subject offsets, heteroscedastic-by-group variance, and an\n", " inflated/under-identified group contrast. The correct pattern is a **fixed** group\n", " effect with a **random intercept**.\n", "* Random slopes are appropriate for **within-subject** covariates, where the\n", " per-subject offset is actually identified.\n" ] } ], "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.11.15" } }, "nbformat": 4, "nbformat_minor": 5 }