{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3a339c3f",
   "metadata": {},
   "source": [
    "# Power for a four-arm savings-reminder trial\n",
    "\n",
    "**Worked example on synthetic data, by Ahmed Beddi.** No real programme, institution or dataset is\n",
    "involved. Parameters (baseline rate, effect size, attrition) are planning assumptions chosen for\n",
    "illustration and are stated where they are set.\n",
    "\n",
    "The notebook does two things:\n",
    "\n",
    "1. computes the sample size analytically, with the same formulas as the `PowerCalculator` demo on the\n",
    "   website (same defaults, same rounding);\n",
    "2. checks that number by simulating the trial many times and running the **pre-specified analysis\n",
    "   model** (linear probability model with stratum fixed effects, HC2 robust standard errors, Holm\n",
    "   correction across the three treatment–control contrasts) on each simulated dataset.\n",
    "\n",
    "Design: individual randomisation of active mobile-wallet holders to four arms (control, plain\n",
    "reminder, reminder + goal framing, reminder + loss framing), stratified on sex and baseline activity.\n",
    "Primary outcome: *regular saver*, an indicator equal to 1 if the person made a savings deposit in at\n",
    "least 3 of the 6 months after randomisation, measured from administrative transaction records."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "47be36fe",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:45.696774Z",
     "iopub.status.busy": "2026-11-03T14:49:45.696527Z",
     "iopub.status.idle": "2026-11-03T14:49:46.679534Z",
     "shell.execute_reply": "2026-11-03T14:49:46.677952Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from math import ceil, sqrt\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e102a5ee",
   "metadata": {},
   "source": [
    "## 1. Analytic sample size\n",
    "\n",
    "For each treatment arm $k$ compared with control, with control-to-treatment allocation ratio $r$\n",
    "($n_c = r\\,n_t$), significance level $\\alpha^\\ast$ after multiplicity adjustment and power $1-\\beta$:\n",
    "\n",
    "$$n_t^{\\text{raw}} = \\frac{(z_{1-\\alpha^\\ast/2} + z_{1-\\beta})^2\\,V}{\\delta^2},\\qquad\n",
    "V = \\begin{cases} p_1(1-p_1) + p_0(1-p_0)/r & \\text{binary outcome} \\\\ \\sigma^2\\,(1 + 1/r) & \\text{continuous outcome}\\end{cases}$$\n",
    "\n",
    "with $p_1 = p_0 + \\delta$. We then inflate by the design effect $\\text{DE} = 1 + (m-1)\\rho$ (cluster size\n",
    "$m$, intra-cluster correlation $\\rho$) and by expected attrition $a$:\n",
    "\n",
    "$$n_t = \\left\\lceil \\frac{n_t^{\\text{raw}}\\cdot \\text{DE}}{1-a} \\right\\rceil,\\qquad n_c = \\lceil r\\,n_t\\rceil,\\qquad N = (A-1)\\,n_t + n_c .$$\n",
    "\n",
    "With $A$ arms there are $K = A-1$ confirmatory contrasts against control. For planning we use\n",
    "$\\alpha^\\ast = \\alpha/K$ for both Bonferroni and Holm: Holm's first step uses exactly $\\alpha/K$, and every\n",
    "later step is more lenient, so a design powered at $\\alpha/K$ has *at least* the stated power for each\n",
    "contrast under Holm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3831c185",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:46.681658Z",
     "iopub.status.busy": "2026-11-03T14:49:46.681392Z",
     "iopub.status.idle": "2026-11-03T14:49:46.690572Z",
     "shell.execute_reply": "2026-11-03T14:49:46.689271Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " alpha_star: 0.0167\n",
      "    z_alpha: 2.3940\n",
      "     z_beta: 0.8416\n",
      "          V: 0.4375\n",
      "      n_raw: 1832.0950\n",
      "         DE: 1.0000\n",
      "        n_t: 1929\n",
      "        n_c: 1929\n",
      "          N: 7716\n",
      " clusters_t: 1929\n",
      " clusters_c: 1929\n"
     ]
    }
   ],
   "source": [
    "def plan(outcome='binary', p0=0.30, delta=0.05, sd=1.0, alpha=0.05, power=0.80, arms=4,\n",
    "         mc='holm', m=1, icc=0.0, attrition=0.05, r=1.0):\n",
    "    K = max(arms - 1, 1)\n",
    "    a_star = alpha / K if (mc != 'none' and K > 1) else alpha\n",
    "    za, zb = norm.ppf(1 - a_star / 2), norm.ppf(power)\n",
    "    if outcome == 'binary':\n",
    "        p1 = p0 + delta\n",
    "        V = p1 * (1 - p1) + p0 * (1 - p0) / r\n",
    "    else:\n",
    "        V = sd ** 2 * (1 + 1 / r)\n",
    "    n_raw = (za + zb) ** 2 * V / delta ** 2\n",
    "    de = 1 + (m - 1) * icc\n",
    "    n_t = ceil(n_raw * de / (1 - attrition))\n",
    "    n_c = ceil(r * n_t)\n",
    "    return dict(alpha_star=a_star, z_alpha=za, z_beta=zb, V=V, n_raw=n_raw, DE=de,\n",
    "                n_t=n_t, n_c=n_c, N=(arms - 1) * n_t + n_c,\n",
    "                clusters_t=ceil(n_t / m), clusters_c=ceil(n_c / m))\n",
    "\n",
    "def power_at(n_t, outcome='binary', p0=0.30, delta=0.05, sd=1.0, alpha=0.05, arms=4, mc='holm',\n",
    "             m=1, icc=0.0, attrition=0.05, r=1.0):\n",
    "    \"\"\"Analytic power for one contrast, given enrolled treated-arm size n_t (same approximation as the demo).\"\"\"\n",
    "    K = max(arms - 1, 1)\n",
    "    a_star = alpha / K if (mc != 'none' and K > 1) else alpha\n",
    "    za = norm.ppf(1 - a_star / 2)\n",
    "    de = 1 + (m - 1) * icc\n",
    "    nt_eff = n_t * (1 - attrition) / de\n",
    "    nc_eff = r * n_t * (1 - attrition) / de\n",
    "    if outcome == 'binary':\n",
    "        p1 = p0 + delta\n",
    "        se = sqrt(p1 * (1 - p1) / nt_eff + p0 * (1 - p0) / nc_eff)\n",
    "    else:\n",
    "        se = sd * sqrt(1 / nt_eff + 1 / nc_eff)\n",
    "    return norm.cdf(abs(delta) / se - za)\n",
    "\n",
    "P = plan()\n",
    "for k, v in P.items():\n",
    "    print(f'{k:>11}: {v:.4f}' if isinstance(v, float) else f'{k:>11}: {v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c08b957d",
   "metadata": {},
   "source": [
    "**Headline (defaults):** baseline regular-saver rate 30 %, minimum detectable effect 5 percentage points,\n",
    "$\\alpha = 0.05$ two-sided split across three contrasts, power 80 %, 5 % attrition, equal allocation.\n",
    "The printout above gives the treated-arm size $n_t$ and the total $N$ for four arms; these are the\n",
    "numbers the website demo shows with its default settings.\n",
    "\n",
    "### Sensitivity\n",
    "The table shows how the requirement moves with the two assumptions that matter most: the baseline\n",
    "rate and the effect we insist on detecting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f0d1d89d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:46.692548Z",
     "iopub.status.busy": "2026-11-03T14:49:46.692375Z",
     "iopub.status.idle": "2026-11-03T14:49:46.897911Z",
     "shell.execute_reply": "2026-11-03T14:49:46.888493Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total N (4 arms, Holm planning alpha, 5% attrition)\n",
      "baseline    0.2    0.3    0.4\n",
      "MDE_pp                       \n",
      "3         16512  21116  23760\n",
      "4          9436  11968  13404\n",
      "5          6128   7716   8596\n",
      "7          3216   3988   4400\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "rows = []\n",
    "for p0 in [0.20, 0.30, 0.40]:\n",
    "    for d in [0.03, 0.04, 0.05, 0.07]:\n",
    "        q = plan(p0=p0, delta=d)\n",
    "        rows.append(dict(baseline=p0, MDE_pp=round(100 * d), n_per_arm=q['n_t'], total_N=q['N']))\n",
    "sens = pd.DataFrame(rows).pivot(index='MDE_pp', columns='baseline', values='total_N')\n",
    "print('Total N (4 arms, Holm planning alpha, 5% attrition)')\n",
    "print(sens)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242d0205",
   "metadata": {},
   "source": [
    "Two further levers, for reference: without a multiplicity adjustment the requirement would fall\n",
    "(but the trial would no longer control the family-wise error rate), and allocating more people to\n",
    "control ($r = \\sqrt{K} \\approx 1.73$, the classic choice when every arm is compared with one control)\n",
    "lowers total $N$ for the same per-contrast power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "83f0c2f7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:46.900668Z",
     "iopub.status.busy": "2026-11-03T14:49:46.900412Z",
     "iopub.status.idle": "2026-11-03T14:49:46.908766Z",
     "shell.execute_reply": "2026-11-03T14:49:46.908110Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no adjustment        n_t =  1446  n_c =  1446  N =  5784\n",
      "Holm, r = 1          n_t =  1929  n_c =  1929  N =  7716\n",
      "Holm, r = sqrt(3)    n_t =  1538  n_c =  2664  N =  7278\n"
     ]
    }
   ],
   "source": [
    "for label, kw in [('no adjustment', dict(mc='none')), ('Holm, r = 1', {}), ('Holm, r = sqrt(3)', dict(r=sqrt(3)))]:\n",
    "    q = plan(**kw)\n",
    "    print(f'{label:<20} n_t = {q[\"n_t\"]:>5}  n_c = {q[\"n_c\"]:>5}  N = {q[\"N\"]:>5}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd6a0fc5",
   "metadata": {},
   "source": [
    "### Power curve (same curve as the demo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0409ecde",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:46.910481Z",
     "iopub.status.busy": "2026-11-03T14:49:46.910304Z",
     "iopub.status.idle": "2026-11-03T14:49:47.129324Z",
     "shell.execute_reply": "2026-11-03T14:49:47.127331Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAFeCAYAAAD0Vnn1AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAedpJREFUeF7t3QWYVOXbBvBnu+hYujslpVOQDhWRFhQERZFSultKRRAUDCQUCUGQkE7p7mbphl2W7fNd97v/M9/Z2dmc3WV25v5d11w7O3PmnPO+J+aZN500TdOEiIiIiByCs0OkkoiIiIgUBn9EREREDoTBHxEREZEDYfBHRERE5EAY/BERERE5EAZ/RERERA6EwR8RERGRA2HwR0RERORAGPwRERERORAGf0REREQOhMEfERERkQNh8EdERETkQBj8EVGyOnv2rBw5ciRZt2FO0zRp3769zJkzx6rtPnv2TF68eGHVOhLryZMnEhQUlOCPBwYGyu3bt+P1ubt370pERES8liVKSu3atZPZs2cn5SopARj8kUXFixcXNzc30yNDhgxStWpV+frrryU8PJy5lsplyZJFhg0bliKpaNmypdSrVy9Fg4yVK1fK77//LqVLl05UGo8dOyYVK1aUTJkyqXO/Zs2acuXKlUStK6F27dql9jtz5sySPn16adiwody5cyfO1ezfv19q1KghadKkkVy5cqnP43q15JtvvpHcuXNLjhw51HK9e/eWsLCwOLdhrmfPnuoeUadOnRg/i0C8SJEiarlRo0aZlhs0aFCUe0zatGmlSpUqap3Xr1+Ptj581nhPMj4vXLhwgvc9IWbNmhXjtm3xfrhw4UK1v6dPn7aYzK+++kq9/+jRo4Rkg1ofPvfPP/8k6HOWZM2aVYYPH65+YFHKc035TVJqgBta2bJlpU+fPmp3/f39Ze/evdK/f385evSo/Prrr6khGRQDHN+UKvHBFyfOH2fnlPutOXfuXMmfP7/UqlUrwefAmTNnVDBTqlQpWbFihQQHB8u4ceNUYIUSTARMyeW///6TN998UwWbkydPVqV/I0eOVNtGQJouXTqLmz5//rwKEosVK6ZKU7Dvf/75p/Tr108Fg927dzd9btq0afLFF1/IZ599Jq1atVJpwpfww4cP5Y8//khQ0vRzaOfOnaqEt0SJEtE+v3XrVrl8+bI4OTlFOefwWQScv/zyi/pMQECAnDp1StasWSNLly6VBQsWyLvvvmtan778jz/+qAIQI6QxOV27dk1tw1IwnZLndXwh4I4tmI/r/dhgvfi8tbp06SLfffed+pGGgJ9SmEZkQeHChbXWrVtHe2fs2LG46rXTp09bzLfg4OBky8+goKA41x0REaE9fvw4zuXiC9sMCwuL7+ImL168iPUz4eHhmr+/f6zLQEzpiU8asd8hISEWt5EhQwZt8ODBFt+LTz5bEhgYqPY3qYWGhiboGDx8+FBzdXXVRo4cGeOuIF9iOq7dunXTcuXKpT19+tT0+Vu3bmk+Pj7a8OHDkzp5UbRo0UIrUaKE9vLlS9Pr58+f11xcXLSvv/46xm0PGTJES5s2bbRzqnr16lrVqlVNn8MxypQpk9a9e/co65o7d666rk+ePJmg9GE9pUuX1rJmzaoNGDDA4mfbt2+v1ahRQ/Pw8IiSfwMHDlTbNIdzu0mTJlq6dOm0u3fvmt5G3mP5xJ6f1njnnXe0ihUrWrMKxXhcYxPb+Rkfv/zyi8qrU6dOWVx88uTJ6n1cK5bElMdYHz63du1ai58LCAiI9ll8J8R0r8O5/sYbb8QnSZTEUu6nONmFt99+W6UDpRC6e/fuSdu2bVUVmY+Pj5QrV079etctWbJEVekYqyAuXbqkXuvUqVOUfMF68HkdqiU6dOigqqa8vb1Vicbq1aujfKZFixaqfdf3338vOXPmVCWWlqA0CNtECYURfsWixKJNmzaml7ENVL15eXmpR5MmTeTGjRuxHsPHjx9L165d1TaQD6hWwzaNsO3GjRuLh4eHWi5Pnjwyb968KMs0bdpU5ctvv/2m1oF8bd68uSr5mD9/vhQqVEi9li9fPtm2bZvps6iWxDo3bNggH3zwgdoHlBTVr19fLl68GOf5h2P22muvmT6HfUCbsLj89NNPap9wfJBXOIbICx1Kb1CapcNz7Kelh7EKCCXMKLnT8xMlWzFVYxlt375dlU4Yj+e+ffvUelCy9s4776g0okoV546xfRza2KHUCWnH+zqcV6i+Rlpjg+pZlIrH9IitDd/9+/dl3bp16th5enqaNlO0aFF1DGPbNkrubt26Fa0ELFu2bKr0UIfzGscG1bxGnTt3Vp/9+eefY02fJShFRp4izSEhIVEWwbZXrVqljkV8q5UzZsyoznMsP3HixATvjyW47mI6Jih5iguqofPmzRvXYtF8++23qokFSj5RpY1rBNWdn3zyibx8+dK0/O7du9X5efDgQWndurVaDucfSmbjU+WfFHB/mTBhghQsWFCdf7i/DBkyREJDQ2NcvX7PQTUw2vDhPo1mEosXL1bp01/DfQjXz/Pnz6OsC+cNSo1Ruk4pLImDSbITMZX8HTx4UP3yW7VqlUopftHlyZNHy5s3rzZv3jz1+nvvvaeWWbhwoVrm/v37mrOzszZt2jRT7syePVtzc3NTpRAoBQP82kXpRf/+/U2/kgsVKqQeCxYs0NasWaNKEbDuDRs2mNbVoEEDVVJTpkwZ7bvvvtNWrFhh8Sg8ePBAc3d314YNGxbl/Z07d6p1/vnnn+r1Q4cOqZKjDh06qO38+uuvKj9ee+21WEu2mjdvruXIkUP7/vvvta1bt2q9evVS68V+A0oxkF6kZ/r06dpff/2ltW3bVi2zadMm0z7VrVtXK168uMrH9evXq9JWlDp16dJFq1y5svrc/PnztaJFi2rZsmUzlSZcunRJrQv5gJIT5MNPP/2k9h0lfTgOOvOSvyVLlmhOTk7ahx9+qP3999/aDz/8oPK0fPnyquQtJigBwDZ79+6tbdmyRZUgYZ+aNm1q+gjypWbNmqb/sf9YTn/MnDlT8/T01CpUqGAqjTt69KgqKWrUqJG2fPlybenSpWpfMmfOrN27dy/Wq2zQoEGat7d3lJK93bt3m/IG+Yp9wHHC8cL5ixIxuHjxoloO75v79ttv1XuxleoiD3DuxPQwlmSZ27dvn1o/9tUcSjG9vLwSdHfZu3evygdjCe+4ceNUXluCUkKUPCYESv5wHur5+8cff0T5+KxZs1Sp5e3bt9X5ZSw5jankT4eaB2OppTUlf998802MxwTna1x8fX21du3aqX2uV6+e1qlTJ2316tVxfUzd83AuoHQU5+W6devUMcC53axZM9Pnt2/fbjo/cd/BunF+Yt/y588f7xJDXWJK/nDtY19xf8R1rZ8rrVq1Mq3XvORPv+fgOE2dOlXtN5bH/eWtt97S+vbtq+6hOAfxHaDf23VYHuf8gQMH4sxLSloM/sgiS8EfggB8ceLLA9VggC8l3DAuX74cZT1vv/22+mLVqw/wxdKwYUPTMi1bttQ++ugj9cWAYAv0G+C2bdvU/7iZIBjUt6VDsIcAx/g/bjaPHj2K82i2adNG3UyN1ZPYDwQVepX1lClT1A3JWO135swZFaDFVFWI9eGG/sUXX0TZBwSUfn5+6rXnz59rP//8c5S8Qv5g341VZgj+kHfGKnTcRLFP586dM60fQRpeO3z4sHpNvxEjSNIDakD+me+bMfhD0I0AqGfPnlH2HQEYvrAR/Mbks88+09KnTx9lezieO3bsMH3EPPgz98EHH6j9MZ5DCPrwhWI8TqgORPDcr1+/WI8zvnwQSBrpwQnWa4RqTqQRwTjg3EOeIhAzt2zZMvUeqmGTA84VrB8BqLk5c+ao9+JzjqMaDcEGflzNmDEjSh7iGOfLl8/i7uOaxbmTEAj+cK+AIkWKRLnGAetr3Lixeo79T0jwh0ArY8aMpt3Rgz8EFbiOjQ/9x2hywA8D7DvuVW+++aYKkkqVKqVeM/8haQ7BH5YbNWpUlLfw4w2v79mzR72u3/uMASEcO3ZMLYfgNSH04A8/FnE/MX8ULFgwSrUvtoPrAOeL0eLFi9Vy+FEDMQV/X375peljuEfiejb+4ANcl+bnF655Y0EBpRxW+1KMUO3Wo0cP9UC1KnoAoxH5lClTVDUYoDF3gwYNVFWBERrworoCjecBRf4o3scwFKhGQHUlqpDLly8vmzdvVsuguhJVPmjsDlu2bFE9LtFZAA3a9Qcav2O9xiomVFeiKjQuqFJD4+09e/aoRbEOpAnVfO7u7uq1MmXKqL9YFtWFqNZCtTCqHV1cXCxuAo3ZUSWNxuuoxn3w4IFaDtVdqLoFVOugWhh5hepN5C86FOCzqBYywrr0/YHs2bOrv2jQr9PXq29Lh7QYG6HjWKGqVU+zuQsXLoifn5/KV2M+owoX2zBW8ZtDXiEtH330kVoO3/E4ZrVr147rUCio3kN1JvJNP4dQ/YSqW1R1Yt/0fUK1aKVKlWLdH0BTAVStWYJG5kao2sc69bzRqyYtHWf9tdiqwayRVNvGOVKyZEmVB+j8getIh23EdA7jdWvShvsEruWrV6+qzaEjiX4PSQzsj6WqYlSdmz8sLZdUcB3gXoaq340bN6pz9vjx4+panjRpkpw8eTLOTZmfd6hmd3V1jXZNmi+H+xrukTFdu3HBtY+OT+YP3GeNcL2BsWMQvPfee6oa13gOWVKhQgXTyzhuqOo23qsA9xLze5V+naKzEaUs9valGGF8M7TNA7RBQfsTtOdC2xUd2nyg3ZQ5/YscQQ1uXgj+Bg8erAJAtC3CzRtBHoYAwRcGhn34999/VXs33BQB68aXP4JOS/Alo99gzHv/xQRBEIbBWLRokdo+2lihXVK3bt1MH0H7vpkzZ6relhgyBEFbs2bNZMyYMar9VUyWLVsmvXr1Uu158BdtDxEU4X8EeIA2QOjhhjZ4CNAwBIal3oJ6HsRGX6c5pM8c2isdOnTI4ur0IUzMv3h0sbUXRICMtpAIMtA7EwE42iyNHTvW4n4YHT58WD799FPV8xTnlg5t19AGCO29LLX50oPemOCz5l9uupjyRs8D/X0EmubQthVi2z5+LCBPYoJhXPDFaIlx2wUKFIiyCLaNc9zX1zfGdev0sdOQD7he0a4K1wqODbZhKW2AbcSVt7F5//331fBBCObROxp/sb+W7g/xgXuHeQABaEOINrMJgaFPEKhZguAGxy02+o8vHQIcrA8/WvADWP/BGBPz8w4/7NAeUw+U4zo/zZeLL7Thww9Jc/gBj+tPh/Mf5wfudUZIJ9r+mf84NWfpfmV+f7J0v9KPY2LGsyTrMPijGCE4wo02NriBGRv36/RfcvpNEyVnGIsLN0o0tq9evbr6i2AMwRA6FqCU4MsvvzStA18c+ExMDd3RuDihcDPDlxQ6h2C76IyCkir8wjbq27evGuYG+7Rp0yYVtGGfEQzjy8ISdHpAAItftwhy0cAcwQ1KLhH4ogPH559/rtaLBtHYJr4EMH5iUrL0KxqvmX+B6fSA4u+//5bKlStH2xVjCaQ55Ce+6DEcCRqro+E3go8dO3bIuXPnYgxicc6gVBTbMw/wEBzhiwKvG4NyXVxDa+BLLKbxy+LKG3TAwfpRiomA3wilPWiEb+wIYg55aexoYs7YkcMcvmQB2zb+wAJsGwF8TAG/JfhixdBM6Mhz4MAB1dEI28CQKvgyx/mqQ8kZOtPEtu9xQdoR+KPDDMbkQ4k6rrX4/jAzwg8ylDbpHcyshR9ZMaUNP2wTA6VWSFtMwbQRzjtjAIgScpyj5tdkXOdncsF9HCX4qOUwLxmO7d5hLf06jekHESUfVvuSVVDcj5I781kQ/vrrL/XlY/zVidI/lHygZ5v+xYrqQXzZYiw43HTwBaVDiSF+naKkEDcn/YGbEapr41M6ZgkCCv3LRe8Za4RqYVSDYL8QGKLHG4I/3KhOnDhhcZ0oZUFgi8APXwoobcGXLkou8Dog3fjyxQC7CHoQVKGaLa5f1QllPgAret3hGMXUCxqBOY4VSgaN+Yx0oMo1ti9HBHjoQYsvQQTH48ePV4Eg0hRT72h88aFqGk0AMK6c+XFENRtKe7Fe4/7gOdYbV6kPvujRu9PSWGTmeYPjhcGR9bzBuYZzEyXDxipQBEzLly9XvRdjg4AKwXBMj9jGo0NwV61aNVWaZNx37CNKqGPbdt26dVWJtTmUooIeZKCEFQGoea9eBP44v+NKX1xQ0o1jhJlVEBTFVgoaE1zb+JGE8wPXXlJAvsZ0TOLaBmofcC2gZN8I1zPOkZiuKyPz8w7XI0q7zD9rvhx+FONHVXy2YQ3cx/EDYP369VFWg/svziHcB5OD3nsf1yylsJRrXkipSUy9fc2hwS46E6BR/5UrV1TvX/TMxWvm48ihUTMaoadJkybKeGLo/IGG3cYOIYCOEuitiEbQaASPMbI2btyoGhO//vrrpl1BBxA8EqJOnTpa2bJlVeeVJ0+eRPkoOjGgcffKlSvVNjGmIcb5wrLoMWwJOlWg4wt6S167dk179uyZ9vvvv6semujlB2j8jfSjcTo6eqDXKnoNYlvGMRXROFtvJK/Tx+Uy0huD6z2f9Q4f6IAxceJE1WHlxo0bqiE/etoZx3Az7+2LceKQPhw7HEPk/SeffKIaY6OjQ0zQ+B7HE72b0UHlyJEjav/RgUbvJWze4WP06NFqvePHj9f++++/KA90igH0GsUyI0aMUGORoaMDxrnDa0OHDo1xf2DRokVqOeMYZ3qHD3QYQS9KpBHnFPYVx8jYqQjnKfLr3XffVetAh5r69eur/IlpfMukgnMO+45OBejcg/xBo310fLp586ZpMz169FCdqPROQegchc+hAwY6reD4obE+xt/DtWLsqISe2Tjn0Hv5+vXr6nzE9We8puCrr76KdWxBMHb4AHTQQWcCjNFn3uA/pg4fODZ4oIcptonOXNg/804OeocPdD7QP6P/Rc/m5IQORNmzZ1e9wHGuoLMVRjjAOHWxjTuod/hAnuA6xbLYV3Q6w2f146J3+MB1g1ETMF7ehQsXtFq1aql74J07d0zJQ8cgdIaJaQxPSExv32rVqqke/tgX7Of+/ftVT26MPKBfyzF1+ECPfCN8f5iPJfnpp59quXPnjrIc7gW4byK9lLIY/JFF8Q3+AL1gMSQBbu7oMYa/CKDMb4q4geDLFz1LjfQBZjEshDn0GsVN1rhufBniS0uXmOAPPVixTgyrYA7BIIZzwPsIAvA3S5YsMQ4ho8MNF4GQ8XPo4aYPI4JAUl8vbnjIKwRYCESTMvhDfuKLRT8W+PLAUC5G5sEfgggMw4DP6J9DsINhZmKDgKRcuXJR0ozjiy8QnXnwh7xEHll6GIc5wbA9CMz0Y49j8f7778c5yC96BWNInzFjxpj2QQ/+EBhi//Q04nzEcDrmEIQh4NL3EUGUsQdzcsLwPHq6sX3sL4JqIxxf/MAyDtKLL1L9/NP3G+efeW95BBwIHvVzG8siCDYOHI4gDsML4fyMjXnwB/jhgXXiejCKKfgzngc4HuitjCGlzOnBn6XzBulOTrjfII+M20Zghh9csdGDv3/++Uflp37e4fidPXvW9FE9+EPAjgBMXw7BIH7wGqHHMX64xiYxwR/OEz2N+vZxrzX2wk/q4A+90jHcFaU8J2wyhQsbKRVAlR2qJWNqJ2YO1Z7o9YaqI1TXxtQwHVUIaLek9xYGVH/gdbyGKj9zqI5ANR4aPaPXmnkViD5Ir3GdcUF1DdKItiYxteHCVFOowsIySFN82gah3QyqT5EmVF0a21XpUI2D6hy09UN1EvYfVcx6XuN/5JFxGrGnT5+qvDWuD9Vj6KWL5bBv2Fe0kUT1JBrZ6z2V0ZvVfFow5CVeM283iW2gzRl63KIqKL7tKtGzE/uCNCCvjO0E0VMS69PbPGHbMU0th2WM7eJQ3YqqdqQfjerRJi8+MLA02smh2hp5id6SaMOKqjrsH/IG1cdIY0zHFe2fsG1UaaP5QkLa21kL56c+nZulTg9oR4q8MZ9qDs0vkGa8j8/Fdv1iwF3kEarczJfDsUQ1NJpvGDvjmEOVNJoVGAdA1q9nvGZs74fzEx1x9F75ONdwXHV6846YoJ2occBqIxwb8xEHkgPOXTQLwXUYn0Gfp0+fLgMHDlT3MOw/rn0cM5zLxuYOaCOLqns0dUAPdJyfuA4snZ9oOoKOPFOnTo0xiTj+qHbH9WKpza5+P0HHImMbWoQD6PyB5hy4l6CDm/G8N7/n6PdRnD9ow63DvRX3cmOvezTXwTmLezigPTXSh7QmdbtnihuDPyI7YQz+0ObQkeGLBe2U0JYNbfiMwR+GtKHY6T3vEaDFNJ8wxc0Y/MU0xA4Ygz/zzj5GCO7RIxed0DDEVmqGzkAIMtHmllIeO3wQkd1B6R46k8TWwYJihi9ldJJg4Gdb0PkEJYL6WKipGUo6zae/pJTDoV6I7ASGoEEVU2xjEToSVP3qUJ2MvLE0zh9FhwGbzQf8pYRDkwucd3E1GUD1LJYzr8Y3h+pijMcY25BBqYWlYZwo5dhFtS/aKcR1cRERERGR2Ga1LxohoyEofgWZTwdjhEFz0TgWDU1RhIxOAUREREQUM5tr84eZD9AgG0EdekTG1CsQA+hiQFHMCIGehuiNh3YQlmabICIiIqJINlfti5Hd33rrLdVFXA8AMbK/OQR6derUUXMXAqp+MQQIJsw2ThGmw1AkeOgQVCJQxFAWrDImIiIiW4K4BsP2YBizuKa1TLCUH1owduHh4WoBfTDJu3fvRvsABg/GIJSHDh2K8t6wYcO0Ro0aWdzAqFGjYhxYNqYBZ/m65YF4mS/MF54DPAd4DvAc4DmQMucAZuxJajbX2zc+0S0GwUVEbD7gK/7HoKSWYP5GTHKuw2C8GKQTg5lyOAMi0gUGB8vrnw9Q/x74Zrp4e3gwc4goxWEgdsQ1GNsxqdlc8BcfGFEezEcux4j9+nvm8J6lCeER+DH4IyKda3CwOLu7q39xb2DwR0SvUnI0TUuVwZ8+Vhemr8HYZjqMf5Q7d+5XeYyIKJXzcneXM/Nmq1TgORHZD9QahgcHSejLAAl7+ULCggIlFH//99z0f9BLCQuKfD/ygf/154FqHfi/Uq/RkrNS3VSXQaky+MNcrAjytm3bFmVAW/yPTiJERImFX9ks7SOyLRHhYRIa6C8hAc8lNDBAQl9E/g154a9eV6+Z/iKww///e6jgDs8RzL0QLSIiyRIX/PxJkq0rJaWa4A+TaWN4l08//VS1C+zVq5dMmTJFTTqOSaXXr18vmzdvlq+++iol84+IiIjigBI3Faz5P5WQgGemR7D/Mwl9of//PPJ1BHYvIgO9kBd431+VsiU1V09vcfXyETcvn8jnnj7i6oW/kQ83L29x8cBfvO9leh3vuXjgfy9Jn6dIUu9Wikg1wd+iRYtkxIgRapL2QoUKqalwDh8+LPnz51cNIq9evSrffvutGu6FiCixQsLC5Ju/1qiPf966pbi7pprbJFGKCA8NkeDnjyX42ZPIv8+fSLD/Ewl+9liCEdz5P1V/8br+HEGdFh5m9f4h6HL3SStuPunEzTuNuHmnVX/Va95pI4M59XoacfOK/Ouq/iLIw/P/D/acnG1uqOMUY3Pj/OmCgoLk9OnT8tprr4mrq6v6H5Namwd3eO3KlStSsWJFNTZgfKEXDaqP0euXHT6ISIfeviV79lb/ou0fq4DJ3kWEhUrQ04cS9PTR//99Fvl/8LNHEoSgTv3F/09U9WpiObt5iEfa9OKeJoO4p8HfdP/7G/kcQZ16zcfwv09atQwCOWdXt8RuOtVJzjjFZn/SYuJqBHQ6/G+pVA+T2HMieyJKKi7OztKtYQO1OjwnSo3Qrg0lby8f35eXj+/JyycPJOjJA3n5BP9HPkfAh9dROpdQKDVzT5tRPNNnEve0Gf73N6N4pIt8ILjzSIcAL/KvR9oMajlXD6+EboqSgc0Gf0REr4KHm5uM6tjuVWyaKE6Rbef8JfDBHRXUBT68K4GP7srLR/ci/8dfPJ48SFA1q5Ozi3ikzyyeGfDIEvk3fWbxSJ/J7C/ez6RK5hy52jS1Y/BHRERkI9Ap4sX9W/LiwW0JfHA78u/DO5FBnvp7Rw1LEl8I6Lwy+YpXxqyRfzP5RgZ3GbOIV0b9eVZVMsdgznEw+CMiIkoBKLVDp4iAe36RAd79m/Li3k0JuHczMtC7f0t1jIgPtIHzzpJdvDJnF+/M+JtNvDNnU3+9Mv3vkTGLQ7WRo/hj8EdEZMAOH2QNDAYccPeG6eF/J/L5i3t+EnDXL15DlqjALmtO8cmaU7yz5oj8myXyuXeWHCrYwzAjRInF4I+IiCgBMLCw/+2r4n/rqvjfvib+d66rvwF3rqsOFnFB6ZyPb25Jky23+uuTLVfkXwR8vrlUr1ai5MTgj4jIAFO6Hf52pnqF07s5rvCQYBXUPb95WQV5z29eUQ8EfeglGxsMV5I2Rz5Jkz2vpMmeR9LkyCtpsuVR//v45hQXd0/HzViyCQz+iIgMML1b5nRpmScOAm3snt24JM/8Lqq/CPae+11SVbWxTQOGThNpcxaQtDnzSdqc+SVNjsi/+N8jbUYHyT1KrRj8ERGR3UOQ9/TaeXl6/YI8w+PGBXl6/aIEPYm5mhYzRqTNVUDS5S4k6XIXlHS5Cqj/EfRh4GGi1IrBHxGRAaZ3++GfDeqVj5o25vRuqQyqa1F69+TKWXl67dz/HufVECkxQa/Z9HmLSLo8RSR9nkKSLk9hFfBhWBSUBBPZGwZ/REQGYeHhMm3lX+qVDxo1ZPBnwzD92JMrp+Xx5TPq79OrZ+WZ3+UYBzdGT9kM+YupQC99vqKSIV8xSZenkLj7pLPhVBIlPQZ/REQGmNKtXe1a6hVO72YbMD4exsF7fPmUPL4U+Xhy+ZQa+NgSDJWSoUBxyZBffxRTQR+DPKJIDP6IiAwwvdvkD95nnrwiKtB7eEceXTgujy+elEcXT8jjSyfV4MiWoA1exoKlJGPBEqa/KOFjdS1RzBj8ERHRKxPs/1QenT+mgr2H6u8Ji50wnFxcJUO+opKpcGnJWKi0ZCqEgK8kx8QjSgQGf0RElCIiwsNUj9uHZw/Lw3NH1QOdM8w5ObuoNnmZi74mmYuUkUxFykrGAsU5Ph5REmHwR0RkgOndKn7WT71yeNZM8fbwYP4kUsiL5yrQu3/6UGTAd/6YhL18EW1tGCMvS/HyKtjLUqycqr7l9GVEyYfBHxGRmZchIcyTREBbvXsn98uD0wfl/umDapgV0bQoa8LYeVmKl5MsxStI5mLlVdDnmT5TIrZGRInF4I+IyMDTzU12TZ2sXsFzsgwdMzALxr0T++T+yf1y//QBCbhzI9rCKNXzLVVZspasJFlKVJD0eYuKs4sLs5XoFWLwR0Rk4OzsLHmyZmGeWBBw10/uHt8j9078p4I+DL9i5OTsLBkLlRLf0q+Lb6nXVcCHgZKJyLYw+CMiIotePnkgd4/tUY97x/eqkr4ogbKrm2QuWk58y7wu2cpUkSwlKnHaM6JUgMEfEZFBaFiYLNyyTb3S5Y164ubqOLfJsKBA1Wbv7tFdcufobjVjhhGGW0GnjOyvVZdsZatJ1pIVxdXTm+cPUSrjOHc1IqJ4CA0Pl3FL/1BLtq9b266DP7TbQ4B3+/AOuXN4h+qkEREatbMLxtTLXq6GeqDtnpt3mnjkIhHZMgZ/REQGzk7O0qpqFfUKntubkIBncufILrl9aJvcPrRdXj6+HyWJPr65JHv5WpKjQi3J/loN8cyQ2d6ygMjhMfgjIjLwdHeTb3r1sJs8QenesxsX5daBLXJr/2Z5cOawaBHhpvS5eHipatwcFWpLjop1JF3ugpwajcjOMfgjIrIzEWGhqkfuzf2bVcBn3lEjXZ7CkqtyPclZqZ7qmevizoGsiRwJgz8iIjuA2TRuH9wmfvs2qerc0BfPTalydvOQ7OWqS67K9SXX629Imux57CDFRJRYDP6IiAwwvVvNgYPUK7unTbHp6d3QXg/Bnt/eDWooFpT46TwzZlWBXu6qDSVH+ZrslUtEJgz+iIjMPPYPsNk8efHgttzY/Y/c2PWPPDh7KMr0aajOzVOtkQr4MG0aBl0mIjLH4I+IyABTum2aMEa9YivTuwXcuyk3dq1VQd/Dc0ejHK/MxcpJnuqN1SN9nkI8lkQUJwZ/REQGmN6taK5crzxPVAnfrnVyfeffUQM+JyfVSSNvzWYq4PPJmuOV7ysRpS4M/oiIbETQs8dyY/c6ubbtL7l/6sD/75WTk2QrU1Xy1momeWs04Xy5RGQVBn9ERAaY3m357r3qlTY1qyf7DB+YUs1v30a5uvUvuXNkp2jhYaa9yVqqsuSr3ULy1mwq3pmz8TgRUZJg8EdEZAz+wsNlyC8L1SutqlVJluAvIjxc7h3fI1e2rBS/PetVAKjLVLiM5K/bUvLVaSk+WXPy2BBRkmPwR0RkgCndGpYvp15J6undnl6/IFc2L5erW1fKy0f3TFtNkyOfFKj/luSv25qdNogo2TH4IyIywPRuP37+aZLlSbD/U7m2fbVc+fdPeXThuGm97mkzSP46LaRA/XckS4kKnFKNiFIMgz8ioiSmRUTI3WN75PKmP+TGng0SERqstuDk4qqmVSvY4F3J9Xp9TqtGRK8Egz8ioiQS+PCOXNr4h1zetExe3PMzrTVDgRJS6M22UqBea/HMkCWJtkZElDgM/oiIDF4GB0uDoSPVK5snjhUvD49Y8wedN24f2iaX1i+RWwe2qFI/cPNJJwXqtZJCjdqpThxOTk7MZyKyCQz+iIgMNBG59eiRegXPYxL46J5c3vi7XFy/RAIf3DYthgGYCzfuIHlrNRVXDy/mLRHZHAZ/REQGHm5usnrkMPUKnkcJDDVN7h3fKxfW/qbG5tPH5EPnjUIN35XCjdtL+rxFmJ9EZNNsKvi7e/euLFiwQM6ePSsVKlSQnj17io+PT4wZeOXKFZk7d65cv35dihcvLp9++qlkzZrVpjOciGybi7OzvFawQJSdDA0MkCtbVsiFv3+VZzcumt7LWrKyFGnWSfLVaiou7p62nTAiov9x0vBT1gbcunVLatWqJYUKFZJ69erJihUrxNXVVbZt2ybe3t7R9vD06dNSpUoVqV69utSsWVPWrVun1nH06NF4BYDPnz+X9OnTy7NnzyRdunQ2kANEZGue37oq51f/LFc2/6kCQHD19JaCDd6RIk07S8aCJWxtl4nITiRnnGIzwd+AAQPk4MGDKthzcXERf39/KVu2rAwePFiVAJrr06ePnD9/XjZu3KjeCg8Pl2LFiknfvn1VCWBckjNTiSj1wvRui5ctktuHd0iuM1vFRYvswJE2V0Ep1rKrCvzcffiDkYiSV3LGKTZT7fv777/LlClTVOCnbrRp00qHDh1k6dKlFoM/lNphGR0+hypivG5JcHCweuiQqWB83dnZWdzc3CQ0NFQi/tdjD7BulEKGhISIsaAUr+E989exDqzLuD3A6+jxh+WN3N3d1XqxXSMPDw+1H8bX8Xksj2A3LOz/5wDVX8dreE/HNAmPE889Ja7rKTgwQK7t+FtOrvlNRofkwZUvI8RZcr5eX4o07STZX6suLq6uvEfwvsd7eTyuJ34/WR9HJCebCP4QDN25c0dKlSoVJa0lS5aUJUuWWEz/+++/L40bN1YlhnXq1JGVK1eqQKxFixYWl580aZKMGTMm2nszZswQT8/Itjrly5eXli1byvr161X1sQ7rr1u3rixbtkwuX75seh3bQtvE+fPny4MHD0yvd+zYUQoXLqzWbQz0Pv74YxWcTp48Ocp+oHQT1c/ff/+96XUEeEOGDFHtGhcvXmx6HW0aP/nkEzl+/Lj8/fffptdRXd6pUyfZvXu37Nixw/Q60yQ8TiI892K5nj7o0FY2LVkgRx8EqesmLEdVyX37tmT0chPfVlNl74nTsnfjbpGNu3k98XpSeC/n91NKfOf6+vomW/xnE9W+CPxy5sypAp0CBf6/ofU///wjbdu2lYCAyLY2RoGBgdK/f38VeCEgun//vowaNUqGDx+uSrvMWSr5y5Mnj/qc3uaPpWQsJQOW0MZdSmYPv+qf374mF9b8Itc3/ynhYaGiObuId5acqgMHeu6myZiZJems8WAtDmumXlltG2Kf5OqbYBPBH6o1USVz6NAhVZKmW7RokYwYMUKuXr0aLZhDKRdeX7VqlYqOL126JI0aNVKvWyrhM5ecdelEZLsenjsqp//8Xvz2bhD5X3+3zEVfkxLvfCR5azYVZxebqBAhIgeXnHGKTdzlUNKSI0cOOXnyZJTg79SpU5IvX75ohx/t49AbGNWwerEoqll79+4tP/74Y7yCPyJyHLhn3Dm8Q079MVvun/zPlPBcVRpIyTY9xbd0Fc7AQUQOI3r96CuCzh0Ys0+vTnr69Klq64b2czo/v8i5MtG5IUuWLKq3rxH+T846ciJKXTD12vVd62T9Z81k6/DOKvBzcnGVgg3aSPO5/0q9MT9LtjJVowR+mN6t4dCR6oHnRET2xiZK/gDt91CFi44V6GDx119/qTZ5qMYF9PpFgIihYLDMhAkTpEePHmqol0qVKsnOnTvl8OHDsmHDBns7RkSUQBHhYXJt219y6vfv5PnNyE5aLh5eUqRJB1W965M1Z4xrxMCnF29HTtdmE4OgEhElMZto86dD54tffvlFleCVK1dOPvzwQ9MAzxjAec6cOfLFF19IhgwZ1EfQIxdDxKDDSP78+aVr165SsGDBeGVRctalE9GrER4aIlc2r5DTy76TgDs31E64p0mvxucr1uoD8UyfKc4dC4+IkAPnL6jlXi9WVM34QUSU0pIzTrGp4C8lJWemElHKQtB3edMyOf3Hd/Li/i21cY/0maTE2x9J0eZdxN3n/8cEJSJKDZIzTrGZal8iooSyFPR5ZvRVnTiKNuukpmIjIqKoGPwRUaoTERaqqndPLvlGXty/qfbfK5OvlHrvUyncuJ24englOk1h4eGy5dgJ9fk3ypUV1//NOkREZC8Y/BFRqoHeu9e2r5YTi2ZIwJ3r/x/0te0thZu0tyro04WEhUnPWbPVv2fmzWbwR0R2x6rgD3PI3r59W/XKNXf37l01tVnevHntLtOIKGVhnD6/Pevl+MLp8uxGZGcMj/SZpfR7vdWMHNaU9JlzdnKSioULqZfxnIjI3ljV4eP69etqmBXjvLY6zH2IMffiM9vGq5CcDSmJKOncObpbjv08WR5dOK5Wit67Jd/9WPXgdfPySboNERHZkOSMUxJd8jd06FC5d++emnuuX79+UbIL8+4uX75cDc1CRJQYCPaO/jxF7h7dpT6Ozhsl3u6hHggAiYgocRId/CHwQ4lfRESEquI1cnFxkbFjx0rbtm0Tt1dE5LD8b1+TY79+Jdd3/K3ywNnNXYo07SRl2n8mnhmyOGy+EBElFauqfTG4cufOnWXz5s1JtT8pJjmLU4ko4YKePlK9dy/+s0j15hUnJyn4xjtStlN/SZM9ervi5BIUEiJtJ32lVr9syJfi6e6eXJsiIopRcsYpiS75gxw5cqjAD1W/adKkUa+dPXtWLly4IPXq1ZN06dLxsBJRrMJDguTsqgVqgObQwAC1bM5KdaX8B0MlY8ESKZ57EZomJ65eU9vFcyIie2NV8AfffvutrFy5UrZv367m161fv76q9s2XL5+aazdtWo6sT0TRaRERcm3Hajn28xTTrByZCpeR8h8OlRzla76yLHN3dZWf+vZR28dzIiJ7Y1W17+PHjyV79uyyd+9e1ev3nXfeUfPuzpw5U+rWravm5u3du7dN5llyFqcSUewenDkkh+aNkUfnj6kFvbPkkHLdBkuBeq3FiXPpEhEla5xi1c9aPz8/Ffwh8EPHD5T+rV69WlX3Nm7cWM6fP8/DR0QmmHf3yIKJcn3HGvWaq5ePGquv+Fvdk3SsPiIiiplVwR8GcEZP30OHDql2fmFhYVK5cmW1tWPHjkn16tWZ90QkYUGBcvrP7+XM8rkSHhykOnMUbtROXnv/C/HKmNWmcig8IkL2njmr9ql6yRLiwpJIIrIzVgV/GTNmlIEDB8rrr7+uRuAfN26ceHh4yMiRI2XPnj0ye3bkFElE5JhwX7ixa60c/nG8BD64rTLBt0wVqdRztGQqXNomMyU4NFQ6T5up9g3Tu3l7eNjkfhIRJZZVwR9MnDhR2rVrp6Z6K1++vHqtRYsW0qlTJylQoEBi94uIUrknV8/Koe9Hyb0T+1RKfHxzS4UewyVvzabiZMPTpmFKtxL/m7KS07sRkT2yqsNHasYOH0TJIyTgmZqD98LaX1WPXhd3Dyn13qdSsk1PtusjIoqn5IxTrC75e/bsmcydO1eOHz8uoaGhUZLUpUsXVQpIRPYPVbxXt6yUIwsmSNCTByrBKOWr0GOEpMmW2/4zgIgolbA6+GvatKnq9JE7d265dOmSGu5l9+7d6osAw74Qkf1DFe/B2cPl/qkDKrHpcheSyp+MkxwVatl/4omIUhmrgj8EewcOHFDz/Do7O0v+/PnVoM8Y9qVWrVrizmmRiOxa6MsXcmLRTDm3ar5oEeHi4uElZTv2VUO3uLilzmnRML3b+9O/Vsft1wF9Ob0bEdkdq4I/wLRu6dOnV7N6uLm5qfl+Me1bgwYNZM2aNVKlShW7yzQiEvHbt1EOzhlp6sWbp0YTqdRzlPj45krV2YMp3fafv6DSwOndiMgeOVuTqMKFC6ugb9GiRWo1derUkT/++EOV/B08eFBV/RKRfXnx4LZsH/Oh7BjTXQV+PtnySN0xv0idET+k+sAPMKXb7E96qQendyMie2R1b98FCxaoDh8I9lasWCFt2rSRbNmySWBgoJrbt0iRIjaZb8nZi4bIHkWEh8uFv3+VY79+JWEvX4iTi6vqwVum/efi6ullj0kmInplkjNOsTr4A39/f0mbNq16js4e586dk4YNG0q+fPleWabFJTkzlcjePLlyVv775kvTXLxZS1aWKn0mSYb8xewtqURENiE54xSrgj90+Jg8ebLMnz/fJjIqIZIzU4nsRXhIkJxY/LWcWT5PtPAwcfNOK+U/HCJFmnQUJzud9gzTux29fFkdwvKFCnF6NyJ6JZIzTrEq+EPnDszv6+fnJ9mzZ38lmZNYyZmpRPYAw7bsm/mF+N+6opKDMfsqfTxGvDOnrms9oQKDg6Vkz97qY5zejYheleSMU6zq7YtevV988YWa03fKlCmq5y8RpW6hgQFy9OfJqn0feGXyldc/nSB5qjdO3QmLJycRyZ/NVy2N50RE9saqkr8bN25IqVKlJCAgwGK+TJs2TQYMGGCTeZacETVRanXnyE757+sv5cX9WyoJhRq1k4o9hot7mvSpNUlERKlScsYpVpX8oap38+bNMWYqqoSJyPaFvHguR+ZPkEvrl6idTZM9r1T5fIrkKF/T9neeiIgSxKrgD1xdXaVixYrRNnrhwgV5+vSpqhomItt1+9B2VdoX+PCO2sliLbtJuW6DxM3Lx3Z3moiIEs2q4A8dPho3biwPHkRO4m40ffp0VTI4ZsyYRO8cESUftO07/OO4/y/ty5FPqvWbKtnKVku+jaYCQSGh8vF3c9Sefv/pJ+Lp7pYK9pqIKP4SHfx9+OGHKuhDnXT79u2jbBEDPG/ZskWWLl0a/z0hohRz9/g+2TdjgLy456e2WazVB1K+2yBx9fROsX2wVRFahGw7cVLtHp4TEdmbRAd/mMUDnJ2dow3zginfunTpIi1atLC3/CJK1cKCXsqxX6bIub8WqHRgarZq/adL9tccu7TPyM3FRaZ+2E29hOdERPYm0cHfxIkT5f79+xIUFCQzZ860t3whsjuPLhyXPV99Ls9vRg5gXLhJR9WT182bQzRFCf5cXeXdWjXs7vgTEekSHfyBr6+vLF68WD1H547g4OAoOYsp37y9WY1E9CpFhIfJqd+/k5NLvlGzdGDcvqr9pkmuyvVe5W4REdEr4mztdleuXCm5c+eWjBkzqupf4+P777+3dvVEZIXnt67KpgFvy4nfpqvAL2+t5tJ87mYGfrHA9G6nr99QDzwnIrI3VpX8+fv7S+fOneXtt9+WTz/9NFrbPwSERJTyNE2TSxt+l8PzRktYUKC4+aST13uPl/z1WouTE+etiE1waKg0GzVWLcLp3YjIHlkV/GFOX3zJzJkzR1XxEtGrF/z8iRq3z2/vBrUzGLql+hdfi0/WnK9+51IBhMbZMmRQe8owmYjskVXBX7FixSR9ek77RGQr7hzZJXun9ZWXj++Ls6ubvPb+l1LynY/EydnqFh4Ow8vDQ/Z/Pc1h0ktEjseq4C88PFzGjh0rbdu2lY8++kgyZcoUJQcLFSqk2gMSUfIKDw2R479OlTPL56oNpctTWGoOmiWZCpdO3g0TEVGq46Sh3jaRbty4oUr/YjJp0iTp27dvgtZ+8eJFOX36tLz22mtSoECBOD979+5d2b9/v2TNmlWqVasW7/ZMyTlhMlFKQqeO3ZM/lccXT6jNFmnWSSr2GCmunl4puRtERJSEkjNOsarkL2/evPLy5cskSSpKEXv27KmGjilRooScOnVK+vfvL5MnT45x/d98840MHDhQihcvruYSxhzDGzduZPtDchhXNi+XA98NU5063NOkl6r9pkreGk0cJv3JAdO79f9hvlr1jI+6c3o3IrI7SdYQ6NmzZ3L48GFVGpiYwsSFCxfK33//LWfPnpUjR47IoUOHVEeSdevWWcz0vXv3yujRo+XgwYNy8uRJNc8wSgHxGpG9C335QvZM7St7p/VTgZ9vmarS7Pt/GfglAUzp9s+hw+rB6d2IyB5ZVfIHDx8+VMO8/PHHH6b8KVy4sBrjr0GDBvHOs3nz5snHH38s+fPnV58pW7asdOjQQb3erFmzaOuZPXu2KiksV66ceg/tDdevX6/2xxIMQG0chBrFqUSp0ePLp2XXxE/E/9YV1ZGjTMd+UrrdZ+LMqciSBKZ0G9upg1oXp3cjIntkdcnf+++/r0rdUFKH4ArPu3btqsb+u3w5chqp+EA7vzp16kRZtHbt2qr61xK088P7c+fOlTZt2qjAMSIiQmrUsDwtE9ofomey/siTJ098dovIZqBE/fzfv8qGvq1U4OedJYc0/GqZlO3Yl4FfEsL0bl0a1FcPPCcisjdWdfhAKRtm9MAcv+Y9fdu1aycVKlSQL7/8Ms480xs1ItArVaqUafnNmzerUj/zaePA09NTlfphIOkmTZrI7t27VbXxf//9pzqLmLNU8ocAENXV6dKli3MfiV6lkBfP5b+vB8mNXWvVbuSq0kCqD5ghHuk4kDoRkT2y2Q4fjx49Eh8fH4szeWCIlwcPHsTreCCQw6wD5kFeUFCQeHlF77GIziFhYWEq4Pznn3/UNvr06SPvvfeeavO3atWqaNv18PBQD6LU5tHFk7Jr4scScOe6Gruv/IdDpXjrD+Pds50SBjUI1+9H3rvy+WYVZ46RSER2xqpq3yJFioi7u7vMnx/ZM053+/ZtWbZsmVSvXj1e2YV1ZMmSRa5fvx5leXQesTROoIuLi+pp3LJlyyjLN2zYMMZqYqLUBtW8F9YulI39W6vAz8c3t7w5fYWUeKs7A79kFBQaKvUGD1MPPCcisjdWNWjBL+KffvpJDfK8YMECKV26tDx+/FgNt9KqVSt566234p1fb7zxhixfvjzKZ1asWKFetwRVvidORI5rpjt+/LgKColSO/Tm3f/NILm2fbVKSu6qDaXagOnikZbVvCkhrYUaByIie2FV8ActWrRQw7MsWrRILl26JAULFlSlfpZ66MYGY/rVrFlTtflDG74lS5aoTh3o7QtHjx6VXr16yc8//ywlS5aUkSNHSpUqVVTvYJT4bd26VS27cuVKezk25KCeXr8guyb0kmc3LoqTi6uU/2CwlHj7I5b2pRBvDw85+f2sFNoaEVHKszr4AwRgw4cPt2rvK1euLGvXrpURI0bIlClT1FAvCOgwbAygPeCTJ09M7QJR8rdmzRoZPHiwmmIOVdAIGJs3b27VfhC9Sle3/aVK/DB2n1fmbFJr6PfiW6ryq9wlIiKyM1b19oWbN2/KJ598okr70HEDMKVbrVq15J133rHZ7ErOXjRECYW5eQ//MFYu/P2r+mj2cjWl5uBZ4pkhS0JXRUREdiA54xSrS/4wxh56/IaEhJiCP5QEdu7cWYoWLSplypSxg0NAlHwCH96RneN7ycNzR9RGSrfvI2U79efYfa9IcGioDP3lN7X1iV07i4eb2yvaEyKi5GFV8Hf16lU1vdq9e/eijJWHkj+Mu4fSQAZ/RDG7e3yf7J70iQQ9fajm5q3+xTeSu4rlTk6UMsIjImTFnr1qY+O6dEyZjRIRpSCrgj+Mteft7R1tgGdAr9tQDpNAZBGGcTm7Yp4c/WmyaBHhkrFgSak94gdJmyMfc+wVc3VxkSFt26i9wHMiIntj1Th/6IyB2TFmzJgRJV8uXryoOl9g+jUiigqdOVDad2T+BBX4FWzQRhrN+IuBn41wd3WVnk0bqweeExHZG6vubJiV47ffflNj82GgZ32cP1T5fvjhh9K0aVN7yy8iq/jfviY7xvaQp9fOqWFcKvUaLUWbd+EwLkRElGKs7u0LmNsX7fswzh9m6qhXr57UqFEjxRKRGMnZi4bIktuHtsnuyZ9JSMAz8czoK7WHfS++pV9nZtkYTO92/9kztVe+6dNzejcieiWSM05JkuAvNUrOTCUyQvu+08vmyLFfpohommQpXkFqj5gn3pmzM6NsUGBwsJTs2Vvt2Zl5s9Wgz0REKS054xSrqn2JKHZo37dvxgC5vnOtWrBwkw5S+eOx4uLOgMKWsaMHEdkzBn9EySTgrp/sGNtdnlw5o9r3Vf5krBRt1jmZtkZJBSV9lxZETitJRGSPGPwRJQOM37drQk8Jfv5EzdJRe/hc8S1dJRm2RERElDBWDfVi7sGDB6qxNJEju7B2oWwZ2kEFfpmKlJEms9Yx8CMiIpthdfD35MkTNaxLpkyZxNfXV/0/e/ZsmTBhgs0kkiglRISFyoHvhsqB74aJFh4m+eu2kjenrRCfrDlTYvOURDC924iFi9UDz4mI7I3VwV+HDh3k1KlTKuDDbB+AwZ8nT54s27Zts7f8IrIIpXxbhnaUC2t/E3FyknLdBkuNQbPE1cOLOZbKYHq337ZuUw88JyKyN1a1+bt9+7YK8DDOH2b6+Oijj1T+NGrUSN577z3ZuHGjGvOPyJ49vX5Bto/uJgF3boirl4/UHDRLcldtaM9Jtmvo6ft5qxYqjez1S0T2yKrgLyAgQNKkSaMCP3MZMmQQF86LSXbu1sFtsntSbwkN9Jc02fNK3dELJEP+4naeavuGKd36vdXKvhNJRA7NqmpfVO+6ubnJwoULo2Sin5+fLFq0iHP7kt3CwM3n/log20d1VYEfZupo/M3fDPyIiMjmWRX8OTs7qyDvs88+k5o1a0pQUJC0a9dOSpYsqeb7RfUvkb1RHTtmDZVDc0eLFhEhhd58T96YtFQ802eyt6Q6JAT2z14EqgeeExHZmySZ3g1DvPzxxx9qbt/MmTOrdn4IBm1Zck6bQvYr2P+p7Jrwsdw9tlt17KjQfbiUeLuHODk52W+iHQyndyMiW5CccYpVbf50WbNmlU8//dQW8ooo2fjfvibbRnWT536XxNXTW2oOmS25qzRItu0RERElB6uDv8DAQPnpp5/k9OnTqtrXCD1+GzdunBz7TZSi7p/aLzvG9lBDunhnySH1xv4iGQuWTNF9oJTh5e4uF+fPVRtjb18iskdWBX9oD9OwYUM5evSoat+XPXv2KHnk6mrV6olswpUtK+W/r7+QiNAQyVSkrNQd/ZN4Z85mE/tGSQ9V+G68dxGRHbMqOkMbv71798qhQ4ekYsWKdpxN5Ijw4+bk4q/lxKIZKvl5azaV6gO/FldPDtxMRESpl1XBH8byy5gxIwM/sjvhIcHy39dfytWtK1XaSr77sZTvNlicnK3qIE+pQEhYmExbvkrt6cA2b6lx/4iI7IlV32To6PHuu+/KsmXL7ClPyMEF+z+RLcM6qsDPydlFqvSZLBU+HMrAz0GEhYfLDxs2qgeeExHZG6t+0t67d08uX76sBnkeMWKEZMoUdZwzjP+HuX+JUgv/O9dl24j35fnNy+LmnUZqDZsrOSvWSS27T0kAnTw+ahw5Rik7fBCRPbIq+PP29pYWLVqohyXFi3OaK0o9Hp4/KttGdpPgZ4/EO2vOyB69BUqkngRQkkA179B27ybJuoiIbFGSDPKcGiXn4ImU+tzYs172fNVHwoODJGOhUirw884ctfc6ERFRSknOOMWqkj/dkSNHZP78+XLx4kXVDrB+/frSrVs3cXFxSak8Ikq0c3/9JIfmjRbRNMlZub7UGjJbVfmSY0Ivb72tH6p9OXsLEdkbqzp8wK+//irVqlVTU7zVrl1bfH19ZfTo0Wrcv4iICHvLL7IjmJf38I/j5NDcUSrwK9K0k9QdvYCBn4N7GRIiRbr3Ug88JyKyN1aV/IWFhalOHejt26pVK1PejB07VipXrqzm+23fvr295RnZgfCQINk7rZ9c37lWpaZct8FSqu0nLOUhIiK7Z1Xwh2peZ2fnKIEfpEuXTlq3bq0Gf2bwR7YGQ7nsGNNd7p86IM6ublKt/3QpUP8tW9tNekUwvdvx2d+qreM5EZG9saraN0eOHOLv7y+3bt2Kli9nzpyR3Llz21t+USr34v4t2TTgHRX4uXmnlfrjf2PgR1GgjV96H2/1YHs/IrJHVgV/mOGjXbt20rhxY/nrr7/UdG8HDhyQXr16yb59+6Rt27b2mGeUSj25clY29Gslz25cFO8s2eXN6Sske7kaqTQ1REREiWNVtS/MmzdPDfDcuXNnCQgIUNXA1atXl02bNkmuXLkSt1dESezusT2yY2wPCQ30l/T5iqoSP5+sOZN4K2QPML3b7L/XqaT0btGM07sRkd1JsnH+MDwCZvzAXL8eHh42n1HJOX4O2Zar2/6SfdP7S0RYqPiWqSp1R80X9zTpbWsnyWYEBgdLyZ691f6cmTdbvFPB/YyI7E9yxilWl/zp0DYme3YOiku25cyKH+TIj+PUTuWt1VxqfDFTXNw9bWsnyaa4ODtL5/r11D7hORGRvUmy4I/IlmAMvyMLJsrZFfPUbhVv/aFU/GikOPHLnOLg4eYm47p0jGsxIqJUi8Ef2Z3w0BDZN2OgXNu2SqWt/IdDpWSbXuy5SUREJMLgj+xLaGCA7BzfU+4c2SlOzi5Srf80KdigjX0lkoiIyArJ1qAlNDQ0UbvFKeEosYKePpJ/B72nAj8XDy+pO+ZnBn6UYOjwUfjDnuqB50RE9ibZgr88efJI8+bNZdasWXLhwoU48+3atWvSoEED8fHxkUKFCsnXX38d52d0mzdvlqxZs8rw4cPj/RmyLwF3/WTjgLfk8cUT4pEuozT86g/JVTmy0T5RQoWFh6sHEZE9SrY2f99++62cPn1adu/eLf/884+sX78+xvwLCgqSWrVqyRtvvCHHjx+Xs2fPSpcuXcTLy0t69uwZa76/fPlSDSodHBwsL168sMdjRHF4cvWsbB3WSV4+vi8+vrmk/oTFkj5Pobg+RmSRp5ub/DdzqnoPz4mI7I1Vwd+JEyeke/fualYPcwmZ3WPJkiVqnMAFCxaIi4uLFC1aVMaOHSvTp0+Xjz76KNaG+lgOg0lXrFgx1mOD4BAPHcbPodQP07RtH/2BhAQ8kwz5i6nBm72z5Ej9CaNXBgPVZ8+Y8ZVtn4gouVlV7ZsvXz45evRovKp1Y7N9+3Zp0aKFCvx0rVq1kosXL1qcN1iH4PO7776TH374Qd2wYzNp0iRJnz696YFqaUrdbv73r2wZ2kEFfllLVpaGU5cz8CMiIoqDVcEf2ueNHDlS+vbtq6p3UWVrfDx48CCu7Ss3b96UvHnzRlk2d+7cqsTPz8/P4jrQMQSlgoMGDZJixYrFuZ0hQ4bIs2fPTI+Y1kupw+V//1TTtYWHBEuuKg3kjYmLxSNthtSx82TTML3bvH82qAeeExHZG6uqfW/duqWCP7DUpm/atGkyYMCAOPMM1bFuZm1rUAqIB9oDWjJnzhw1l/CXX34Z5/oBU86lhmnnKG5nls+VI/MnqAUxjEvVvl+JsyvbZlHSQEePScuWq5V1fqMe5/YlIrtjVfCHat+YgjNwdY3f6nPmzKnmBTZ69OiRhIWFqRJAcygpRM9edCRxd3e3u4NClqFd6LGfJ8vpZXPUAiXe+UgqfDiMs3ZQksKUbu/UqK7WyendiMgexS86i4VemoZSOLTRK1OmjGp/F1cbPCN8Zs2aNVFe27Ztm6RLly5adTCsXr1aVd/WqFEj2p798ssv8uTJE3s8Vg4tIjxM9n87RC5v/F3lA2btKPXuxw6dJ5Q8ML3b9B4fJM/KiYhsgFVt/gBt5xo1aqQCtQoVKqigbMqUKfLxxx+rkpr4QI9hdN749ddf1eKPHz+WMWPGSLdu3UzB5ezZs+XUqVPqfQwDg3EBjY8mTZrIBx98oNZD9iU8JEh2TfxEBX6Ym7dq36kM/IiIiBLJ6uDvvffeU+31duzYoTqAAMbrW758uaxduzZeu4Vq38WLF6sx/VCVjP9z5Mgh48ePV59HQPn555/LvHnz1P9p06ZVyxkf2DYCUPbitS+Yrm3byK7it2e9OLu5S62hc6Vw43b2lUgiIqIU5KTFt3jOguvXr0vJkiVVSR2qfxGUoTQuc+bM0rt3bzVIMzp9xBfG3jt06JAat8+8B+/58+dVYOft7W1xddgu2v8hcIwPbAtDvyCwRNBItifo2WPZNqKLPLpwXFy9fKTuqAWSvVz0qn6ipIQp3ar0HahWuf/raeLt4ZGUqyciipfkjFOsavMXHh6uSv0sdbpAT10EgwmBIKx+/foWPxLXcC758+dPyKbIxr14cEe2Dusoz25cVNO11Ru3ULIUK2fje032wv/lS3tJChFRNFZV+xYsWFCyZMkikydPjrLiY8eOycKFC9VcvUQJ9fzWVdk04C0V+GG2jjenrWDgRykGU7ptmzxBPTi9GxHZI6uCP1i2bJnMnTtXChQoIIGBgVK7dm2pWrWqDBw40GJvXKLYPL58WjYNeFte3L8laXMVkDenr5T0eYsw0yjFYKSCAtmzqUdCRi0gIkotrKr2BfTwRXu8DRs2yKVLl1R7v1q1aknhwoVTSx6QjcA8vdtGdZPQF88lY6FS8saEReKZIYuN7B0REZF9sCr4w3h6f//9tzRt2lRat25tHzlCr8Stg9tk5/iPJDw4SHxLvy51x/ws7j7siEMpLzQsTJZu36k23L5ubXGL52D1RESphVXVvpiWDXPrZsuWTZX2TZ06Vc6dO5da0k424tqONbJ99Acq8MtZub7UH7+IgR+9MqHh4TJy0RL1wHMiIntjVfCXPXt2uX37thw+fFiaN2+uqn7LlSsnRYoUUXP6Xrhwwd7yi5LYxfVLZPfkT0ULD5N8dVpKnZE/iqunVxJvhSj+nJ2cpWmliuqB50RE9saqcf7MoSQQ07J98803KhDEGH8IAm1Rco6fQ/Fz+s+5cnTBBLVwkWadpPIn48XZxSV+HyYiIrJjyRmnWNXmDzCd2qZNm9Rj9+7datBnTLU2YsQIqVKlih0fFkosjCt+7JcpcvqP2WoVpdp+IuW6DRYnJ6fErpKIiIjiyargDzN8vPbaa1KoUCEZPny4LF26VPX2JYqJFhEhB+cMlwtrf1OLlP9giAr+iIiIKGVY1aDF19dXVe0i+OvTp4+0bNlSRo8eLQcOHJCIiIiUSQGlGhFhobJnat/IwM/JSV7/bBIDP7I5L/83vRumeMNzIiJ7Y1Xwh7l70aZv48aNcv/+fRX4BQUFSb9+/VRnkBUrVthbflEihYcEyc7xPeXatlXi5OIqNb78Voo265TItRElH01E7j19qh54TkRkb6yq9jW6cuWKnDp1Sk6ePClnzpyRHDlySJo0aewtvygRQgMDZPuYD+Xe8b3i4u4htYbNldxVOPUf2SYPNzdZN2ak2jk8JyKyN1YFf+iJgpI/9Oy9c+eO1KxZU1X9orcvZ/ggCPZ/IttGvC8Pzx0VVy8fqTv6Z8n+WjVmDtksF2dnKZUvr83uHxGRtawK/gICAuTp06cyceJEadasmWTKlMna/SE7Evjonmwd1lGeXjsv7mkzSP3xv0mWYuXsKIVERESpT5KN8/fy5UtV9Zs1a1bVEcTWJef4OSQScNdPNg9pLwF3rotXJl95Y+ISyZC/GLOGbB6md/tr3361n62rVeH0bkT0SiRnnGJVhw/A/L5du3ZV7ftKly6tpnrD3507I+fGJMfzzO+SbBr4jgr8fLLlkTenr2TgR6kGpnT7YsHP6sHp3YjIHlkd/HXr1k1N44ZgD1Hq1atX5a233pIWLVqo5+RYHl86JZsGtpHAh3ckfd4i0mj6CkmbI59jZQKlapjSrV7ZMurB6d2IyB5ZVe37+PFjVc177949yZIlS5T8effdd+X111+XL774wibzLTmLUx3V/dMHZdvIrhL64rlkKlJG6o9fJJ7p2Q6UiIgooZIzTrGq5A9j+6G619KsHvnz51fvk2O4c2SnbBnaUQV+vqVflwaTf2fgR0REZIOsCv6KFCkirq6u8ttvkVN16RD0LV++nHP7Oogbu9fLtlHdJDz4peSsVFeV+Ln7pHOQ1BMREaUuVg314uLiIj/++KO0b99efvrpJ9XRA1XB69atkzfeeEPeeeed1JUblGCX//1T/ps5UM3Zm7dWc6nx5Tfi4uae4PUQ2QpM6dZk5Bi1O+vHjhIvDw9b2TUioiRhVckfvP3222pmDwzw/PDhQzXWHwLBlStXipOTU5LsJNmmc6t/ln3T+6vAr1CjdlJz8HcM/CjVw5Ru1+7dVw9O70ZE9siqDh+pWXI2pLR3mqbJqaXfyvGF01RSi7/VXSp+NJLBPtmF8IgIOXr5skpL+UKF1IwfREQpLTnjlETf1fbt2yeNGjWSfPnySe3atVVVL9k/BH5H5k8wBX5lO/Vn4Ed2BcFepSJF1IOBHxHZo0QFf5jJA236MMzL8OHDpWzZsmpsv927d9tjHtH/RISHy/5vBsnZFfPUK5V6jZaynfqxxI+IiCgVSVSHj8WLF6s2fosWLTIlFe375s+fr14n+xMeGiJ7p/WV6zv+FidnZ6na9ysp9OZ79pdQcnhh4eGy8fBRlQ+NKpYXVxcXh88TIrIviSr5u3jxolStWjVKTlSrVk3N9EH2JyzopewY20MFfs6ublJzyBwGfmS3QsLCpPecueqB50RE9iZRJX+hoaHiYTb8Af7H62RfQl48l+2jusn9UwfExcNT6oz4QXJWqmdfiSQycHZykirFiqpX8JyIyN4kKviDM2fOyOrVq035ceDAAXny5EmU10qVKiWFCxe2tzxzGEFPH8nW4Z3l8aWT4uadVuqN/UXN3kFkzzzd3eWPIV/acxKJyMElKvhzdnaWZcuWqYe5Nm3amF6aMmWK9O/f38GzOHV68eCObBnaQZ77XRKP9JnljQmLJFPh0qkzMURERGTCcf6ePZN06TgVmdHzW1dly5AO8uL+TfHOmlPemLhE0ucpxMuGiIgohSTnOH+JKvkj+/XkyhnZMrSjBD19KGlzFZAGk5aKj28u+00wkZmgkBB5a9wk9eqqEUNUNTARkT1h8Ecm908flG0ju0roi+eSsVApVdXrmSELc4gcSoSmyVk/P5VmPCcisjcM/ki5fWib7Bj3kYQHB0nWUpVV5w53H1aHk+PxcHOT3wb2UwnHcyIie8Pgj+TajjWyd2pfiQgLlZyV60ntYfPE1dOLOUMOCVO61SpdyiHTTkSOgcGfg7uwdqEcmD1cRNMkX52WUn3gTHFxYxsnIiIie8Xgz0Fpmianln4rxxdOUzlQtHkXqfzJODV1G5Ejw/RuO0+eVllQu0wpTu9GRHaHwZ8D0iIi5PCPY+XcqgUq9WU6fC5lOw9Q8zMTOTpM6fbB19+qbDgzbzaDPyKyOwz+HAza9e2b+YVc3bJCpbxSr9FSvPWHDpYLRDHDlG5lC+RXC3B6NyKyRzYV/J04cUIWLVokt2/fltq1a0v37t3VbCIxwXLffvutXLp0SYoVKyZ9+vSRbNmy2eNxShJhQS9l18SP5daBLeLk7CLV+k+Tgg3+f0YWIhI1rt+aUcOZFURkt2ymgdfevXulRo0a8uDBAylRooRMnjxZOnbsGGPG37p1S8qWLStXrlyRKlWqyK5du6RcuXJy//59uz1Y1gj2f6oGb0bg5+LuIXVGzWfgR0RE5IBsZnq3+vXrS6lSpWTWrFnqMFy7dk2V5u3YsUOqVq0a7dCMGjVKtmzZIrt371bvoQNDxYoVVcA4YMCAOA9lck6bYmsCH92TrcM6ydNr58Q9TXqpO/on8S39uq3tJhEREf1PcsYpNlHti9K6bdu2mQI/yJ8/vzRs2FB+//13i8FfmzZtpG3btqaTBJ0VvLy85MmTJxZPnODgYPUwZqojeH7zimwZ1kle3PMTr0y+Un/8IslYsIQjJJ0oUTC9W8evpqvPLv5yAKd3IyK7YxPB3/Xr11XbvuLFi0fJ4JIlS8rFixctZnqZMmWivH7s2DHZv39/lADSaNKkSTJmzBi7O4CxeXj+mGwb+b4EP3ssaXPml/oTFknaHPkcKg+IEgpTuh2+dFl9jNO7EZE9song7+HDh5ImTRpxcXGJkscZMmRQbQDjgpLDVq1aqereChUqWFx8yJAh0r9/f9N7KE7NkydPXKtOte4c2Sk7xvaQsKBAyVSkjNQft5Dz9BLFg7urq8z7rLdaEs+JiOyNTdzZsmbNKv7+/hIWFiauhpvt48ePxdfXN9Y8DwgIkGbNmqmgD6V7MfHw8FAPR3B121+yb3p/NaxL9nI1pc7IH8XNO40jJJ3Iaq4uLtKoYnmr10NEZKtsordvvnyRVZGnTp2Kkk8nT540vWcJ2vC1bt1alRouXbo01mFhHMXZlT/KnimfqcAvX50WUm/sLwz8iIiIyMQmgj+U/KFzB8bs0505c0a2bt1qGu4FpYJ+fn6m9yMiItR7qL5ds2aNeHp6OvRhjZy1Y5wc/mGsyodirT6QmoO+U8O6EFH8hUdEyL6z59QDz4mI7I1NVPvCxIkTVQCI9n+FCxeWJUuWyPvvvy+VKlVS7w8dOlRmzJghd+/elSxZssjw4cNlxYoVquTP2JYPw8MMHDjQ3o5TrMJDQ2TfjAFybdtfarnyHw6Tkm16cro2okQIDg2V9lMi57zG9G7eDtJchIgch80EfxijD711MbQLAryZM2dK+/btTUeiQ4cOqgNI5syZ1WulS5dWAaC5TJkyOc7RE5GQF/6ya0JPuXNklzi5uEq1fpi14x2HygOipIQZrovkzKlWydmuicge2cwgzyktOQdPTCmBj+7KthHvy5MrZ8TV01tqD58nOSvVTanNExERUTJJzjjFZkr+KGGeXjsvW0d0kcAHt8UzY1bVsSNzkbIJWwkRERE5HAZ/qdC9E/vUGH4hAc8kXe5CUn/8QkmTPW8qTAkRERGlNAZ/qczVbatk34yBEhEaIllLVlLz9Hqky5jKUkFkuzC9W/dvIucYn//5Z5zejYjsDoO/VELTNDn9x3dy7Jev1B7nqdFEanz5jbh6eKWSFBClDpjSbffps2pnOb0bEdkjBn+pAAZs3j9rqFze+Lva2xJv91DDuTibTYdHRNbDlG5ff9RdrYjTuxGRPWLwZ+Mih3LppebqdXJ2lkq9xkixll1tfK+JUi9M79a6etXUmwAiojgw+LNhAfduyvZRXVXPXhcPL6k1ZLbkrtrQhveYiIiIbB2DPxv18NxR2T7mQwl68kA8M/pKvTE/Seair9no3hLZD0zpduradZWg0vnziYuzTcyCSUSUZBj82aDru9bJ3qmfS3hIsGQoUEKN4eeTNXLGASJKXpjerdXYCWojnN6NiOwRgz8bYt6jN9frb0jNwd+Jm3caG9pLIvuGKd1y/W8aSU7vRkT2iMGfjQgPCZJ9M7+Ua9tWqT0q1uoDqfjRSPboJUphXh4esmf6lBTeKhFRymHwZwNePr6vZux4eO6IODm7SOXe46Ros842sGdERERkbxj8vWKPL5+W7aM/UHP0uqdJL7WHz5Ps5Wq84r0iIiIie8Xg7xW6vvNv2Tt9gIQHv1Rz9NYd87Oky1XgFe4REQWFhMpn389TGTHr457i6e7GTCEiu8Lg7xWICA+X4wunyuk/Zqut56hYR43hh5I/Inq1IrQI+ffoMbUTeE5EZG8Y/KWwkBfPZc+UPnLrwBa15ZJtekm5boPZsYPIRri5uMikrl3U3uA5EZG9YfCXgp7duKg6djy/eVlc3D2kat+pUqD+Wym4B0QUFzdXV2lft3ZcixERpVoM/lLIjd3rZe/0fhL28oV4Z8khdUb+yBk7iIiIKMUx+Etmqn3fr1Pl9LLI9n2+ZapKraFzxCtj1mTeMhElRkREhFy6c0d9tHCOHOLM6d2IyM4w+EtGwc+fyO7JveXOkV1qKyXe7iHlPxwqzi7MdiJbFRQaKm8OG6V2j9O7EZE9YhSSjE4smqECPxcPL6nWb6rkr9sqGbdGREklU1pOqUhE9ovBXzIq13WQBNy7KeW6fikZC5RIxi0RUVLx9vCQI7O+TqrVERHZHAZ/ycjNO43UG/NzMm6BiIiIKGGcE7Y4EREREaVmLPkjIjLA9G6DfvpFvTLlg66c3o2I7A5L/oiIDDCl2+r/9qsHp3cjInvEkj8iIgNM6Tai/XvqFU7vRkT2iMEfEZEx+HN1lQ8bNWSeEJHdYrUvERERkQNhyR8RkQGmd7v16LF6JVfmTJzejYjsDoM/IiIDTO9W64vB6hVO70ZE9ojBHxGRGS93d+YJEdktBn9ERAaY3u3sD3OYJ0Rkt9jhg4iIiMiBMPgjIiIiciCs9iUiMggODZVRvy1Rr4zp3EE83NyYP0RkV1jyR0RkEB4RIb/v3KUeeE5EZG9Y8kdEZODq4iID326tXsFzIiJ7w+CPiMjA3dVVPm3ZnHlCRHbLpqp9w8PD5fDhw7Jy5Uq5detWnJmOkfiPHTsmK1asED8/vziXJyIiInJ0NhP8+fv7y5tvvikNGzaUCRMmSIECBWTu3LkxHp/AwEBp0aKF1K1bVyZNmiSFChWSGTNmOPrxJCIraZomj577qweeExHZG5up9h03bpw8evRIrl+/LmnTppUNGzZI8+bNpX79+lK0aNFo+T5t2jS5fPmyXL16VTJmzCg7duxQgWO9evWkfPny9naciCiFvAwJkYp9+qmtcXo3IrJHNhH8obp3wYIF8v3336vADxo3biy1atWSH3/8UaZOnRot7+fNm6cCRgR+UKdOHWnUqJH88MMPaj3mgoOD1UP37Nkz9fT58+f2eFyJKJECg4MlIiTEdH8I8/BI5JqIiBJPj0+SowbCJoK/u3fvyuPHj1WwZ4T/Dx06FC3nELjdvn3b4vJr1661mNOoGh4zZky09/LkyZP4I0NEdi3H0t/sOn1EZPvQLC59+vRJuqM2EfzpnTuyZs0aJXG+vr5y8+bNaAnWl8f7RjEtD0OGDJH+/fubFkdnEQScmTNnFicnpyTN1NQGvy4QBKPTTLp06VLb7icL5kl0zJOomB88R+KD5wnPk/iwdJ6gxA+BX86cOeOzigSxieDPy8tL7TSqZV1d/3+XgoKCxMfHJ1qC9OXxvjEajml58PDwUA+jDBkyJCiz7B0CPwZ/UTFPomOe8ByJC88RXjfxwfMk7vMkqUv8dDbR2zd37txqf9DZwwj/58qVK1ruIApGaV18lyciIiKiSDYR/KHTRoUKFeTPP/80HZewsDBZvXq16sFrDiV4aA9oXB7VuKtWrbK4PBERERFFsolqXxg8eLB07txZcuTIIWXLlpXp06er+u5OnTqp9//66y+ZOHGiGgImU6ZMMmjQIHnrrbckb9688vrrr8vs2bPlyZMn8sEHH/DYJhCC6VGjRkWrFndkzJPomCdRMT94jsQHzxOeJ/GR0ueJk5YcfYgT6bfffpOvv/5abty4IVWrVlUBHYI7QCkfhnbZunWrZMmSxfQahoG5cuWKVK5cWWbNmiWFCxdO5NaJiIiI7J9NBX9ERERElLxsos0fEREREaUMBn9EREREDoTBHxEREZEDYfBHRERE5EAY/BERERE5EAZ/durQoUNSvnz5aI/79++bUvzs2TPp3r27Gh6nZs2aaqgdc2fPnpXWrVtL/vz5pWnTpvLff/+lqhzD/M0dOnRQ+XDs2DHTvgcGBkrv3r2lWLFiUq1aNZk7d260dF29elXatGmj0t6oUSPZtm1btGW2bNkib775phQoUEDeffdduXbtms3nD8bDxPiZyJPDhw+b9rdnz57RzheM/2iE8Tbr1asnBQsWVOu4d+9etPR+9dVXUrFiRSlVqpQajzM0NNSm8+TXX39VY4ViWCmMHXru3Lko+4v9xzikSA/SNWXKlGjpwXWFcUqRL8gf5JO5gwcPSvPmzdX51LJlSzl16pTN5suiRYtUnmDOb1z/Z86cibKvw4YNi3au4HoywvWC6wbpxXWE68ncvHnz1PWH6xCfx3Vpi3DN9O3bV0qUKKHul7169VL3FqO7d+9Kx44d1TlQv359WbNmTbSk7N+/X5o1a6bypFWrVtHyFXAfxiQG2A7uz7hP26L45EmDBg2inScLFy6MkpzvvvtOqlSpIsWLF5fPP/9cTdNqhGlf+/Xrp/Ie5ySGg0sN9O+d8PDwKLt7/vx5dZ/Rv1P37t0bLTnr1q2TN954Q51L7du3l9u3b0dbZubMmWqIu5IlS8qAAQNUPiUIxvkj+7No0SItd+7c2vLly6M8AgMDVWIjIiK0SpUqafXq1dP+/fdf7YcfftC8vb21n3/+2ZQZfn5+WqZMmbSePXtqO3bs0AYPHqyWOX78eKrIsIMHD2p58uTRmjRpoomItnv3btN+169fX6tSpYq2ceNG7ddff9XSp0+vzZw50/T+o0ePtBw5cmidO3fWtm/fro0dO1bz8PDQ9uzZY1oG68Nr48ePV8t06NBBfebx48c2mz+HDx/W8uXLZ8oT7LeuePHi2ueffx7lfEEe6v7++2+VXuTT1q1btRYtWmiFCxfWXr58aVpm6NChWvbs2bVly5ap5UuUKKG99957NpsfixcvVsce5wCuA0vHEK8hHUgP0oVjPGjQIFOagoKCtCJFiqj8QL58/fXXKp9Wr15tWubMmTNamjRptC+++EJdS5988omWIUMG7erVqzaXN3/88YeWLl06dS/YvHmz1qlTJy1btmzqmtDVrVtXe//996OcK7t27TK9v2/fPs3T01MbM2aMOsdwHeG8ePDggWmZb7/91pT3uA5xPWK9tqh27dpa9erVtd9//12luVy5clrDhg1Nu4proFChQlrr1q21bdu2adOnT1fnwLp160zLnDx5UvPx8VHnDs4B3FczZsyo3bhxw7QMzkPcY3E/xvmI+3OFChW08PBwm8uWuPIkJCREc3Z2Vnlh/B46f/68KS1fffWV+o7B99X69eu1ihUrao0aNYqS1pYtW6r8/ueff7QlS5ZomTNn1saNG2dz+WGEa8fd3V1974SFhZneun37tpYlSxatR48e6hzA/dLLy0s7cuSIaZlNmzapc2fq1KnqXHr77bfVPTsgIMC0DL6PsmbNqi1dulSdY2XLltVatWqVoDxh8GenJkyYoL3xxhsxpg4nTNq0abXnz5+blsGXFr7MERjCkCFDtBo1akRZR5s2bdRNPzUYOXKkShNuzMbgDwGcm5ubdu/ePVMyfvnlF/WljhsWTJkyRStTpkyUm2737t2jXGD4sscNXIcbdMmSJdVFa6vwZYybMdKJPDEGf/jSwc0mJtWqVdOGDRtmejs4OFjLmTOn6QfDs2fPVICDm7ju7NmzmpOTU5Qbvi3BF40xTTj3cR7gywouXryovsAQvOk2bNigvsSfPHmiXsIXNvIB+aEbPny4CmZ0uNmb35wR6AwYMMCWskNp1qxZlOAWeYIfUQh8dAUKFFBfPDF56623tA8//DDKtYEvqEmTJqnXQkNDo5w7gOsRX5jGINIWIGBt0KCB+jGsQyCC6wfnPCxYsEDlkX7/APxYNt4/u3Xrpr3zzjtRklSzZk1TXiOfixYtGuVHKO7PuE/jh4ctiU+eXL58WeWRMeA3wo8mPYDR3bx5U3NxcTH96ERBA66/a9eumZbBtYmg+cWLF7aUJSb37983Bajmwd+IESO0qlWrRtnvdu3aaR07djS9hoDfeF/AtZI/f37t+++/V8sg3fjhuGrVKtNnrly5ovLpxIkT8c4TBn92CkEJbjYXLlzQtmzZok5Io08//VRr27ZtlNdu3bqlTlb9lxl+beHXuRFKBXDBpga4uQD+GoM/fDG/+eabUZLg7++vLp79+/er1/HFjF9XRiidwJcTbtJ44DlKeoxwcePitVV6nuDXqDH4ww0a/6N0Ar9Ijx49GuUX69OnT9X75jeXjz/+WHv33XdVcteuXatuysZfuoAg+ptvvrHJLMGXqzFow3mAEqs1a9ao/f3uu+9U0GKEIB83X71kDzdvlOQZ4UsL+aWXlqEUHqUWRrNnz1YlirYGeWDME3zZ4IfBypUr1a4i/fjxhB8KKP0+cOBAlOUBpRkowTDCD49atWqplw4dOqS+5LEtI5Sm4UenrUMJMI6p/kMZP4r79OkTZbdRym4MEPUScSNcF7g+4NKlS2p5lA4Z4frCdWbrzPME5wfOA1wD+A46d+5clCToNSfGmgPAOYJzBSZPnqxKXI1wbuIaxTptEUrKUVuAc9w8+ENt24wZM6Ls9ooVK9R9E/SCClxTRv369VOFDYDSePzINv7QgMqVK6tCi/iymbl9KWlhijy0+1u7dq2aMxBtkjAVXp8+fdSGbt68qdrZGOXMmVMte+vWLSlatKhaBm3ZjPD/w4cPVfsCLGvLYto/S+lKkyaNmjYQadfzx1LaQ0JCVPojIiLUc0vL/PHHHzabLTHlCc4XQBstLy8v1cYEbXBWr16t2n0hz8BSeo8fP27KM7Sbc3FxiXbO6Plqa9KmTRtll6ZNmybZs2dX+RDTeeDs7KzSaTxX0BbQ/DoBLJM+fXq5c+eOxbzT89WW4FowmjFjhro2mjRpol5G2za0g0RbJG9vb3Vvwb0D7RzRLuvRo0fy8uVLi+k15hnWab4tnF+2eq4Aphzdt2+f7NmzRz13cnIynSc1atSI8RzANYX2sZauH2Oe4PrE/Pbm55J5O1RbElOe4J6C+yTOiXTp0qn20LiucH/08fFR5z6uNU9PzxjvF5bu1e7u7up8s8Xz5N9//5V//vlHHS/9nmoU03cq2k+ivSvuEzHdZ3fs2GE613LlyiVubm5W3WfZ4cNOjR8/Xs2HjBuzn5+f/PTTT9K/f385evSoSnFAQIC6cZvDTer58+cxLoPPYDpofD61iintSFtcaQcso6ffPA+N60hNChUqpH4ooAPIpUuX1I0V54LeiF9PL14ziivPILXkyYoVK1RnjiVLlqgvmNiuk7jSrecTzhU0YEejb0vnyosXL9QXpK1CQDdhwgRZvHix6Us6Y8aMsn79ehUEXb58WQWDRYoUkffff9+UZ/pxt7dzBWm+fv26Cl6MP3Is3VP0441zAF/suG9aOgf8/f3jfa7ZopjyBJ09cE3hfnLx4kXV0eHkyZPqfLLm2rLVewp+8KAj0PTp0yVr1qwWD1Vs5wnOg/h8r8Tn+ys+GPzZqQoVKkidOnVMqUMvNJRObNiwQb2GX5fGnr+AkqynT5+qXxUxLYPP4Ndp5syZU23OWUo7IG1xpR3wq1P/dW6eh8Z1pCYonUIvRPwSB3zBjx49WjZu3ChhYWGm9D548CBKsuLKM/N8tVX4Vd2lSxfV8xe9T3WJPVf0fML5hEABJYyWzhWsHyWJtmjXrl3qvrFgwQI1GoAOgW3jxo1Vb0VA+iZNmqR6MyON+jkU27URn3y1RegFjXxBnqDHpt5T01J69PTjHMD1hXyzlCe4n+j3W9x/zXvH2/o9JaY8QbpwT9FL19FzFT2m0ZM1pu+X+NyHbfWeMmbMGHVNdO3aNcZTN6bzBKV4CBjj872SVNeObd51yGq4KZnfaPCrU7+xoErPvIs5fsG5urqqixQsLYPPmFcXpzaW0oWqS/xyQ3V3bGlHdR9u4shLVIea52FqzR9UN2AYCiOkEaVSKLXCjRw38djSi3zFcB561QWgeQBKE235nEFpOIZewS/2tm3bRskDpAn7bxxGAVXiqMLS0xTTdYIqTf2LPbVdSydOnFB5gpJQDFlhhGpdBEDm5wrghwJ+HKJ6M7ZzBdcZSsOwHSPcg2ztXEHJJqp5jTA0CYZiQcAb2/0C1Z16MBzXOYD7Lu6/yAMjW7ynxCdPUOJnPpyR8TsIxxnnEkoEdbjfYDgx47WF+5JxuBSsF99ttnaezJo1S44cOaK+F/DA0GiQL18+Wbp0aaznCUrO8SMQAWCmTJnivM/q9yAdrrsDBw4kLE/i2ziQUhd0w0evPfQUAgy9gEba+nAn6NyB/9HYFNB4FL2D0XhdhwbeaNiudwDBZ3LlyqUawacm5h0+0EsTDdjRQw/QgB1DNDRu3NiULHTkwDLo+AAPHz5UPaExrIsOHULQO08f2gUNfNEQ2diD1laZd/jYuXOnapyMv4DzBo3YjZ1X0KAdDZb1XnbIIzTaN3YCwbAU6NmqQ69zX19fU0cTW4PevNg/NCy3BI3L0VAf6dB99NFHqjOU7vTp0yof9M4/yB/kEzpV6TB0B3oR6z1GsV0McWHsQWsr0EvTPM1GuB/gXPnzzz/Vy+j0gA4Jxs4rEydOVNcLrhvAdYTrCY3VdU2bNlU9oPVhTHA96h0EbAk6gSG9xs5duI/i/qlfL+jg4+rqqjpLATqylC9fXuvbt68pKXPmzFH3T9xHAfmI+6vesxzQUQDDUOmN+XF/xnqNPY1tQXzyBOlFxwR9KBt0rkKefPbZZ6YkIK3oeKh3nEEnKOMoFHrPVr3DGJZDHqGXtK3BKBI4L/THvHnzTJ3q7t69q3YXHckwjJI+esCdO3e0vHnzqlEpdF9++aXqZKaPxKGfa8Zht9AJpkuXLqbPYASHhPaAZvBnp3CzxZcaxkTCcAMYNwg9UY3mz5+vXkfXc/RcwzAlxpsMbsr4osOwFggmMR4ahnDQe4ymFubBH6DXMgK1119/XXWjx9AV+EI2GjhwoPoyQu8zfFGjJ6KxdyKeYwgI5DGWwfpw4aYG5sEf4Ascw7JgrC2MK1WwYMEoQ7TgixzDl2DMN5xT6O1sHjThvMP7xYoVUwESbv7oJW2rEKShl3fp0qWjPIzBHsZbwxcS0oN0IX3GcbkAw/sgP5AveB/nlXGICwSR6LWJGz+uJeQLeuOb94y2BbgfIE9KlSoVJU+Mvd9Hjx6tlsGXFMY4RGBr7KGIMckwjA6uG/3a6N+/f5TkIcjEOYbrD/mF6ye24WNeJfRARnpx/JFePDfvlYxhbHA/RWCCcwD5aBwbEfchjNmG+yjOAdxXMXyUcQw/DHWCfMf9GJ/HOYUfDrYorjzBOY/7I4JXBCsYJQJpwsgBOvQARvCDMRJxLeJ+qxdI6BAw4YcD7kv4QYG8wQ8uW/fff/9F6+2r/1BCenAOIPjHOIbGHs8onMA5hO9v/MU5pfd+1p06dUr9kMC5iB/cOJcSOhyQE1aW2qrtKH7Q0Hjz5s3qL6op0OvK3JUrV1SVRrZs2VS7HkudQFDtheoZfYR1W22jFBNUJaCBevXq1VVbNh06wqD6Cq8h7eY9PwHpRvpRLF+1alVVLWOE4na94TvaVJYpUyZ+B+cVw2WPXmlo34ZqBh1mQUF6UVWFmSrMzwdUf6IaCj3Z9FH5zaHn2s6dO9WytWvXNlV72SLMQmFpVgm03cGMHjpUcyFNqNLE7AvGPNOh+gpVVmgagHPNvGc1zkNUE2K2DJwnlSpVssUske3bt6uOKOaQLuP5jdl/kF5cPzhX0K7N/NrA+2jsj2ujbNmy0daJRu67d+9W5wyuQWzDVqG6EenBtYN7Ae4J5tDLE9WUqOrDtWXpHED1HPIO+WHeSxxwPiJP0DsYPYj1Zji2KD55gusG14Y++4359we+n3Afxl9cW7lz546WVHRAwzJoSoFlzM81W4T2m/huRfWv3itch6YmuNeiKQBmLTEfIQHt7/G9guYleN/SdzfWjzzBtYr7rN7EJL4Y/BERERE5EHb4ICIiInIgDP6IiIiIHAiDPyIiIiIHwuCPiIiIyIEw+CMiIiJyIAz+iIiIiBwIgz8iIiIiB8Lgj4iIiMiBMPgjIiIiciAM/oiIiIgcSNSJSonIIXz77bcSGhoaLa2dOnVS8zwnhzt37siSJUtkwIABavULFixQc+BamrcyPjC/8F9//SV9+vSJz+JERPQ/LPkjckDDhg2T33//XU0gb3wEBAQkW24gWBs4cKBp/ZMnT5YjR44kensXL15U6XiVMKl6w4YNxc/PL1l3Y/HixfLll18m6zbigvPFePyIKPViyR+Rg8IX+XvvvWd16iMiIsTJyUk97JWmaerh7Bz193JYWJhs3rxZAgMDLSY9PDxcXFxcrM6Wa9euydGjR2NcD46Bcd+wX66uSXt7v379uhw+fNjqtEBM+WmeX8mRDiISlvwRUXStW7eWuXPnyvvvvy8ZM2aUAgUKyPLly00Ljho1Snr37i09e/aUdOnSyT///KPeW79+vVSsWFG8vb2lbNmy8ueff8Y7e7ds2SKVK1cWLy8vKVmyZLTPIsgqV66cer9WrVpy69atWNedNWtWtV9VqlRRn8G69+7dG+UzsW3z119/lebNm8u4ceMkS5YsMmvWrCifPX78uBQsWFC9VrVqVVMgjYD6q6++kmbNmqntImiCOXPmSPHixVXeYP8PHTpkWl9wcLAMGTJESpcuLe7u7lKhQgX577//1PufffaZTJw4UXbs2KHStGnTJrl69ap6vnHjRilSpIh6jvWfP39eqlWrJp6enmodly9fNm0jKChIVZHnypVLMmTIIO3bt5e7d++a3sfyOMYtW7aUtGnTqmOAbUK/fv1k7NixsmfPHrUtHGdLkGc1atRQ6c6XL5/8+OOPpsVQzY/zCudO5syZ5fvvv1fv5c6dW1auXKnOG5xL/fv3lydPnkiTJk3Ex8dH8ufPr0qniSgJaUTkcNKkSaMtXrxYCwsLi/LQ1a1bVytcuLC2fPly7dSpU1rXrl01b29v7eHDh2qRgQMHqv/Hjh2rnThxQnv69Km2Z88ezd3dXZs1a5Z26dIl7aefftK8vLy0devWqc/8999/WmShTySsf9GiRer5sWPH1PoWLFig3bhxQ72Ode3atUu9f/78ec3NzU0bP368WvfSpUu13Llzq3TEBNvCNlatWqVdvHhRGzp0qNqfmzdvxmubc+fOVfvUq1cv7ejRo9q9e/eibCo4OFilHdvZsmWL6f3GjRtrmTJl0ubPn6+dPHlSLffjjz9quXLl0jZs2KBdv35d7UuWLFm0+/fvq3VOmTJFvb927Vrt8uXLaps5cuRQ7yHPBwwYoNWoUUO7du2aFhgYqPIA2+3Ro4faBvLD09NTq1atmspv7Ff58uW1d99917TPHTt21GrVqqXSgvxs0qSJVqdOHS0iIkItkz9/fvWZTZs2aUeOHNGaNWum5cyZU50f2IdBgwZpVatWNe2DuYMHD6pjNnv2bO3q1avqWDo5OanzB7777juVn71791b7oKcd+920aVNt37592rZt27QMGTJoZcqU0RYuXKhduHBBa9WqlVaiRAmHu0aJkhODPyIHhKAJwYPxgS9/HYI/fEnr8OXv6uqqghxA8FezZs0oOYcv6e7du0d5bciQISrggNiCv06dOmmfffZZlM+2bNlS69Kli3qtT58+WsOGDaO8jwAqruBv6tSpUT5TpUoVbfTo0eq1uLaJ4A8BmB4cWYKgF3l47tw509sI/hC8GSEIRbCtQ0CFdc+ZM8f0mr+/v+n5oUOH1HoR7AGC3gYNGpje14M/BEe6SpUqaX379jX9j3UXLVpU/e/n56c5OzubAl84ffq02saZM2fU/zj+xvw6fvy4eh/BKEyePFmdF7ExpiEkJETz9fVVQTAg+MuTJ0+0/ETwt3LlStNq8UPDmNadO3eqIDIoKCjWbRNR/CVtoxAiSjVQjdmiRQvT/pq3EUN1nA7VdKj69Pf3N72WPXv2KGk9ffp0tA4YqN5E9XFcTp06pdrN1atXz7QoOohkypRJ/X/u3DlVfWuEaty4mC+D6tmzZ8+qj8W1TUDP58S0ZTTmzcuXL1X168yZM6NUg6K39ZkzZ9R28BxVqajaRts+vbo2rg44ev5AmjRpomwX74WEhKj3cGzQJhC9uY1wzLEPeo9r4zEvVqyYWtR4zOOCDjxIB6rX0SHn8ePHUdIQU37GlQ60EUT7Pw8Pj7h2gYjigcEfkYNC2y20y4pJQoMefDGbDx+D4CM+X9ho54a2aq1atYqyO2jzBWjDhnZxRub/W6IHPzq0e8O6IK5tJhUEWOj00blzZ1NApUP7O+jRo4ds3bpVDYPTrVs3yZkzZ6zHJqGQVuzHoEGDoh3XMmXKmFaX0GNutGjRIvnggw9Ur2S060P7RrQfJCLbw+CPiJJEpUqVZN26dSqQ0a1duzZeJXT47O3bt6VRo0amz6KzgT7mIDp6rFq1SqZMmWIKUNDxIS4bNmxQQ7EASo7wGX1cwLi2mRBYd0zc3NxU5xeUABrTt2bNGlP6/v33X1Vq+sknn6jV7Nu3L9rqYttGXBDgoUcwOu/oJagoVURJo7HULS6x7QPSgA4y48ePV6tBoI1eykRkexj8ETko9LC8dOlSlNR37NhR9a5MjJEjR0r58uVV8Ieemjt37lSDOpv3sLVk6NChKkDq3r276iWLqkj0Lp03b57qzYoer7Nnz5Z3331X9VLdv3+/HDx4MM7dRPCH0kOs+6efflIlk9gGxLXN+EifPr2qap02bZp06NDBFGiamzRpkrzzzjuqahuBLPYLPWORDpTAIiD74YcfVCCGoMnYSxZQBYvAF/nRuHHj+OxaFKiyHzx4sNoHBJkIArEN9Bo2BqSxwT4gsMM+4DPm5wnSgB7L6MWLPJk/f36C95OIUgYHeSZyQAgC0AYM7d6MD72NWf369aPNvIEqWb2aEgEMhvQwwrAn+qDNaN/27NkzFaC99tprajG0G2zXrp3pIwjy9GpnrBdDp6CKeMaMGWqYk6VLl0rXrl3V8gheMDQKgi2sGyVQv/32m0pHbL755ht59OiRfP3115I3b145cOCAahsHcW2zUKFC8QqMMMsI8g1/oU6dOmrIFqM333xTDSuDkrCpU6eqdnS7du0y5Q0CMeQn9hdpx//IK6QXkE4EVqtXr1btH5EGvI/qXJ35MUN6jW06x4wZo4agQWkf2nuWKlVKDaGCYVkAy+IzOpwf2AaGhQEc/+HDh6t90NsqGn388ccqsMSA1Dh2KGH99NNPpWjRomoxBPHIB3MI6H19fU0vY8gZtM3UIQ+wH0k9biGRI3NC3xBHzgAisk9ov4ZAyryjCBGRo2PJHxEREZEDYfBHRHYJ7dOMQ5cQEVEkVvsSERERORCW/BERERE5EAZ/RERERA6EwR8RERGRA2HwR0RERORAGPwRERERORAGf0REREQOhMEfERERkQNh8EdERETkQBj8ERERETmQ/wPvKLZaDtdzDAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 650x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "power at n_t: 0.8001\n"
     ]
    }
   ],
   "source": [
    "ns = np.linspace(200, 2 * P['n_t'], 120)\n",
    "pw = [power_at(n) for n in ns]\n",
    "fig, ax = plt.subplots(figsize=(6.5, 3.6))\n",
    "ax.plot(ns, pw, color='#a3471c')\n",
    "ax.axhline(0.8, ls='--', color='grey', lw=0.8)\n",
    "ax.axvline(P['n_t'], ls=':', color='#1f5f63')\n",
    "ax.set_xlabel('Enrolled per treatment arm'); ax.set_ylabel('Power, one contrast')\n",
    "ax.set_ylim(0, 1); ax.set_title('Power vs sample size (p0 = 0.30, MDE = 5 pp, Holm)')\n",
    "plt.tight_layout(); plt.show()\n",
    "print('power at n_t:', round(power_at(P['n_t']), 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63a84f91",
   "metadata": {},
   "source": [
    "## 2. Power by simulation of the pre-specified analysis model\n",
    "\n",
    "The analytic formula is an approximation. The test that matters is whether the *analysis we will\n",
    "actually run* rejects as often as planned. For each replication we:\n",
    "\n",
    "1. create $N$ participants in four strata (sex × baseline activity), randomised in equal shares to the\n",
    "   four arms within each stratum;\n",
    "2. drop 5 % at random (closed accounts, ported numbers: outcome missing);\n",
    "3. draw the outcome; 4. estimate\n",
    "\n",
    "$$Y_i = \\alpha + \\sum_{k=1}^{3} \\tau_k T_{ki} + \\sum_{s=2}^{4} \\gamma_s S_{si} + \\varepsilon_i$$\n",
    "\n",
    "by OLS with HC2 standard errors; 5. compute two-sided p-values for $\\tau_1,\\tau_2,\\tau_3$ and apply\n",
    "Holm (and, for comparison, Bonferroni) at family-wise $\\alpha = 0.05$.\n",
    "\n",
    "The OLS/HC2 step is written in numpy for speed; the next cell checks it against statsmodels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "70e3208f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:47.131686Z",
     "iopub.status.busy": "2026-11-03T14:49:47.131438Z",
     "iopub.status.idle": "2026-11-03T14:49:47.548758Z",
     "shell.execute_reply": "2026-11-03T14:49:47.546820Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max |coef diff|: 3.608224830031759e-16  max |SE diff|: 1.7277845820728999e-15\n"
     ]
    }
   ],
   "source": [
    "def ols_hc2(y, X):\n",
    "    XtX_inv = np.linalg.inv(X.T @ X)\n",
    "    b = XtX_inv @ X.T @ y\n",
    "    e = y - X @ b\n",
    "    h = np.einsum('ij,jk,ik->i', X, XtX_inv, X)\n",
    "    meat = (X * (e ** 2 / (1 - h))[:, None]).T @ X\n",
    "    V = XtX_inv @ meat @ XtX_inv\n",
    "    return b, np.sqrt(np.diag(V))\n",
    "\n",
    "def holm_reject(p, alpha=0.05):\n",
    "    order = np.argsort(p); K = len(p); rej = np.zeros(K, bool)\n",
    "    for step, idx in enumerate(order):\n",
    "        if p[idx] <= alpha / (K - step):\n",
    "            rej[idx] = True\n",
    "        else:\n",
    "            break\n",
    "    return rej\n",
    "\n",
    "def simulate_once(rng, n_t, effects, p_strata=(0.30, 0.30, 0.30, 0.30), attrition=0.05):\n",
    "    arms = len(effects) + 1\n",
    "    n = n_t * arms\n",
    "    strata = rng.integers(0, 4, n)\n",
    "    arm = np.empty(n, int)\n",
    "    for s in range(4):                       # complete randomisation within stratum\n",
    "        idx = np.flatnonzero(strata == s)\n",
    "        lab = np.resize(np.arange(arms), idx.size)\n",
    "        arm[idx] = rng.permutation(lab)\n",
    "    keep = rng.random(n) >= attrition\n",
    "    strata, arm = strata[keep], arm[keep]\n",
    "    p = np.asarray(p_strata)[strata] + np.r_[0.0, effects][arm]\n",
    "    y = (rng.random(arm.size) < p).astype(float)\n",
    "    X = np.column_stack([np.ones(arm.size)] + [(arm == k).astype(float) for k in range(1, arms)]\n",
    "                        + [(strata == s).astype(float) for s in range(1, 4)])\n",
    "    b, se = ols_hc2(y, X)\n",
    "    t = b[1:arms] / se[1:arms]\n",
    "    pval = 2 * norm.sf(np.abs(t))\n",
    "    return pval, (y, X)\n",
    "\n",
    "# cross-check one replication against statsmodels\n",
    "import statsmodels.api as sm\n",
    "rng = np.random.default_rng(1)\n",
    "pval, (y, X) = simulate_once(rng, P['n_t'], [0.05, 0.05, 0.05])\n",
    "fit = sm.OLS(y, X).fit(cov_type='HC2')\n",
    "b, se = ols_hc2(y, X)\n",
    "print('max |coef diff|:', np.abs(fit.params - b).max(), ' max |SE diff|:', np.abs(fit.bse - se).max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fb0dec72",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:49:47.552554Z",
     "iopub.status.busy": "2026-11-03T14:49:47.552255Z",
     "iopub.status.idle": "2026-11-03T14:50:01.450821Z",
     "shell.execute_reply": "2026-11-03T14:50:01.449588Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A. all three reminders +5 pp\n",
      "   rejection rate per contrast, Holm      : [0.862 0.855 0.861]\n",
      "   rejection rate per contrast, Bonferroni: [0.808 0.792 0.788]\n",
      "   any false rejection (FWER), Holm       : 0.0\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "B. only goal framing +5 pp\n",
      "   rejection rate per contrast, Holm      : [0.029 0.793 0.024]\n",
      "   rejection rate per contrast, Bonferroni: [0.024 0.792 0.016]\n",
      "   any false rejection (FWER), Holm       : 0.047\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C. global null (no effects)\n",
      "   rejection rate per contrast, Holm      : [0.024 0.015 0.018]\n",
      "   rejection rate per contrast, Bonferroni: [0.024 0.014 0.016]\n",
      "   any false rejection (FWER), Holm       : 0.046\n",
      "analytic power per contrast: 0.8  Monte Carlo SE ~ 0.009\n"
     ]
    }
   ],
   "source": [
    "def power_sim(effects, reps=2000, seed=20261103, **kw):\n",
    "    rng = np.random.default_rng(seed)\n",
    "    holm = np.zeros((reps, 3), bool); bonf = np.zeros((reps, 3), bool)\n",
    "    for i in range(reps):\n",
    "        pval, _ = simulate_once(rng, P['n_t'], effects, **kw)\n",
    "        holm[i] = holm_reject(pval); bonf[i] = pval <= 0.05 / 3\n",
    "    return holm, bonf\n",
    "\n",
    "scenarios = {\n",
    "    'A. all three reminders +5 pp': [0.05, 0.05, 0.05],\n",
    "    'B. only goal framing +5 pp':   [0.00, 0.05, 0.00],\n",
    "    'C. global null (no effects)':  [0.00, 0.00, 0.00],\n",
    "}\n",
    "res = {}\n",
    "for name, eff in scenarios.items():\n",
    "    holm, bonf = power_sim(eff)\n",
    "    res[name] = (holm, bonf)\n",
    "    print(name)\n",
    "    print('   rejection rate per contrast, Holm      :', holm.mean(0).round(3))\n",
    "    print('   rejection rate per contrast, Bonferroni:', bonf.mean(0).round(3))\n",
    "    print('   any false rejection (FWER), Holm       :',\n",
    "          round((holm & (np.array(eff) == 0)).any(1).mean(), 3))\n",
    "print('analytic power per contrast:', round(power_at(P['n_t']), 3), ' Monte Carlo SE ~', round(sqrt(0.8 * 0.2 / 2000), 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45e64ca3",
   "metadata": {},
   "source": [
    "**Reading the simulation.**\n",
    "\n",
    "* Scenario B is the planning scenario (one framing works, the others do not). Simulated power for\n",
    "  that contrast should sit at about the analytic 80 %, within Monte Carlo error (about ±2 points at\n",
    "  two standard errors). Holm and Bonferroni give almost the same power here: with only one real\n",
    "  effect, Holm's first (and binding) threshold is the Bonferroni threshold.\n",
    "* Scenario A shows that Holm's power per contrast rises when several arms work: after the first\n",
    "  rejection the remaining thresholds relax to $\\alpha/2$ and then $\\alpha$.\n",
    "* Scenario C checks the error rate we promise: the probability of declaring *any* reminder effective\n",
    "  when none is should be at or below 5 %.\n",
    "\n",
    "## 3. Covariate gain: why the analytic number is conservative\n",
    "\n",
    "In the scenarios above the strata are not prognostic, so stratum fixed effects add nothing. In a real\n",
    "wallet base, people already active at baseline are far more likely to be regular savers. If the\n",
    "regular-saver rate is, say, 10 % among the baseline-dormant and 50 % among the baseline-active (still\n",
    "30 % on average), the same design and the same model have more power. We do not plan on this gain;\n",
    "it is a margin of safety against an optimistic effect size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1e6c8860",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:50:01.453139Z",
     "iopub.status.busy": "2026-11-03T14:50:01.452915Z",
     "iopub.status.idle": "2026-11-03T14:50:06.042867Z",
     "shell.execute_reply": "2026-11-03T14:50:06.042308Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scenario B with prognostic strata, Holm power for goal-framing contrast: 0.878\n"
     ]
    }
   ],
   "source": [
    "holm, bonf = power_sim([0.00, 0.05, 0.00], p_strata=(0.10, 0.10, 0.50, 0.50))\n",
    "print('Scenario B with prognostic strata, Holm power for goal-framing contrast:', holm.mean(0)[1].round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4861cbac",
   "metadata": {},
   "source": [
    "## 4. Minimum detectable effect for a fixed budget\n",
    "\n",
    "If the partner can only include, say, 1,200 people per arm, what is the smallest effect the trial can\n",
    "reliably detect? Solve $\\text{power}(\\delta) = 0.8$ for $\\delta$ by bisection (the demo does the same)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "51e9d1df",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:50:06.044901Z",
     "iopub.status.busy": "2026-11-03T14:50:06.044669Z",
     "iopub.status.idle": "2026-11-03T14:50:06.086587Z",
     "shell.execute_reply": "2026-11-03T14:50:06.085618Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_t =   800:  MDE = 7.83 percentage points\n",
      "n_t =  1200:  MDE = 6.37 percentage points\n",
      "n_t =  1929:  MDE = 5.00 percentage points\n",
      "n_t =  2500:  MDE = 4.38 percentage points\n"
     ]
    }
   ],
   "source": [
    "def mde_at(n_t, target=0.80, **kw):\n",
    "    lo, hi = 1e-6, 1 - kw.get('p0', 0.30) - 1e-6\n",
    "    for _ in range(80):\n",
    "        mid = (lo + hi) / 2\n",
    "        if power_at(n_t, delta=mid, **kw) < target:\n",
    "            lo = mid\n",
    "        else:\n",
    "            hi = mid\n",
    "    return hi\n",
    "\n",
    "for n in [800, 1200, P['n_t'], 2500]:\n",
    "    print(f'n_t = {n:>5}:  MDE = {100 * mde_at(n):.2f} percentage points')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18e2d835",
   "metadata": {},
   "source": [
    "## 5. What the trial is *not* powered for: differences between framings\n",
    "\n",
    "Comparing two active framings (goal vs loss) is a secondary question. If both raise the regular-saver\n",
    "rate and differ by 2 percentage points (35 % vs 37 %), the per-arm sample needed for 80 % power at the\n",
    "same planning α is far larger than the trial. The pre-analysis plan therefore reports framing\n",
    "contrasts as secondary and does not let the programme decision rest on them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "00b208cf",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-11-03T14:50:06.088343Z",
     "iopub.status.busy": "2026-11-03T14:50:06.088158Z",
     "iopub.status.idle": "2026-11-03T14:50:06.092911Z",
     "shell.execute_reply": "2026-11-03T14:50:06.091937Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "per arm for a 2 pp framing difference: 12690  vs planned per arm: 1929  (ratio 6.6x)\n"
     ]
    }
   ],
   "source": [
    "q = plan(p0=0.35, delta=0.02)\n",
    "print('per arm for a 2 pp framing difference:', q['n_t'], ' vs planned per arm:', P['n_t'],\n",
    "      f' (ratio {q[\"n_t\"] / P[\"n_t\"]:.1f}x)')"
   ]
  }
 ],
 "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.13.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}